Uniform processing and analysis of IGVF massively parallel reporter assay data with MPRAsnakeflow

MPRA design

Arbitrary sequences can in principle be tested via an MPRA, constrained only by technical limitations like construct or vector length. Typically, sequences are generated by DNA oligo synthesis, which currently is still most economical for sequence lengths of ∼300 nucleotides (nt) and covers a majority of open chromatin peak lengths. Longer sequences can be derived either from combining shorter sequences or from selections or amplifications of genomic DNA fragments. When barcodes are not directly synthesized with each sequence but added randomly (e.g., by amplification with overhanging barcode primers), the resulting MPRA library must be sequenced to determine barcode-to-test sequence associations. Paired-end short-read sequencing is usually sufficient to cover the entire test sequence with a small overlap (we recommend >10 nt for high-complexity sequences). For example, 2 × 150 nt paired-end sequencing can cover a 200 nt test oligo + 20 nt barcode + 30 nt spacer with 50 nt overlap. Alternatively, long-read sequencing can be employed. However, we note that various MPRA designs have different requirements on the precision of these sequence read-outs to obtain an unambiguous association, specifically if sequence variants need to be distinguished.

Some MPRA experiments test synthetic or engineered sequences that possess particular activity properties, such as tissue-specific enhancer activity or optimal spacing between two elements (de Almeida et al. 2022; Georgakopoulos-Soares et al. 2023; Gosai et al. 2024). However, most MPRAs are used to test cCREs derived from a genomic or inferred (ancestral) genomic sequence, and identified through epigenetics maps, motif scans, or sequence models. For longer fragments, tiling might be necessary (Kosicki et al. 2025). These choices will impact the average element activity in the assay and thereby experimental sensitivity. Within IGVF activities, variants are one complementary focus of the consortium. Again, variants may be synthetic and prioritized using sequence models or selected from various sources like population data, ancestral alleles, or prior functional mapping from genome-wide association studies (GWAS) or expression quantitative trait loci (eQTL) studies. For variants, the reference and alternative sequence is designed. Usually one variant is centered within the oligo, but alternative positioning as well as the introduction of variant combinations is possible. We note that a sufficient reference activity is generally important to detect variant effects. Therefore, shifting the sequence toward known transcription factor binding sites (TFBSs) or using the cCRE is beneficial.

In addition to test sequences (e.g., variants or cCREs), control sequences are required to validate experiment performance. By allowing the quantification of baseline activity, they help to identify significantly activating or repressing sequences and establish that the activity distribution of tested sequences exceeds that of negative controls. The most common approach is to use (di-nucleotide) shuffled sequences as negative controls to estimate background activity. Positive element or variant controls are more challenging to design. Ideally there exist previously tested sequences with high activity in the same cell type, or variants with known properties in the same system. If such sequences are not available, sequences with activity across many cell types or from a related cell type should be used. Generally, only a small number of positive controls (around hundred) might be required, whereas a larger set of several hundreds of negative controls is necessary for later statistical testing. However, if previously validated positive controls are unavailable, more putative positive controls are recommended.

We recommend avoiding the use of sequences with certain undesirable properties. Sequences containing long homopolymers or specific restriction sites (e.g., EcoRI or SbfI) used in MPRA library construction should be removed. Restriction sites may also be eliminated by sequence modification. Simple repeats and overlap with transposable elements may cause unwanted strong activity or construct silencing. The inclusion of sequences with binding sites for factors impacting chromatin structure (e.g., CTCF binding sites) should be excluded, because short MPRA constructs cannot model chromatin looping structures. Finally, sequences containing additional transcriptional start sites may interfere with the minimal promoter, resulting in inconsistent background activity of the core element.

To facilitate this process, we developed a workflow called MPRAOligoDesign, which can generate a design file from sequences, variants, genomic regions, and paired variants with genomic regions. It extends regions to the desired sequence length, tiles longer regions, designs reference and alternative sequences for variants, and can be configured to perform recommended property checks, removing sequences that do not pass the filters.

Standardized file formats

Within the past 10 years, more than 100 MPRA experiments have been published. Various efforts have been made to integrate different experiments into common databases, like MaveDB, MPRAVarDB, and MPRAbase (Jin et al. 2024; Rubin et al. 2025; Zhao et al. 2025). However, because there are no defined standards for metadata, such as the genomic region tested or sequence design, setting up similar data resources becomes very challenging. One major resource of IGVF will be a catalog of element and variant effects. With various different MPRAs being conducted, they must be integrated into a single catalog as well as a data portal for reuse by other projects. We see a tremendous need to harmonize metadata regarding the design and various stages of output files for MPRA data. We identified four categories for which standardized files are necessary—design, counts, effects, and genomic effects—and provide a detailed documentation of the file formats in Supplemental Note S1.

