MAGIX: Model-based Analysis of Genomic Intervals with eXponential enrichment.

MAGIX is a generative model that explicitly connects the enrichment of TF-bound genomic intervals to the fragment counts observed across GHT-SELEX cycles. MAGIX, models how TF-bound intervals progressively occupy a higher proportion of selected fragments pool in each cycle relative to genomic background. These fragment proportions, in turn, are treated as latent variables in the model that, together with a sample-specific library size factor, determine the number of observed reads through a Poisson process.

Requirements

• Unix-compatible OS.
• R version 3.0.1 or later.
• R libraries: ComplexHeatmap, circlize, ggplot2, tictoc, Matrix, uwot, stringr, data.table, PRROC, sads, Rcpp, and optparse.
• Python version 3.6.0 or later.

Installation

git clone https://github.com/csglab/MAGIX.git

After cloning, you can add the line export PATH=${cloning_directory}/MAGIX:$PATH to your .bashrc file.

To create a conda environment with all the requirements:

# Create environment named MAGIX_env and install R
conda create --name MAGIX_env r-base r-essentials
conda activate MAGIX_env
# Check that you are using the R installation from conda
which R Rscript
# Install R dependencies
Rscript ./src/_R/install_R_lbraries.R

Aternatively, we have created an Apptainer container image (sif file) for MAGIX. To run commands with this image:

apptainer exec MAGIX.sif MAGIX --help

You can follow the steps here to install Apptainer.

Usage

Step 1: Analyzes every sample for the Target TF separately. Creates a coefficients correlation heatmap, Q-Q plots (for all sample permutations), and estimates library size factors.

Step 2: Combines the samples and fits coefficients for the target TF itself. You can limit the analysis to samples that are selected with the results from step 1 (e.g. *_step_1_coefs_heatmap.pdf) using --step2_batches (default = all).

Step 3: Performs a Likelihood-ratio test with the complete model from Step 2 and a reduced model without the target TF coefficient. It performs the LRT with the entire dataset by default, can be limited to a subset of the peaks with --ltr_sample_size.

Demo

Demo scripts and corresponding dataset can be found in ./demo_scripts and ./data/CTCF_demo/IN. Demo run time for Step 1 and 2 is ~ 5 min each and ~2 hours for Step 3.

For example:

MAGIX --outdir ${out_directory} \
--step 1 \
--count_matrix ${count_matrix} \
--design ${design_file} \
--TF CTCF \
--TF_label CTCF_FL \
--account_covariance TRUE \
--test_depletion FALSE

You can download the demo output here.

Input files

• --count_matrix COUNT_MATRIX Count matrix each column name corresponds to the Experiment_ID column in the design file (--design).
• --design DESIGN Design matrix, must have columns for Batch, Cycle, Experiment_ID, and Target. E.g. ./data/CTCF_demo/IN/CTCF_design_matrix_per_TF.txt.
• --step2_batches For step 2 and 3, sample names to include in the analysis (default = all).

Preprocessing and pipeline

The scripts used to cut adapter, align, obtain the count matrix, run MAGIX genome wide (to estimate the library sizes), and run MAGIX on regions with a minimum signal are provided here. This corresponds to the following part of the the methods section:

For each TF, an aggregate BAM file is created across replicates and cycles. For each genomic region with continuous non-zero coverage in this aggregate BAM file, the maximum coverage is calculated. Regions with a maximum coverage smaller than a threshold are discarded, with the threshold being set to the total fragment count in the aggregate BAM file divided by 2×106. Then, per region, the coordinate with the highest coverage is defined as the summit. The candidate peaks are defined as the 200 bp regions centered on the summits. Once the candidate peaks are identified, their count profiles are calculated from the original (unaggregated) BAM files (per cycle and replicate) and used as input to MAGIX to calculate enrichment coefficients, using prefixed library sizes that are estimated by fitting a similar MAGIX model to count profiles of 200-bp non-overlapping genomic intervals (a total of ~13M bins). For each candidate peak, we also calculate a P-value, representing the statistical significance of the enrichment coefficient (null hypothesis is that the enrichment coefficient is zero). To do so, we obtain maximum likelihood estimate of the model coefficients and perform a likelihood ratio test (LRT) against a reduced model in which the enrichment coefficient is restricted to zero. FDRs are calculated using the Benjamini & Hochberg correction.

Citation

Jolma, A., Hernandez-Corchado, A., Yang, A. W. H., Fathi, A., Laverty, K. U., Brechalov, A., Razavi, R., Albu, M., Zheng, H., Kulakovskiy, I. V., Najafabadi, H. S., & Hughes, T. R. (2026). GHT-SELEX demonstrates unexpectedly high intrinsic sequence specificity and complex DNA binding of many human transcription factors. Nature Methods, 1–11. https://doi.org/10.1038/s41592-026-03177-9

Heatmaps are generated with the ComplexHeatmap and circlize packages. If you use them in published research, please cite:

• Gu, Z. Complex Heatmap Visualization. iMeta 2022. or
• Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016. and
• Gu, Z. circlize implements and enhances circular visualization in R. Bioinformatics 2014.