Design of Agilent oligonucleotide library
HepG2 pilot library
For the HepG2 pilot library, we collected two replicates of DNase I hypersensitivity data derived from HepG2 cells (ENCODE narrowPeak BED files: ENCFF505SRS and ENCFF268DTI, hg19 human genome build)55. For each replicate, we collapsed overlapping peaks using bedtools merge (parameters “-o collapse -c 2,3,7”). Then, we identified peaks that intersected between the two replicates, merged these peaks, and removed the subset of merged peaks that overlap promoters (defined as regions ±2,500 nt around any annotated TSS). The resulting distribution of peak sizes was such that 97% of peaks were ≤200 bp in length. We therefore centred the designed oligonucleotides at each merged DNase peak, consistent with the known region of maximal regulatory activity18, and added ±100 bp to either side. This procedure resulted in a set of 66,017 cCREs. For this pilot library, we sought to evaluate cCREs which overlapped a wide range of putative transcription factor binding sites. We therefore intersected these potential enhancers with wgEncodeRegTfbsClusteredWithCellsV3.bed.gz56 in order to count the number of putative HepG2 transcription factor binding sites intersecting these cCREs. We uniformly and randomly sampled ~1,834 cCREs with 0-1, 1–5, 5–10, 10–20, and >20 TFBSs for a total of 9,172 elements. Including 50 positive and 50 negative controls from each of two previous studies10,19 resulted in a total of 9,372 elements. These 171-bp controls from previous work were linked downstream of a 29-nt random sequence GGTGCTCGATTTAATTTCGCCGACGTGAT to match the 200-bp sequence length of cCREs. For the final oligonucleotide library, each element was linked to the 5′ adaptor AGGACCGGATCAACT and 3′ adaptor CATTGCGTGAACCGA, designing two 230-bp oligonucleotides per element to minimize the impact of oligonucleotide synthesis errors.
K562 pilot library
An analogous procedure was followed for the K562 pilot library as in ‘HepG2 pilot library’, with the following modifications: (1) ENCODE narrowPeak BED files ENCFF027HKR and ENCFF154JCN (hg38 human genome build)55 were used; (2) merging these peaks resulted in 34,367 potential enhancers; (3) after intersecting the peaks with K562 transcription factor binding sites, we sampled ~1,278 enhancers from each transcription factor binding site bin to test a total of 6,394 cCREs; (4) 250 additional negative controls were chosen by dinucleotide shuffling 250 random potential enhancers possessing 1–5 TFBSs; (5) positive and negative controls were chosen from CRISPRi screens20,21, a previous MPRA18, and select loci of interest such as α-globin and β-globin; (6) a total of 7,500 elements were tested; and (7) controls were already 200 bp in length, requiring no addition of a random sequence.
HepG2 large-scale library
Following the procedures outlined in ‘HepG2 pilot library’, we tested all 66,017 previously identified cCREs in both orientations. For human protein-coding gene promoters, we extracted the average signal across cell types in TPM for each CAGE peak listed in hg19.cage_peak_phase1and2combined_tpm_ann.osc.txt.gz from the FANTOM5 consortium57,58. The first exons of all protein-coding gene transcripts were collected from Ensembl v.83 (hg38 genome build)59, transformed into hg19 coordinates using liftOver60, and then intersected with the CAGE peaks to identify a single promoter per gene corresponding to the promoter with the maximal average TPM. To select the final 200-bp oligonucleotide for testing, we identified the centre of the promoter DNase peak on the basis of the HepG2 DNase peaks merged across replicates (described in ‘HepG2 pilot library’). In the scenario in which no DNase peak overlapped the promoter, we centred on the midpoint of the CAGE peak. In the scenario in which neither a DNase nor CAGE peak existed, we centred on the TSS defined by the Ensembl annotation. This resulted in a total of 19,104 protein-coding gene promoters, of which 6,181 were centred on a DNase peak, 9,735 were centred on a CAGE peak and 3,188 were centred on a Ensembl TSS definition. The oligonucleotide tested included the ±100-bp window around this central position in the sense orientation with respect to the gene. A random subset of 12,411 promoters were also tested in the antisense orientation. We tested 102 positive and 102 negative controls from a previous study19 as well as 175 dinucleotide shuffled negative controls in both orientations. These shuffled controls were derived from shuffling a random subset of 175 DNase peaks. This resulted in a library consisting of 164,307 elements, for which we ordered one 230-bp oligonucleotide per element.
K562 large-scale library
To acquire a set of DNase peaks for testing, we used the ‘optimal peak’ calls derived from processing ENCODE experiment ID: ENCSR000EOY through the ENCODE DCC Irreproducible Discovery Rate (IDR) pipeline, available at https://github.com/ENCODE-DCC/atac-seq-pipeline (generously provided by A. Kundaje). Removing DNase peaks overlapping human promoters resulted in 87,618 potential enhancers tested in both orientations. The promoters tested were identical to those described in ‘HepG2 large-scale library’, except that it included all 19,104 promoters tested in both orientations. We tested 50 positive and 200 negative controls from a previous MPRA study18 as well as the same 250 dinucleotide shuffled negative controls as tested in ‘K562 pilot library’. Finally, 14,918 tiles not overlapping DNase peaks, and subsampled from the ±1 Mb region around the following 7 genetic loci: GATA1, MYC, HBE1, LMO2, RBM38, HBA2, and BCL11A, were chosen using our representative subset selection approach (described in the ‘Representative subset selection’ section below) and tested in both orientations. This resulted in a library consisting of 243,780 elements, for which we ordered one 230-bp oligonucleotide per element.
WTC11 large-scale library
To acquire a set of DNase peaks for testing, we used the peak calls derived from applying the hotspot2 pipeline (https://github.com/Altius/hotspot2) at FDR = 0.05 to ENCODE experiment ID: ENCSR785ZUI (generously provided by R.Sandstrom)61. This resulted in two independent replicates, which were merged into a unified set using the procedures described in ‘HepG2 pilot library’. Removing DNase peaks overlapping human promoters resulted in 83,201 potential enhancers. Together with the 19,104 promoters described in ‘HepG2 large-scale library’, these elements were subsampled to select 30,121 potential enhancers and 7,500 promoters using our representative subset selection approach described below, and tested in both orientations. We also tested 100 positive and 100 negative controls from a previous study24 as well as 100 dinucleotide shuffled negative controls, which were derived from shuffling 100 random sequences from our set of 30,121 potential enhancers. This resulted in a library consisting of 75,542 cCREs, for which we ordered one 230-bp oligonucleotide per element.
Joint library tested in HepG2, K562, and WTC11 cells
Given the measured potential enhancers from the forward orientations in each of the HepG2, K562 and WTC11 large-scale libraries, we binned each set of cCREs into ten equally sized bins spanning the range of measured log2(RNA/DNA) values in the selected cell type. We randomly sampled an approximately equal number from each bin, resulting in 19,000 HepG2, 18,958 K562 and 18,946 WTC11 potential enhancers. A similar procedure was followed with sense-oriented promoters, except that the ten bins were established on the basis of the mean log2(RNA/DNA) across all three cell types (that is, instead of performing the procedure independently in each cell type as before), and the top 1,000 promoters exhibiting the greatest variance across three cell lines were also selected for testing. This resulted in the selection of 2,396 out of 19,104 promoters. We also tested 181 positive and 169 negative HepG2 controls from a previous study19, 50 positive K562 controls from a previous study18, and 300 dinucleotide shuffled negative controls. The shuffled controls originated from selecting 100 shuffled controls from each of the three cell types. This resulted in a library consisting of 60,000 cCREs, for which we ordered one 230-bp oligonucleotide per element.
Representative subset selection
Given the limited number of testable elements in the large-scale K562 and WTC11 libraries, we designed a subset selection procedure to more optimally sample a non-redundant subset of elements associated with diverse biochemical features. For K562 cells, we used ground sets of non-overlapping 200-bp windows uniformly covering each of 7 genetic loci; for WTC11 cells, we used ground sets of 83,201 potential enhancers and 19,104 promoters. To perform representative subset selection with these ground sets, we utilized an objective function called facility location. This submodular set function can be optimized using a greedy algorithm, and yields a subset of elements that covers the epigenetic diversity of the ground set62. The facility location function is given as:
$$f(A)=\sum _{v\in V}\mathop{\max }\limits_{a\in A}\varphi (v,a)$$
where V is the ground set, A is a subset of V with k elements and φ is a nonnegative similarity function. Optimization of the facility function was performed using the Python package apricot (https://github.com/jmschrei/apricot/)63. For this study, we set k = 2,231 for each of the 7 loci in K562 cells, k = 30,121 for WTC11 potential enhancers, and k = 7,500 for promoters in WTC11. From the 7 loci, we then filtered out the tiles that overlapped DNase peaks as they had already been tested, and then subsampled to ~2,131 tiles per locus to retrieve 14,918 tiles among all loci.
To assess the pairwise similarity of each element, we utilized hundreds of ENCODE histone and transcription factor ChIP–seq experiments derived from K562 and WTC11-H1 embryonic stem cells (Supplementary Table 6). For each 200-bp tile in the ground set, we computed the mean signal for each ChIP–seq dataset, resulting in a vector of biochemical measurements for each 200-bp tile. We used the Pearson correlation coefficient as a similarity function given these ChIP–seq features.
Generation of MPRA libraries
The MPRA libraries were generated as previously described16. In brief, the Agilent oligonucleotide pool was amplified by 5-cycle PCR using forward primer (pLSmP-enh-f, Supplementary Table 13) and reverse primer (minP-enh-r, Supplementary Table 13) that adds a minimal promoter and spacer sequences downstream of the cCRE. The amplified fragments were purified with 0.8x AMPure XP (Beckman coulter), and amplified for 15 additional cycles using the forward primer (pLSmP-enh-f) and reverse primer (pLSmP-bc-primer-r, Supplementary Table 13) to add 15 bp of random sequence that serves as a barcode. The amplified fragments were then inserted into SbfI/AgeI site of the pLS-SceI vector (Addgene, 137725) using NEBuilder HiFi DNA Assembly mix (NEB), followed by transformation into 10-beta competent cells (NEB, C3020) using the Gemini X2 machine (BTX). Colonies were allowed to grow up overnight on carbenicillin plates and midiprepped (Qiagen, 12945). For HepG2 and K562 pilot libraries, we collected approximately 1 million and 1.3 million colonies, so that on average 50 and 100 barcodes were associated with each cCRE, respectively. For HepG2, K562 and WTC11 large-scale libraries, we collected approximately 8 million, 12 million and 3 million colonies aiming to associate approximately 50, 50 and 40 barcodes per cCRE, respectively. For the joint library, we collected approximately 3.3 million colonies, aiming to associate approximately 55 barcodes per cCRE. To determine the sequences of the random barcodes and their association to each cCRE, the cCRE-mP-barcodes fragment was amplified from each plasmid library using primers that contain flowcell adapters (P7-pLSmP-ass-gfp and P5-pLSmP-ass-i17, Supplementary Table 13). The fragment was then sequenced with a NextSeq mid-output 300 cycle kit using custom primers (Read 1, pLSmP-ass-seq-R1; Index read, pLSmP-ass-seq-ind1; Read 2, pLSmP-ass-seq-R2, Supplementary Table 13).
Cell culture, lentivirus packaging and titration
HepG2 (ATCC, HB-8065) and K562 (ATCC, CCL-243) cell culture were performed as previously described10. WTC11 human iPS cells (RRID:CVCL_Y803) were cultured in mTeSR plus medium (Stemcell technologies, 100-0276) and passaged using ReLeSR (Stemcell technologies, 100-0484), according to the manufacturer’s instructions. WTC11 cells were used for the MPRA experiments at passage 43–49. Cells were not authenticated or checked for mycoplasma contamination. Lentivirus packaging was performed using HEK293T (ATCC, CRL-3216), as previously described with modifications16. In brief, 50,000 cells per cm2 HEK293T cells were seeded in T175 flasks and cultured for 48 h. The cells were co-transfected with 7.5 μg per flask of plasmid libraries, 2.5 μg per flask of pMD2.G (Addgene 12259), and 5 μg per flask of psPAX2 (Addgene 12260) using EndoFectin Lenti transfection reagent (GeneCopoeia) according to the manufacturer’s instructions. After 8 h, cell culture media was refreshed and ViralBoost reagent (Alstem) was added. The transfected cells were cultured for 2 days to complete lentivirus packaging. The lentivirus libraries in the culture media were separated from the HEK293T cells and concentrated using the Lenti-X concentrator (Takara) according to the manufacturer’s protocol. To measure DNA titre for the lentiviral libraries in HepG2, K562, or WTC11, cells were seeded at 1 × 105 cells per well in 24-well plates and incubated for 24 h. Serial volume (0, 2, 4, 8, 16 and 32 μl) of the lentivirus was added along with Polybrene at a final concentration of 8 μg ml−1. The infected cells were cultured for three days and then washed with PBS three times. Genomic DNA was extracted using the Wizard SV genomic DNA purification kit (Promega). Multiplicity of infection (MOI) was measured as relative amount of viral DNA (WPRE region, forward; 5′-TACGCTGCTTTAATGCCTTTG-3′, reverse; 5′-GGGCCACAACTCCTCATAAAG-3′) over that of genomic DNA (intronic region of LIPC gene, forward; 5′-TCCTCCGGAGTTATTCTTGGCA-3′, reverse; 5′-CCCCCCATCTGATCTGTTTCAC-3′) by quantitative PCR using SsoFast EvaGreen Supermix (Bio-Rad), according to the manufacturer’s protocol.
Lentiviral infections and DNA and RNA barcode sequencing
For the HepG2 and K562 pilot libraries, 2.4 M HepG2 or 10 M K562 cells per replicate were seeded in 10 cm dishes or T75 flasks, respectively, and incubated for 24 h. The HepG2 and K562 cells were infected with the lentiviral libraries along with 8 μg ml−1 Polybrene, with an estimated MOI of 50 or 10, respectively. The higher MOI in HepG2 is due to these cells being adherent compared to K562 that grow in suspension. For the large-scale HepG2 library, 15 M HepG2 cells per replicate were seeded in 3× 15 cm dishes (5 million per dish), incubated for 24 h, and infected with the library along with 8 μg ml−1 Polybrene, with an estimated MOI of 50. For the large-scale K562 library, 85 million K562 cells per replicate were seeded in 3 T225 flasks (28.3 M per flask), incubated for 24 h, and infected with the library along with 8 μg ml−1 Polybrene, with an estimated MOI of 10. For the large-scale WTC11 library, 38.4 million WTC11 cells per replicate were seeded in four 10 cm dishes (9.6 M per dish), incubated for 24 h, and infected with the library along with 8 μg ml−1 Polybrene, with an estimated MOI of 10, due to higher MOIs being lethal for these cells. For the joint library, 5 million HepG2, 28 million K562 and 38.4 million WTC11 cells were infected with the estimated MOI of 50, 10 and 10, respectively. For each experiment, three independent infections were performed to obtain three biological replicates. After three days of culture, genomic DNA and total RNA were extracted from the infected cells using AllPrep DNA/RNA mini kit (Qiagen), and sequencing library preparations were performed as previously described16. The libraries were then sequenced with a NextSeq high-output 75 cycle kit using custom primers (Read 1, pLSmP-5bc-seq-R1; Index1 (unique molecular idenitifier (UMI) read), pLSmP-UMI-seq; Index2, pLSmP-5bc-seq-R2; Read 2, pLSmP-bc-seq; Supplementary Table 13)16.
MPRA processing pipeline
Associating barcodes to designed elements
For each of the barcode association libraries, we generated FASTQ files with bcl2fastq v.2.20 (parameters “–no-lane-splitting –create-fastq-for-index-reads –use-bases-mask Y*,I*,I*,Y*”), splitting the sequencing data into paired-end index files delineating the barcodes (I1 and I2) and paired-end read files delineating the corresponding element linked to the barcode (R1 and R2). These files were used to associated barcodes to elements using the association utility of MPRAflow 1.016 (run as: nextflow run association.nf –fastq-insert “R1.fastq.gz” –fastq-insertPE “R2.fastq.gz” –fastq-bc “I1.fastq.gz” –fastq-bcPE “I2.fastq.gz” –aligner “bt2_strand” –design “designed_sequences.fa”). Here, designed_sequences.fa was a FASTA file incorporating all of the element sequences that had been ordered from the corresponding Agilent library, and bt2_strand was used to map elements in a strand-aware fashion to accommodate the existence of elements tested in both orientations. The final output of this utility was a filtered_coords_to_barcodes.pickle file mapping barcodes to elements.
Replicates, normalization and RNA/DNA activity scores
For each of the indexed DNA and RNA libraries, we demultiplexed the sequencing run and generated Fastq files with bcl2fastq v.2.20 (parameters “–barcode-mismatches 2 –sample-sheet SampleSheet.csv –use-bases-mask Y*,Y*,I*,Y* –no-lane-splitting –minimum-trimmed-read-length 0 –mask-short-adapter-reads 0”), where SampleSheet.csv catalogued the correspondence between the index sequence and DNA or RNA replicate sample of origin. In several instances, the “–barcode-mismatches 2” resulted in an index assignment clash, requiring us to instead use “–barcode-mismatches 1”. These commands split the sequencing data into paired-end read files delineating the barcodes (R1 and R3) and a file indicating the UMI (R2) for each DNA or RNA replicate sample. We compiled a table of these files to indicate the 3 RNA and 3 DNA files for each of the three replicates in the file experiment.csv. Finally, we used the count utility of MPRAflow 1.016 (run as: nextflow run count.nf –e “experiment.csv” –design “designed_sequences.fa” –association “filtered_coords_to_barcodes.pickle”) to compute activity scores for each element and replicate as log2(RNA/DNA). Elements with which were measured with fewer than 10 independent barcodes were removed to reduce the impact of measurement noise in downstream analysis. This filter led to the following number of retained elements: (1) HepG2 pilot library, 9,153/9,372 (97.7%); (2) K562 pilot library, 7,323/7,500 (97.6%); (3) HepG2 large-scale library, 139,886/164,307 (85.1%); (4) K562 large-scale library, 226,255/243,780 (92.8%); (5) WTC11 large-scale library, 56,093/75,542 (74.2%); (6) HepG2 joint library, 56,018/60,000 (93.4%); (7) K562 joint library, 56,008/60,000 (93.3%); and (8) WTC11 joint library, 55,983/60,000 (93.3%). To combine the data from all three replicates, the distribution of activity values was normalized to the median activity value within each replicate, and then the activity values were averaged across the three replicates.
Regression modelling
Biochemical model features
We extracted all transcription factor ChIP–seq, histone ChIP–seq, DNase-seq, and ATAC-seq bigWig files available for HepG2, K562, and WTC11 cells for the hg38 human genome assembly under ‘released’ ENCODE status56. To account for the lack of WTC11 data available, we also collected all such datasets for H1-ESCs for inclusion in the predictive model. This resulted in 1,506 bigWig files for HepG2 cells; 1,206 files for K562 cells; and 277 files for WTC11/H1-ESCs (Supplementary Table 6). For each candidate element aside from controls, we computed the mean bigWig signal extracted from the corresponding region of the human genome using bigWigAverageOverBed60. All data was right-skewed, and was therefore log-transformed (that is, after adding a pseudocount of 0.1) to approximate a normal distribution. Finally, for each cell type, multiple replicates corresponding to the same ‘experiment target’ (Supplementary Table 6) were averaged to compute the consensus signal for each target in each cell type. This led to a total of 655 HepG2 features, 447 K562 features and 122 WTC11/H1-ESC features considered by the models.
Sei and EnformerMPRA model features
For the large-scale libraries, Sei37 and Enformer36 were used to predict element activity in both orientations (that is, including adaptors in a fixed orientation to simulate the MPRA experiment). The resulting 21,907 Sei and 5,313 Enformer predictions for each of the two orientations were averaged. For the joint library, Sei and Enformer were used to predict element activity in only the forward/sense orientation, and the resulting human predictions were carried forward as features. As Sei requires an input sequence length of 4,000 bp and Enformer requires one of 196,608 bp, all elements were extended with “N” padding in both directions while centring on the element sequence.
Data pre-processing and model training
For each of the three large-scale libraries, the log2(RNA/DNA) scores for each element were averaged among both orientations in which the element was tested, and then randomly assigned to one of ten cross-validation folds (Supplementary Table 8). All predictive features (that is, biochemical features from the matched cell type, or all Enformer features) were z-score-normalized to scale the features similarly. This enabled a direct comparison of coefficients among features derived from the resulting linear models. As described before10,14,32, for each regression task we optimized the λ regularization hyperparameter using tenfold cross-validation, and then used the optimal value for λ to train 10 lasso regression models, each on 9 of the 10 folds of data, to evaluate the performance of each model on the held-out fold. To evaluate the most relevant features selected, we trained a lasso regression model on the full dataset and visualized the 30 coefficients with the greatest magnitude. A similar strategy was used for data from the joint library tested in all three cell types, ensuring that the same element measured in different cell types was always assigned to the same fold (Supplementary Table 12).
Training MPRALegNet
The LegNet architecture35 was adapted to the training data in the following ways: (1) to account for longer sequences but smaller training set size compared to the original LegNet, we added max pooling layers after each local block; (2) the kernel size and number of blocks were selected to match the model’s receptive field to the sequence length; (3) weight decay during training was increased to prevent model overfitting; (4) gradient clipping was used to avoid gradient explosion during the one-cycle learning rate policy (see Supplementary Methods for more details). For the training-time augmentation, each sequence was provided twice, both in the forward and reverse complement orientations, along with their corresponding measured element activity scores. At test time, the model predicted four scores: (1) the scores for the forward and reverse element orientations relative to the fixed flanking (adaptor) regions; and (2) as a further augmentation, the scores for the full reverse complement sequences (that is, obtaining the reverse complementary sequence of the element and adaptor regions together) from step (1). The final prediction represented an average of these four values. The reproducible code, including implementation and complete parameter settings, is available on GitHub (https://github.com/autosome-ru/human_legnet) and Zenodo (https://zenodo.org/records/10558183 and https://zenodo.org/records/13908857).
Interpreting motifs identified by MPRALegNet
As a step towards motif interpretation, ISM was performed for all possible single nucleotide variants on each 200-bp sequence. Owing to the nature of our cross-validation strategy, for each sequence there were nine models for which the sequence was held out during training. ISM scores were generated for every sequence by averaging the predictions from these nine models. The average reference sequence prediction was then compared with that of the alternative sequence32. We then interrogated our ISM scores to identify the most pertinent motifs associated with changes in variant activity using TF-MoDISco-lite v.2.0.4 (https://github.com/jmschrei/tfmodisco-lite), a more efficient version of TF-MoDISco38. The TF-MoDISco-lite algorithm was used with default settings and similar seqlet patterns were matched against JASPAR 2022 CORE vertebrate non-redundant database64 using Tomtom65.
Modelling dose-dependent and combinatorial motif effects learned by MPRALegNet
A non-redundant set of positional weight matrices (PWMs) from each cell type, as ranked by TF-MoDISco-lite (Extended Data Fig. 6), were extracted and scanned (using FIMO 5.5.466, parameters “–text –thresh 0.001”) along each promoter and potential enhancer that was tested bidirectionally. The motif scans were summarized into a matrix of counts for each transcription factor and element tested, as well as log-likelihood (sum of the log(probabilities)), reflecting the likelihood of a given transcription factor binding the element while considering all TFBS instances in both orientations and their respective binding affinities. We then performed an analysis of homotypic (that is, dose-dependent effects for a single transcription factor) as well as heterotypic (that is, combinatorial effects among pairs of transcription factors) for the top 10 activating transcription factors of each cell type.
For homotypic analysis, we plotted the median element activity (that is, both predicted and observed) for elements possessing 0, 1, 2, 3, 4, or 5 motifs, filtering away elements with ≥1 site to any of the other top 10 transcription factors to reduce the chances of a confounding effect. Groups with a sample size of <10 were also filtered out to minimize the impact of noise. The expected dose-dependent responses (for example, dashed lines in Fig. 3f) were computed using linear regression models examining the relationship between either the observed or MPRALegNet-predicted MPRA activity and the number of TFBSs, given log-transformed and untransformed space to model either multiplicative or additive effects, respectively. The expected trend for multiple sites was extrapolated on the basis of the slope and intercept terms of these linear models.
For heterotypic analysis, we evaluated every pair of the 10 activating motifs, isolating cases in which the element possessed 0 counts of both transcription factors, 1 count of one transcription factor or the other, or 1 count each of the first and second transcription factor. Again, all elements were filtered to those with ≥1 site to any of the other top 10 transcription factors other than the transcription factor pair considered. To further account for confounding effects that could be attributable to all other transcription factors (that is, including those beyond the top 10), we computed the residuals from a linear model which considered the log-likelihood values for all other transcription factors besides the pair of transcription factors under consideration. We call these ‘adjusted log2(RNA/DNA)’ (for example, y axis in Fig. 3h) because they removed variability explained by the binding affinities and occurrences of other transcription factors. Finally, a regression model was fit independently to the predicted and observed activity scores. This model sought to predict activity as a function of the presence of TF1, TF2 or an interaction term (TF1 × TF2). The coefficient for the interaction term represented the strength of the super-multiplicative effect (that is, if the coefficient was positive) or the sub-multiplicative effect (that is, if the coefficient was negative)41,42.
Prediction using MPRALegNet
To generate predictions on an arbitrary sequence, we recommend generating predictions using all 90 pretrained models (considering test-time sequence augmentations such as orientation and shifting for extra precision), and then averaging the predictions to achieve the final prediction. We recommend replacing the fixed 15-bp adaptors with the surrounding natural genomic sequence context whenever available, to reduce the chances that artefactual motifs occurring at the adaptor-sequence boundaries could bias the results.
Calculation of element specificity scores
To compute element specificity scores (ESSs) using activity scores from the joint library, log2(RNA/DNA) values from each cell line were first z-score-transformed. Then an ESS for each element was computed by subtracting the element’s score in each cell type by the mean element score across cell types. A full table of ESSs is provided (Supplementary Table 10).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