For “design,” we defined a reporter sequence design file that contains the main metadata of designed oligos such as identifier, sequence, sequence origin, and class and provides key information about the purpose of the assay. Despite the “designed oligo” terminology, such design files should also be created for MPRA experiments not derived from oligo synthesis, and equivalent information should be collected. In the design file, each oligo has a unique name or ID and the nucleotide sequence. It further includes the category (variant, element, synthetic, scrambled), an assigned class (test, variant positive control, variant negative control, element active control, element inactive control), and additional information. Chromosome and start and end positions, along with the reference genome version, must be specified if applicable. For variants, the position within the tested sequence, a unique ID using the canonical SPDI notation (Holmes et al. 2020), and whether the reference or alternative allele was tested within the sequence must also be included. This file enables reanalysis of the data as well as the generation of element and/or variant effects from the data.

“Count” files report the number of observed reads (i.e., integer counts) counted per barcode or aggregated to elements per modality (i.e., DNA or RNA). Multiple counts and modalities can be described in the same file. Barcode-level files contain the barcode sequence and the associated oligo identifier. Element-aggregated files describe the replicate ID, oligo ID, aggregated and normalized DNA and RNA counts, the number of aggregated barcodes, and an inferred log-fold change (normalized RNA/DNA ratio). For IGVF, these are generated via the standardized pipeline MPRAsnakeflow (see below). However, existing published pipelines, like MPRAflow (Gordon et al. 2020), also generate similar types of data and could be converted to our standardized format.

“Effect” files report the results of differential activity analysis. This includes normalized counts, activity estimates (i.e., log-fold changes), and significance estimates (i.e., P- and q-values). This can be performed using tools like MPRAnalyze, BCalm, or mpralm (Ashuach et al. 2019; Myint et al. 2019; Keukeleire et al. 2025), with the count data and information from the sequence design file as input (e.g., to map oligos to reference or alternative alleles for a variant). Because differential activity analysis differs between elements and variants (e.g., statistical testing for elements always requires a control set of oligos), we separate them into element and variant effects files. In addition to the values listed above, element files contain aggregated modality counts. Variant files contain those for both the alternative and reference allele, a posterior probability of the regulatory effect, and a 95%-confidence interval of the effect.

Finally, “genomic effects” are those that can be mapped to a genomic location, and we have created standardized BED files for variants and elements. These contain information from effect files with genomic coordinates for visualization of MPRA results in genome browsers such as UCSC or Integrative Genomics Viewer (IGV) (Robinson et al. 2011; Perez et al. 2025). Tools for MPRA data from IGVF will support these standards. In our MPRAlib library (see below) a command set called validate-file is implemented to compare actual files with their schemas to ensure files are formatted correctly.

A unified MPRA pipeline: MPRAsnakeflow

MPRAsnakeflow is a further evolution of MPRAflow (Gordon et al. 2020), which established basic functionality for the analysis of MPRA data. We focus on six topics to enhance and extend its functionality: (1) speed and memory improvements to support large complex MPRA assays; (2) interoperability by containerization of software dependencies; (3) standardized outputs, quality reports and additional statistics to help uncover potential issues with the MPRA experiment; (4) support for multiple assays (incl. raw read structures); (5) optimized barcode assignment to support designs with small edit distances, like multiple variants for the same cCRE; and (6) new features like strand awareness, outlier removal, and support of multiple tools for statistical quantification. For a detailed comparison of MPRAflow with MPRAsnakeflow, see Supplemental Table S1. MPRAsnakeflow is implemented in the workflow language Snakemake (Mölder et al. 2025). It is organized into two subworkflows: the assignment workflow and the experiment workflow (see Fig. 1).

Figure 1.

Overview of an MPRA experiment with UMIs (top) and the MPRAsnakeflow pipeline for analyzing MPRA experiments (bottom). MPRAsnakeflow consists of an assignment workflow (left), which maps barcodes to designed oligos, and an experiment workflow (right), which counts barcodes and computes oligo activity. Inputs are shown in yellow, main processing steps in blue, and outputs in purple. Notably, MPRAsnakeflow supports different input read structures, including separate reads for oligo and barcode (A), as well as combined structures in which the barcode and oligo are present within read pairs (B). Reads can also be single-end reads.

Open in new tabDownload PowerPointLink to figure

The assignment workflow aims to associate barcodes with designed oligos. It supports multiple read structures typical for certain assay types (Gordon et al. 2020; Gosai et al. 2024). For example, the barcode can either be sequenced in a separate index read or positioned at the beginning of the first paired-end read, separated from the start of the oligo by a linker. In brief, the assignment workflow first verifies that the design file (in multi-FASTA format) is appropriate, ensuring unique headers and no duplicate sequences. Next, paired-end reads of the sequenced oligos are merged using NGmerge (default) or fastq-join (Aronesty 2013; Gaspar 2018), with a minimum recommended overlap of 11 nt. The merged reads are then mapped to the design file using a read mapper (BBMap [default], BWA-MEM [Li and Durbin 2009; Bushnell 2014], or exact string matching). For long-read data, we implemented the Pacific Biosciences (PacBio) minimap wrapper pbmm2 (Li 2018). Barcodes are assigned based on the best unique sequence match. Finally, barcodes are aggregated and filtered based on a minimum number of required observations per oligo (default: three), and a majority of barcode occurrences must map to the same oligo (default: 0.75) to minimize ambiguity while allowing barcode collisions or best alignment errors. This workflow can be skipped if barcode associations are already known (e.g., by design or with custom assignment scripts).

The experiment workflow processes the sequenced barcodes extracted from DNA and RNA samples of the MPRA experiment in a specific tissue or cell type. It counts the observations of each barcode, merges DNA and RNA counts, and assigns barcodes to oligos using an assignment file, which can be generated by the assignment workflow. An optional step removes outlier barcodes. Barcode counts per oligo are then aggregated, and activity is calculated as the log₂ ratio of normalized (counts per million reads) RNA over DNA counts for each oligo. When multiple replicates are available, the workflow computes correlations of DNA counts, RNA counts, and log₂ ratios to assess data quality. Additionally, the experiment workflow supports the use of UMIs to disambiguate unique transcripts from PCR duplicates to reduce PCR artifacts.

Both subworkflows produce a variety of statistics and plots to help users inspect the data and identify potential errors. A unique feature of MPRAsnakeflow is the quality control report generated for both workflows, which summarizes key metrics and plots along with additional explanations. The quality report is provided as a single structured and image-embedded HTML file, enabling quick and efficient data inspection and making it easy to share results with collaborators. We have comprehensively documented MPRAsnakeflow; full documentation is available at https://mprasnakeflow.readthedocs.io. The documentation provides installation instructions, detailed descriptions of the available configuration options and workflow components, and a step-by-step tutorial. In addition, we present analyses of more than 10 example data sets, including recent, independently generated MPRA studies available from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) that were not part of the IGVF consortium, representing a broad range of library designs. Together, these examples illustrate the flexibility of the workflow and provide practical guidance for users seeking to adapt it to their own assays.

Mapping strategies for barcode associations

To mitigate the impact of sequencing errors, sequenced oligos are typically mapped (Alser et al. 2021) to the originally designed sequences of the MPRA experiment. Approaches relying solely on exact string matching, as implemented in some other pipelines like esMPRA (Li et al. 2025), fail to account for sequencing errors and consequently result in a reduced number of assigned reads and therefore less observed barcodes. The MPRAflow tool (Gordon et al. 2020) uses solely BWA-MEM and does not allow different mapping configurations without changing the software. Our observations indicate that the choice of mapping strategy can substantially influence downstream analyses, depending on the library design. For this reason, MPRAsnakeflow offers multiple configurable strategies: exact matching (requiring the full read sequence matching the designed oligo, in either the forward or reverse-complement orientation), BWA-MEM, or BBMap. We observed that BWA-MEM reports MAPQ 0 values despite distinguishing oligos with minimal sequence differences in alignment (as apparent from alignment score differences reported in the output). This causes barcodes to be lost from associations in variant-rich MPRA libraries. To address this, we implemented an additional filtering step for BWA alignments, allowing users to still consider MAPQ 0 results based on the sequence similarity, number of mismatches, and alignment length.

Using different mapping strategies, we report in Figure 2 the assignable oligos, identified genomic regions without alternative alleles designed, the assigned variants (for which both reference and alternative oligo must be present), as well as the number of barcodes per oligo. In short, BBMap was able to assign more oligos on average (35,249) compared with BWA-MEM (33,736) and exact matching (33,471) across all MPRA sets used. The same holds for assigned regions and especially for variants. We see a drastic decrease in reported variants using the BWA mapping strategy especially for the 80K-neurons data set, confirming the outlined issues with MAPQ values when two oligos are in close edit distance (e.g., 1 nt for SNVs).

Figure 2.

Comparison of four different mapping strategies for the association step for three different MPRA libraries (20K-HepG2_P, 80K-neurons, 8K-neurons). From left to right, we present the number of assigned oligos, the number of assigned elements (candidate cis-regulatory elements, reference sequence only), and the number of retrieved variants. A variant is only considered if both the reference and the alternative allele oligos are assigned. The rightmost bars show the number of barcodes per oligo (mean across replicates). The best result for each MPRA library is highlighted in bold in each panel.

Open in new tabDownload PowerPointLink to figure

We further inspected the design on 80K-neurons because of the difference in variants for BWA-MEM and noticed that most reference alleles had at least three alternative alleles each (62%). Within this subset of oligos, the proportion of assignable variants dropped to <40% for BWA-MEM (see Supplemental Fig. S1). This shows that small edit distances play a crucial role in retaining variants and should be carefully considered for elements with many designed variants or for saturation mutagenesis (Kircher et al. 2019).

To compare the impact on downstream analyses, we utilized BCalm (Keukeleire et al. 2025) to determine significant variant effects for all four mapping strategies using two different MPRA experiments (8K-neurons [Koesterich et al. 2023], 80K-neurons). Supplemental Figure S2 shows the overlap in all variants and significant variants. BWA-MEM detected the lowest number of variants, followed by the exact matching, whereas most variants were detected by either the additional filtering applied to BWA output (8K-neurons) (Supplemental Fig. S2A) or BBMap (80K-neurons) (Supplemental Fig. S2B).

Variants detected by the default BWA-MEM were detected by all other mapping strategies, and only 0.1% (n = 2, 8K-neurons) or 0.2% (n = 96, 80K-neurons) of the variants were only in the exact matching strategy for the 8K-neurons and 80K-neurons data set, respectively. Overall, the concordance in identified variants is very high; for example, the proportion of variants detected across at least three out of four methods is 97.6% (n = 3406, 8K-neurons) and 94.7% (n = 37,956, 80K-neurons). When looking into the significant variants, concordance between the mapping strategies drops, with only 48.7% (n = 416, 8K-neurons) and 59.2% (n = 1011, 80K-neurons) of variants being significant for at least three out of four mapping strategies. We observe larger numbers of variants only significant for one of the mapping strategies. This occurs for 14.4% (n = 123, 8K-neurons) of variants (Supplemental Fig. S2C) and 24.3% (n = 416, 80K-neurons) of variants (Supplemental Fig. S2D), respectively. Notably for 8K-neurons, BWA-MEM with additional filtering (n = 721) and BBMap (n = 710) approaches yielded sets that are more than 1.5-fold larger than those of the other strategies (BWA n = 456, exact matching n = 420). When restricting to the two sets with the highest number of significant variants (BWA-MEM and BBMap), the concordance on the significant variants increased to >90%, which suggests that the higher number of barcodes assigned leads to a more stable variant effect quantification. Apart from the different number of mapped oligos, we could not identify large effects on replicate correlation (n = 3) across mapping strategies (see Supplemental Fig. S3).

In summary, the choice of the mapping strategy for the association significantly influences variant analysis outcomes on the count level, particularly for complex and large libraries. Exact matching approaches fail to account for sequencing errors, leading to fewer barcode and oligo assignments, whereas error-tolerant methods like BWA struggle with confident distinction of related sequences, necessitating custom scripts for mitigating this issue. MPRAsnakeflow implements BBMap as standard, because it showed an overall good performance (Fig. 2; Supplemental Fig. S1), but we recommend checking for other mapping strategies if the number of assigned barcodes or retained oligos is low.

Complexity and sequencing depth estimation

A common question is whether deeper sequencing of an MPRA library improves data quality, such as increasing correlation across replicates. To address this, the pipeline implements two key indicators. First, we use pairs of replicates to estimate overall library complexity using the Lincoln–Petersen method (Petersen 1896; Lincoln 1930). The difference between the observed number of barcodes and the Lincoln index provides an estimate of the potential improvement achievable with additional sequencing. The Lincoln index can also be used to identify replicates or modalities with decreased barcode complexity. Second, the pipeline includes a built-in downsampling option to assess whether quality metrics (e.g., replicate correlation) already reach saturation. If not, this indicates the potential benefit of further sequencing. Downsampling can be performed either on individual barcode counts or on barcodes within the assignment file, using either a fixed number or a specified proportion. Multiple downsampling configurations can be applied within a single experiment run using named configurations, ensuring that only the necessary parts of the workflow are rerun, which saves disk space and computational resources.

For the small library experiment, 8K-neurons (Koesterich et al. 2023), we observe a median of 1,444,480 assigned barcodes across replicates after assignment and RNA and DNA merging, requiring barcodes to have at least one count of DNA and RNA in the full data set. The median of pairwise replicate Lincoln index values is 1,523,572. Consequently, we are missing ∼5% of the barcodes from each replicate within this data set. We assume that this might be a common loss from transfections or transductions of MPRA libraries into cell populations.

The medium MPRA library (80K-neurons) (see Supplemental Table S2) presents a median of 5,459,247 assigned barcodes across replicates and a median Lincoln index of 6,243,618. Thus, we are missing ∼12% of the library in each replicate. The Lincoln index for the large 240K-HepG2 library (see Supplemental Table S2) is 20,408,916, and we observe a median of 13,038,694 assigned barcodes across replicates. This is ∼36% of missing barcodes.

Depending on the size of the designed library and the estimated complexity of the data, we see an increasing gap between missing barcodes in our libraries and a higher variability across replicates. For the medium and low complexity libraries, we see a good concordance between replicates, but the question arises whether more sequencing data could help for the higher complexity library. Therefore downsampling of RNA counts was performed for fixed proportions with MPRAsnakeflow (before assignment) and also for the final assigned count file (barcode reporter experiment) using MPRAlib. For all libraries (for 8K-neurons, see Supplemental Fig. S4A; for 80K-neurons, see Fig. 3A,B; for 240K-HepG2, see Supplemental Fig. S4B), we see a saturation in number of assigned barcodes as well as the observed oligos and the Pearson's correlation of the oligo activity (log₂-fold change). Both approaches show very similar trajectories and similar values for the observed number of barcodes. MPRAsnakeflow downsampling has a smoother curve and is more accurate but takes more time and computational resources, whereas using MPRAlib can be done within minutes, depending on the library size. Based on these results, 8K-neurons and 80K-neurons are already saturated, and new sequencing data will not substantially improve these metrics (Fig. 3; Supplemental Fig. S4A). For 240K-HepG2 (Supplemental Fig. S4B), we still see a slightly increasing line but with no improvement in the number of detected oligos and Pearson's correlation. Therefore, the benefits of resequencing will be neglectable.

Figure 3.

Downsampling analysis for 80K-neurons data set. RNA is downsampled in equal proportions to the overall number of counts. Complexity (Lincoln index), the number of barcodes, as well as retained oligos and Pearson's correlation of oligo activity across replicates (median values across replicates or replicate comparisons). (A) RNA count downsampling in MPRAsnakeflow. Sampling is done before the assignment on the raw counts. (B) Downsampling using the MPRAlib on the overall output (after the assignment).

Open in new tabDownload PowerPointLink to figure

A library for interactive MPRA count data analysis: MPRAlib

We developed a Python library, MPRAlib, to provide a simple and accessible framework for analyzing MPRA count data (Supplemental Fig. S5A). MPRAlib uses the AnnData object (Virshup et al. 2024) as its core data structure, with variables representing barcodes or oligos and observations corresponding to experimental replicates. The AnnData object is designed to support both barcode-level and oligo-level count data, enabling functionalities such as barcode outlier detection and data aggregation per oligo.

The main class, MPRAData, is an abstract base class that implements most core functions. Two subclasses, MPRABarcodeData and MPRAOligoData, inherit from MPRAData to support analyses at different data granularities. AnnData layers are used to store different modalities (e.g., RNA or DNA counts), normalized counts, and activity measures (such as log₂ ratios). The flexible AnnData structure also allows the addition of metadata and pairwise data to variables or observations, which is leveraged to store information such as the number of barcodes per oligo, masking arrays for filtered data, and pairwise correlations across replicates.

Built-in functions are provided for reading MPRAsnakeflow outputs (at both barcode and oligo levels). We implemented a dynamic approach for adjusting the minimum required number of barcodes per oligo, allowing users to change this threshold at runtime without losing access to the original raw data. This flexibility enables users to recompute metrics, such as correlations, on different thresholds, without reloading the data. At the barcode level, MPRAlib offers various options for detecting outliers, setting minimum required counts per modality, and subsampling counts or barcodes. The library includes a comprehensive suite of plotting functions, such as replicate correlations, DNA versus RNA scatter plots, and histograms of barcodes per oligo (for example plots, see Supplemental Fig. S5B–D). To facilitate integration with differential analysis tools such as BCalm, mpralm, and MPRAnalyze (Ashuach et al. 2019; Myint et al. 2019; Keukeleire et al. 2025), MPRAlib supports a standardized metadata format and includes functions for exporting data in formats compatible with these tools. Additionally, MPRAlib can import results from these tools and generate further outputs, such as BED file tracks for genome browsers.

MPRAlib is extensively documented (https://mpralib.readthedocs.io) and is accompanied by several tutorials using Python notebooks, demonstrating its use in various scenarios. For convenience, we also provide a command-line interface (CLI) implementing commonly used workflows, such as correlation analyses and plotting routines. The CLI allows users to validate input file formats against standardized schemas, which are defined using JSON Schema (see Supplemental Note S1). The library is written in Python (version ≥3.10) and is available for installation via PyPI and Bioconda (Grüning et al. 2018).

Parallel processing and performance benchmarking

We used MPRAsnakeflow to analyze a very large experiment (240K-HepG2) (see Supplemental Table S2) consisting of approximately 240,000 oligos testing for element-level activity (60,000 oligos) and variant-level activity (170,000 oligos). Supplemental Table S3 contains a detailed list of MPRAsnakeflow rules with number of jobs; maximum, minimum, and average run-time; and memory usage on the 240K-HepG2 data set. In summary for the assignment workflow, the three primary input FASTQ files (paired oligo reads and one barcode read) were split into 30 partitions (assignment_fastq_split) resulting in an association pipeline workflow consisting of 164 jobs, with four steps accounting for 150 of these jobs: assignment_attach_idx, 60 jobs; assignment_mapping_bbmap, 30 jobs; assignment_mapping_bbmap_getBCs, 30 jobs; and assignment_merge, 30 jobs. With a maximum number of parallel jobs set to 24, the workflow completed in under 6 h with a maximum of 28 GB RAM required for a single job. The majority of jobs for the assignment step required <10 GB of RAM and completed in under an hour.

The experiment workflow consisted of 103 jobs and required significantly more time because of only six jobs, specifically the counts_umi_create_BAM jobs for the DNA- and RNA-seq files (three of each). The DNA-seq FASTQ files took ∼12 h to process, whereas the larger RNA-seq FASTQ files required ∼34 h; however, these were all run in parallel. The remaining jobs required only four additional hours in total.

We run the same data set with MPRAflow to compare memory and runtime with MPRAsnakeflow. MPRAflow required a substantial amount of more memory (maximum up to 116 GB for the “map_element_barcodes” step) with a slightly longer runtime of 10 h for the rate limiting BAM creation step of the experiment workflow (see Supplemental Table S3).

Barcode discrepancy detection and outlier removal

We applied three outlier detection methods to identify problematic barcodes across nine MPRA data sets, comprising six lentiviral-based assays and three episomal-based assays (indicated by “_P” suffix). The three outlier detection methods included (1) global outliers based on z-score thresholds across all barcodes, (2) oligo-specific outliers based on deviation from within-oligo expression patterns, and (3) expression outliers with extreme log₂ expression ratios. Outlier detection patterns revealed noticeable differences between lentiviral and plasmid-based MPRA platforms (Fig. 4A). Global outlier detection showed platform-specific bias, with lentiviral assays exhibiting substantially higher rates (0.67%–1.56%) compared with plasmid-based assays (0.12%–0.47%). Oligo-specific outliers are quite consistent across all experiments, ranging from 0.85% to 1.84%. Notably, large expression outliers are very low in lenti-based assays (up to 0.01%) but had higher rates in plasmid-based assays, in particular 0.9% for 100K_GM12878_P (0.38% for 20K-HepG2_P, 0.03% for 60K-A549_P). This represents roughly fourfold to 90-fold increases over the lentiviral data sets. We also detected a difference in the barcode activity distribution (log₂ ratio) between lenti and episomal assays. Global outliers for episomal showed a high activity, similar to large expression outliers (see Supplemental Fig. S6A).

Figure 4.

Comparison of different MPRA experiments (lenti and episomal) using different outlier detection strategies. Episomal MPRAs have an “_P” at the end of their ID. (A) Outlier detection rates on different data sets using three different methods: global outliers (blue), oligo outliers (orange), and expression outliers (green). The percentage shown is the average among the replicates for global and oligo outliers. (B) Concordance of the global outlier detection method across replicates. All lenti-based as well as the 100K episomal assays have three replicates, and the 20K and 60K episomal libraries have six and five replicates, respectively. (C) Pearson's correlation of variant coefficients between log-fold change estimates from original data sets versus data sets with outliers removed using three methods: global outliers, oligo-specific outliers (oligo outliers), and large expression outliers. Number of variants per comparison are within the bars. Only data sets with more than 1000 tested variants are shown. (D) Jaccard indices measuring overlap of significant variants detected in original versus outlier-filtered data sets. Number of overlapping significant variants are written within the bars. Only data sets with more than 1000 tested variants are shown.

Open in new tabDownload PowerPointLink to figure

To assess the stability of barcode-level outlier identification across replicates, we calculated the percentage of barcodes that were consistently identified as outliers across all three biological replicates using the global outlier detection method. Outlier concordance varied across data sets, ranging from 5.6% to 78.6% (Fig. 4B). Although most lentiviral assays demonstrated poor to moderate concordance (5.6%–49.0%), plasmid-based assays exhibited substantially higher consistency rates (56.9%–78.6%). This higher concordance in plasmid assays may reflect more uniform transfection conditions and reduced technical noise from integration-site effects that can vary between lentiviral replicates. We also notice an increased outlier concordance observed in differentiated cell types (cardiomyocytes vs. cardioprogenitors: 43.8% vs. 7.8%, P < 10⁻⁴; neurons vs. WTC11: 49.6% vs. 29.7%, P < 10⁻⁴; for P-value estimation, see Methods) despite expected greater cellular heterogeneity in the former from incomplete differentiation. This may reflect more synchronized transcriptional states and reduced stochastic gene expression noise after differentiation. Additionally, differentiated cells may exhibit more stable chromatin landscapes and consistent regulatory network activity compared with progenitor cells, which are inherently more plastic and variable in their transcriptional responses.

To investigate potential sequence-based biases in outlier detection, we analyzed the GC content distribution of barcodes identified as large expression outliers compared with background barcodes (those with no outliers detected by any method) (Supplemental Fig. S6B). Lentiviral-based assays (240K-HEK293T, 240K-HepG2) showed no significant difference in GC content between large expression outliers and background barcodes, with both populations displaying similar median GC content and overlapping distributions. In contrast, all tested plasmid-based assays showed a significant elevation in GC content among the large expression outliers compared with the background controls. These results suggest that GC content bias in large expression outliers is not universal but potentially depends on specific experimental conditions, library characteristics, or cell type–specific interactions with GC-rich sequences. The absence of GC bias in lentiviral assays indicates that integration-based delivery may be less susceptible to sequence composition effects compared with episomal approaches.

To assess whether outlier detection is driven by specific problematic sequences or random technical noise, we analyzed barcode sequence overlap between data sets and the concordance of outlier calls for shared sequences (see Supplemental Fig. S7). The 15-nt barcode space used in most data sets is enormous (4¹⁵ ≈ 1 billion possible sequences), resulting in minimal sequence overlap between independent libraries. The data sets derived from the same underlying libraries showed substantial overlap: 240K-HEK293T versus 240K-HepG2 (16.2 million shared sequences), 12K-cardiop versus 12K-cardiom (1.9 million shared), and 80K-WTC11 versus 80K-neurons (5.4 million shared). Despite this substantial sequence overlap in related data sets, outlier concordance was low. For example, for 240K-HEK293T and 240K-HepG2, only three out of 16.2 million sequences were identified as large expression outliers in both data sets (Supplemental Fig. S7B). This discordance between sequence overlap and outlier concordance indicates that outliers are predominantly driven by data set–specific technical factors rather than inherent sequence properties.

To assess whether outlier removal improves the quality of MPRA variant effect analysis, we compared results from BCalm for original data sets and data sets with the outliers removed. We evaluated two key metrics across five data sets: Pearson's correlation of log-fold changes and Jaccard indices measuring the overlap of significant variants detected by each approach. The correlation analysis (Fig. 4C) showed consistently high correlations (>0.92) across data sets for all three outlier removal methods, indicating that outlier removal had minimal impact on overall effect size estimates. Similarly, the Jaccard index analysis (Fig. 4D) demonstrated high concordance in significant variant calls between original and filtered data sets, with 15 of 21 comparisons (71%) exhibiting Jaccard indices above 0.7. Notably, an inverse relationship was observed between the number of detected outlier barcodes (Fig. 4A) and the Jaccard index, indicating that more extensive outlier removal leads to greater divergence in the set of significant variants identified.

Mapping strategies

We utilize the BWA maximal exact matches (MEM) algorithm (Li and Durbin 2009) with an increased clipping penalty (default: 80) to discourage soft-clipping of read-ends by the aligner, effectively requiring a global sequence alignment. The mapper calculates a mapping quality score (MAPQ) based on the difference between the best and the next-best alignment. In MPRAsnakeflow, BWA-MEM alignments are then filtered using a configurable MAPQ threshold (default: one) to remove ambiguous mappings. Further, we implemented an additional BWA step using MAPQ 0 and filter alignments based on sequence similarity, number of mismatches, and alignment length.

Additionally, MPRAsnakeflow provides support for BBMap (Bushnell 2014), an alternative k-mer-based aligner. BBMap reports higher MAPQ values, reflecting increased alignment confidence even for variant-rich libraries. For BBMap, a more stringent default MAPQ threshold (30) is applied to remove ambiguous assignments.

To compare the performance of these mapping strategies, we analyzed three distinct MPRA data sets of different library sizes and design choices: 8K-neurons, 20K-HepG2_P, and 80K-neurons (for a description of the data sets, see Supplemental Note S2; Supplemental Table S2). For each data set, we applied the different mapping routines and compared the number of assignable oligos, the identified genomic regions, and the number of variant sequences detected (see Fig. 2). Only oligos associated with at least 10 unique barcodes were considered for these analyses. Furthermore, we used BCalm (Keukeleire et al. 2025) to obtain a variant effect quantification (log₂-fold change) and a Benjamini–Hochberg adjusted P-value per variant, omitting outlier barcodes; to identify significant variant effects (adjusted P-value < 0.1); and to assess the consistency of results across the tested mapping strategies.

Complexity and sequencing depth estimation

We randomly downsampled RNA counts from 5% to 95% of the original counts, in 5% increments. This analysis was performed on libraries of varying complexity: low (8K-neurons), medium (80K-neurons), and high (240K-HepG2) complexity (see Supplemental Table S2). After downsampling, we computed the Lincoln index, number of total assigned barcodes, and the total number of oligos retained with at least one barcode. To measure consistency of activity scores between replicates, we calculated Pearson's correlation across replicates at each downsampling level. We compared these results to those obtained using the downsampling strategy implemented in MPRAlib to assess whether the metrics are comparable.

Parallel processing and performance benchmarking

240K-HepG2 was processed with MPRAsnakeflow (v0.4.4) on a high-performance computing cluster. Individual job execution was managed using the slurm executor plugin for Snakemake: https://snakemake.github.io/snakemake-plugin-catalog/plugins/executor/slurm.html. The number of jobs per Snakemake rule, as well as time and memory usage, was recorded in Supplemental Table S3.

MPRAflow (v2.3.5; Nextlow version 20.01.0.5264) was run on the same high-performance computing cluster as MPRAsnakeflow. Individual job execution was managed by setting the “process.executor” to “slurm” in a configuration file. The reported metrics in Supplemental Table S3 were collated using the “elapsed” and “maxRSS” format options to extract run time and memory usage, respectively, from slurm accounting via “sacct.”

Barcode discrepancy detection and outlier removal

For the outlier analysis, seven IGVF data sets, one ENCODE data set, and one other MPRA data set were selected as follows: for lentiMPRAs, 240K-HEK293T, 240K-HepG2, 12K-cardiop, 12K-cardiom, 80K-WTC11, and 80K-neurons; for episomalMPRAs, 20K-HepG2_P, 60K-A549_P (Gosai et al. 2024), and 100K-GM12878_P (Abell et al. 2022). For more description on the data sets, see Supplemental Note S2 and Supplemental Table S2.

Raw count data were processed to calculate counts per million (CPM) for both DNA and RNA libraries in each replicate. We computed the sum of CPM values across replicates for DNA (dna_cpm_sum) and RNA (rna_cpm_sum), followed by the calculation of expression ratios as log₂(rna_cpm_sum/dna_cpm_sum). Quality filtering was applied to remove low-confidence measurements: Barcodes with zero RNA CPM or DNA CPM below the 5th percentile were excluded. Additionally, we retained only oligos represented by at least 20 barcodes across replicates to ensure sufficient statistical power. We implemented three complementary outlier detection approaches:

Outlier concordance comparison

Given the lack of a suitable parametric distribution describing the outlier concordance metric (proportion ∈ [0,1]), we estimated the P-value of the one-sided test that the larger measured value was greater than the smaller by using random permutations to determine null distributions. For a specific null (the smaller concordance proportion), we assumed the ability to detect all outlier barcodes (“null barcodes”) in all replicates that would yield the null proportion. We then assumed the additional observed outlier barcodes (n observed – n null) represent random noise and randomly sampled these from the total number of barcodes in the experiment less the number of null barcodes. We then calculated the cardinality of the set consisting of the intersection of all sampled replicates. The P-value was reported as the proportion of these sampled consistency metrics that exceeded the observed larger value, namely, the proportion of samples for which we reject the null.

Code availability

All custom software developed for this study is provided as Supplemental Code. MPRAsnakeflow (v0.7.0) (Supplemental Code S1) is available at GitHub (https://github.com/kircherlab/MPRAsnakeflow) and archived at Zenodo (https://doi.org/10.5281/zenodo.18163777). MPRAlib (v0.10.5; Supplemental Code S2) is available at GitHub (https://github.com/kircherlab/MPRAlib) and archived at Zenodo (https://doi.org/10.5281/zenodo.18173084). MPRAOligoDesign (v0.1.2; Supplemental Code S3) is available at GitHub (https://github.com/kircherlab/MPRAOligoDesign) and archived at Zenodo (https://doi.org/10.5281/zenodo.18173304).