Transcription factors operate across disease loci, with EBNA2 implicated in autoimmunity

EBNA2-occupied genomic sites intersect autoimmune-associated loci

We applied RELI to 213 diseases and phenotypes obtained from the NHGRI GWAS catalog²⁹ and other sources (see Online Methods), revealing nine phenotypes displaying strong EBNA2 association in addition to SLE and MS: rheumatoid arthritis (RA), inflammatory bowel disease (IBD), type 1 diabetes (T1D), juvenile idiopathic arthritis (JIA), celiac disease (CelD), chronic lymphocytic leukemia (CLL), Kawasaki disease (KD), ulcerative colitis (UC), and immunoglobulin glycosylation (IgG) (Supplementary Data Set 3). We designate the seven disorders among these with particularly strong EBNA2 associations (P_c<10⁻⁸) the “EBNA2 disorders.” A recent study performed statistical fine-mapping of the variants for six of the seven EBNA2 disorders (IBD was not included)³⁰. Of the resulting 1,953 candidate causal variants in that study, 130 overlap with EBNA2 ChIP-seq peaks in Mutu B cells (RR=8.7, P_c<10⁻¹³²). Notably, this represents the second-ranked ChIP-seq dataset out of the 1,544 considered in our study, trailing only POLR2A ChIP-seq performed in CD4+ T cells (Supplementary Data Set 3). Thus, the overlap between EBNA2 ChIP-seq peaks and loci associated with the EBNA2 disorders is even stronger when only considering statistically likely causal variants.

Consistent with the SLE results (Figure 1a), the same TFs tend to cluster with distinguishing loci for each disorder (Figure 1b–g, Supplementary Data Set 5). Further, there is a stronger association in EBV-infected than in EBV negative cells for many TFs, and the 10 most associated TFs consistently intersect more strongly in EBV-infected B cells than in other cell types (Figure 1b–g, Supplementary Data Set 5). Hierarchical clustering identifies a core set of 47 TFs binding to 142 risk loci across the seven EBNA2 disorders (Supplementary Figure 4). RBPJ, an established EBNA2 co-factor^31–33, has the most similar binding profile to EBNA2 across loci, as expected.

NFκB proteins RELA, RELB, REL, NFKB1, and NFKB2 comprise many of the strongest associations with EBNA2 disorder loci (Supplementary Data Set 3). We therefore collected the 348 loci associated with at least one of the EBNA2 disorders, and removed the 179 among these loci that contain at least one disease-associated variant located within a ChIP-seq peak for any NFκB protein in EBV-infected B cells. Among the remaining 169 loci, 19 still contain disease-associated variants falling within EBNA2 ChIP-seq peaks (2.15-fold enrichment, P=0.00012), indicating that many of these loci may be occupied by EBNA2 independent of NFκB involvement.

In order to identify candidate EBNA2 co-factors, we isolated EBNA2 disorder-associated variants located within EBNA2 ChIP-seq peaks and evaluated them using RELI. This analysis confirms the importance of RBPJ, followed by members of the basal transcriptional machinery (TBP and p300), and NFκB subunits (which are involved in EBNA2-mediated gene activation³⁴) (Figure 2b). Interestingly, predicted EBNA2 co-factors vary with disease phenotype; for example, EBNA2 and EBNA3C are highly synergistic at the disease loci of three of the EBNA2 disorders (IBD, MS, and CelD), but rarely coincide at loci for the other four diseases (Supplementary Data Set 6).

The particular TFs tend to be shared across the EBNA2 disorders, but the loci they occupy are less frequently shared. No EBNA2-bound locus is associated with all seven EBNA2 disorders; most loci are unique to only one disorder (Figure 2c). Thus, the loci occupied by EBNA2 in each disorder are largely distinct from one another. One counterexample involves the IKZF3 locus encoding the Aiolos TF, a key regulator in B lymphocyte activation³⁵, with genetic variants from five different EBNA2 disorders intersecting EBNA2 ChIP-seq peaks (Supplementary Figure 4).

If changes in gene regulation explain these results, then expression quantitative trait loci (eQTLs), ChIP-seq peaks for Pol-II, and histone marks associated with active gene regulatory regions should be relatively concentrated at the risk loci occupied by EBNA2. These predictions are indeed true for each of the seven EBNA2 disorders (Figure 2d and Supplementary Data Set 3). For example, <1% of all common variants in the human genome are eQTLs in EBV-infected B cell lines (Figure 2d). This value rises to 2.3% for common variants located within open chromatin in EBV-infected B cell lines, and rises further to 2.7% for common variants within EBNA2 ChIP-seq peaks (Figure 2d, upper left panel, bars labeled “Common variants”). Thus, there is a slight trend for a common variant located within an EBNA2 ChIP-seq peak to influence gene expression in EBV-infected B cell lines. Strikingly, this relationship is >10-fold increased for EBNA2 disorder-associated variants - 27.8% of EBNA2 disorder variants that are located within EBNA2 ChIP-seq peaks are also eQTLs, a value significantly greater than EBNA2 disorder variants located within open chromatin in EBV-infected B cell lines (20.5%, P<10⁻⁵, Welch’s one-sided t-test) or EBNA2 disorder variants in general (10.4%, P<10⁻⁸) (Figure 2d, upper left panel, bars labeled “EBNA2 disorder variants”). Similar trends hold for the other data types examined (Figure 2d). In aggregate, these results hint at the potential magnitude of the environmental influence of EBNA2 upon host gene expression within EBNA2 disorder loci in EBV-infected B cells.

EBNA2 participates in allele-dependent formation of transcription complexes at disease risk loci

The observed associations (Figure 1) are genetic if and only if they are driven by causal allele-dependent differences. Since EBNA2 imitates the binding of NOTCH to RBPJ³⁶, genetic variants at these loci could alter the binding of RBPJ (or another TF to which EBNA2 binds) or enable allele-dependent binding of a TF that requires the presence of EBNA2 by modulating the local chromatin environment (Figure 3a). Re-analysis of ChIP-seq data provides a means to identify allele-dependent protein binding events on a genome-wide scale - in cases where a given variant is heterozygous in the cell assayed, both alleles are available for the TF to bind, offering a natural control for one another since the only variable that has changed is the allele. We therefore developed the MARIO (Measurement of Allelic Ratios Informatics Operator) pipeline to identify allele-dependent protein binding by weighing imbalance between the number of sequencing reads for each allele of a given genetic variant, the total number of reads available at the variant, and the number and consistency of available experimental replicates (see Online Methods). MARIO is an easy-to-use, modular tool that extends existing methods^37–40 by (1) calculating a score that explicitly reflects reproducibility across experimental replicates; (2) reducing run-time via utilization of multiple computational cores; and (3) allowing the user to directly provide genotyping data as input. To identify heterozygotes for analysis, we genotyped five EBV-infected B cell lines with available ChIP-seq data and performed genome-wide imputation (see Online Methods). We applied MARIO and a related method, ABC³⁷, to a deeply sequenced (~190 million reads) GM12878 ATAC-seq dataset (GEO accession {"type":"entrez-geo","attrs":{"text":"GSM1155957","term_id":"1155957"}}GSM1155957) and observed strong agreement between the 2,214 resulting scores (Spearman correlation of 0.98 (P<10⁻¹⁵)).

Open in a separate window

Figure 3

Allele-dependent binding of EBNA2 to autoimmune-associated genetic variants

a. Theoretical models presenting possible allele-dependent action of EBNA2. See text for discussion. b. Allele-dependent co-binding of EBNA2 with multiple proteins. ChIP-seq datasets from EBV-infected B cell lines were examined for evidence of allele-dependent binding at heterozygotes. Datasets are sorted by the proportion of EBNA2 GM12878 allele-dependent events (MARIO ARS value > 0.40, see Online Methods) that favor the same allele (X-axis). Values (N) indicate total number of variants. c. Allele-dependent binding of EBNA2 and human proteins at the CD44 locus. Top to bottom: chromosomal band (multi-colored bar), location of EBV-infected B cell line ChIP-seq peaks for various TFs, location of rs3794102 variant, allele-dependent binding events (green bars). X-axis indicates the preferred allele, along with a value indicating the strength of the allelic behavior, calculated as one minus the ratio of the weak to strong reads (e.g., 0.5 indicates the strong allele has approximately twice the reads of the weak allele). d. Allele and EBV-dependent expression of CD44. Allelic qPCR of CD44 expression in EBV-infected and EBV negative Ramos B cells (see key). Fold-change in expression is provided relative to the C allele. Error bars represent standard deviation (n=12: three independent experiments of technical quadruplicates). P-values were calculated using a two-way ANOVA with a Tukey post-hoc test. EBV status and variant genotype were used as the two factors.

We next applied MARIO to 271 ChIP-seq datasets performed in one of the five genotyped cell lines, altogether assessing 98 different molecules. Since EBNA2 binds DNA indirectly as a co-factor, we first asked if the variants displaying EBNA2 allele-dependent binding might coincide with similarly altered binding of other TFs. This analysis revealed strong concordance of allele-dependent binding events both within and across cell types. For example, we identified 68 heterozygous common variants located within allele-dependent EBNA2 GM12878 ChIP-seq peaks. EBF1, whose binding is globally influenced by EBNA2³⁶, has a coincident ChIP-seq peak favoring the same allele at 39 (57%) of these loci, as opposed to only 8 (11%) on the opposite allele (P<10⁻⁴, binomial test, Figure 3b). Similar results were obtained when pairing EBNA2 binding in GM12878 with EBNA2 binding in Mutu cells, with established partners SPI1 and RBPJ, or with ATAC-seq chromatin accessibility data (Figure 3b). Analogous results are obtained for EBNA2 ChIP-seq data in Mutu and IB4 cell lines (Supplementary Figure 5). In total, MARIO confidently identified 21 variants associated with 15 different autoimmune diseases displaying allele-dependent EBNA2 binding in at least one cell type (Table 2, Supplementary Data Set 7). We note that the number of heterozygous autoimmune variants for which EBNA2 prefers one allele over the other is not significantly more than expected by chance (see Online Methods). We also note that several variants might involve the HLA genes, and the current view is that coding alleles in the HLA class II in general are likely (though not certainly) causal for autoimmune diseases. Nevertheless, most of these variants also involve allele-dependent host protein binding, chromatin accessibility, or presence of histone marks such as H3K27ac (Supplementary Data Set 8). Together, these results suggest that many autoimmune-associated variants may act by modifying host gene regulatory programs via altered binding of EBNA2 and additional proteins.

To detect potential downstream effects of allele-dependent EBNA2 binding, we measured genome-wide gene expression levels by RNA-seq in Ramos, an EBV negative B cell line that can support an EBV infection. We confirmed the expected presence or absence of EBNA2 by sequencing (Online Methods) and western blot (Supplementary Figure 6). We identified a total of 80 genes with significant EBV-dependent alterations in gene expression (Supplementary Data Set 9), confirming that EBV modulates the expression of human genes. These results are highly consistent with a previous gene expression study and the literature (see Online Methods).

We next searched for autoimmune-associated variants that might impact EBNA2 binding, resulting in allele-dependent expression of a nearby gene. This analysis was dependent on the small subset of genetic variants satisfying four necessary criteria: the variant must be (1) plausibly causal for an autoimmune disorder; (2) immunoprecipitated by EBNA2 antibodies; (3) heterozygous in the cell line assayed; and (4) proximal to a plausible target mRNA that contains a heterozygous variant in Ramos cells (to detect allele-dependent expression). For example, the 21 EBNA2 variants listed in Table 2 satisfy the first three criteria, but only five satisfy the fourth criterion of being within 50kb of a potential target gene containing a heterozygous variant in the Ramos cell line.

Despite these limitations, our approach identified autoimmune-associated variants displaying allele-dependent EBNA2 binding and allele-dependent expression of a nearby gene. For example, rs3794102, a variant strongly associated with vitiligo (P<10⁻⁹ for case/control association), has significantly skewed allele-dependent binding of eight proteins - EBNA2, its suspected co-factor EBF1³⁶, and chromatin accessibility all favor the non-reference ‘G’ vitiligo risk allele (Figure 3c, Table 2, Supplementary Data Set 8). Intriguingly, the proteins favoring the ‘G’ allele are considered activators, whereas the two proteins that prefer the ‘A’ allele are repressors, suggesting that the variant and virus might act synergistically as an allelic switch. rs3794102, which is located within an intron of SLC1A2 (a gene for which we detect no RNA-seq reads), loops to the promoter of the neighboring CD44 gene based on Hi-C experiments performed in GM12878 cells (Supplementary Figure 7). rs3794102 is also an established eQTL for CD44 in EBV-infected B cell lines (P<10⁻¹¹, ‘MRCE’ dataset, RTeQTL database⁴¹), and particular isoforms of CD44 are dependent on the presence of EBNA2⁴². In our experiments, CD44 expression is 6.8-fold higher in EBV-infected Ramos cells compared to uninfected Ramos cells (P = 0.00015, Supplementary Data Set 9). Further, we identified a heterozygous genetic variant (rs8193) in strong LD with rs3794102 (r²=0.87) located within the CD44 gene body with 12 ‘T’ allele RNA-seq reads and only 5 ‘C’ allele reads in EBV-infected Ramos cells, and no detectable reads in Ramos cells lacking EBV (Supplementary Data Set 10). We independently confirmed this result with allelic qPCR, observing a significant increase in expression for the T relative to the C allele in EBV-infected Ramos cells, with significantly lower levels of expression in the absence of EBV (Figure 3d). CD44 is a transmembrane glycoprotein involved in B cell migration and activation. Taken together, these results suggest that the ‘G’ vitiligo risk allele enhances formation of an EBNA2-dependent gene activation complex, resulting in elevated expression of CD44, and consequent increased B cell migration and/or activation. We also identified a variant (rs947474) associated with T1D and RA (Table 2) located near PRKCQ, another gene with allele- and EBV-dependent expression in our data (Supplementary Data Set 10). Intriguingly, PRKCQ plays an established role in activation of the EBV lytic cycle⁴³. Together, these examples establish that multiple autoimmune variants may alter binding events of protein complexes containing EBNA2 and host proteins, resulting in EBV-controlled allele-dependent host gene expression.

Autoimmune-associated genetic mechanisms in EBV-infected B cells

We next used RELI to rank cell types by their relative importance to each of the EBNA2 disorders, based on the intersection between disease-associated variants and likely regulatory regions in that cell type. This procedure revealed a clear enrichment for EBV-infected B cells in SLE. For example, of the 175 H3K27ac ChIP-seq datasets available, the highest ranked 30 datasets are all from EBV-infected B cell lines (Figure 4a). Analogous results are obtained for “active chromatin marks” (a model based on combinations of various histone marks⁴⁴) (Figure 4b), H3K4me3, and H3K4me1, for SLE and virtually all of the seven EBNA2 disorders (Supplementary Data Set 3, Supplementary Data Set 11). Collectively, these results support the EBV-infected B cell being an origin for genetic risk for each of the seven EBNA2 disorders. This analysis also reveals a likely involvement of other immune cell types in these disorders, including T cells, natural killer cells, and monocytes (Supplementary Data Set 3). Although there are limited TF ChIP-seq data available for these cell types, one or more of the EBNA2 disorders are associated with 17 of the available T cell TF ChIP-seq datasets (Supplementary Data Set 3). Further, several EBNA2 disorder loci appear to be specific to T cells. For example, six MS-associated loci are largely T cell-specific, collectively intersecting 67 T cell ChIP-seq datasets, compared to only 12 EBV-infected B cell datasets for these same loci (Supplementary Data Set 12). Together, these results are consistent with multiple shared regulatory mechanisms acting across autoimmune risk loci, some common between cell types (e.g., B and T cells) and others being exclusive to a certain cell type.

Open in a separate window

Figure 4

Cell types and TFs at disease-associated loci

a. SLE variants significantly intersect H3K27ac-marked regions in EBV-infected B cells. H3K27ac ChIP-seq peaks were collected from 175 different cell lines and types. The Y-axis indicates the negative log of the RELI P-value for the intersection of SLE-associated variants with H3K27ac peaks in each dataset. b. SLE variants intersect active chromatin regions in EBV-infected B cells. Same as (a), but instead using “active chromatin” regions, which are based on combinations of histone marks⁴⁴. c. Global view of RELI results – all diseases against all TFs. Columns and rows show the 94 phenotypes/diseases and 212 TFs with at least one significant (P_c<10⁻⁶) RELI result. Color indicates negative log of the RELI P-value (see key). Disease abbreviations are provided in the main text. d. TFs intersecting breast cancer loci. Intersection between disease loci with TF-bound DNA sequences, as in Figure 1. However, here the cluster of TFs and risk loci instead largely may operate in ductal epithelial cells, independently of EBNA2. The top 20 TFs are shown - full results are provided in Supplementary Data Set 3.

Regulatory Element Locus Intersection (RELI)

RELI takes a set of genetic variants as input, expands the set using LD blocks, and calculates the statistical intersection of the resulting loci with every dataset in a compendium (e.g., ChIP-seq datasets) (Figure 2a and Supplementary Figure 1). In Step 1, sequencing data from 1,000 Genomes⁶⁹ are used to identify all variants with linkage (r²>0.8) to any input variant within each major ancestry (European, African, Asian), thereby assigning them to LD blocks. In Step 2, overlapping genomic coordinates determine whether an observed intersection is recorded between each LD block and each dataset. In Step 3, the expected intersection is estimated between each LD block and each dataset. The most strongly associated variant is chosen as the reference variant for the LD block. A distance vector is generated providing the distance (in bases) of each variant in the LD block from this reference variant. A random genomic variant with approximately matched allele frequencies to the reference variant is then selected from dbSNP⁷⁰, and genomic coordinates of artificial variants are created that are located at the same relative distances from this random variant using the distance vector. Members of this artificial LD block are intersected with each dataset, as was done for the observed intersections. This strategy accounts for the number of variants in the input LD block and their relative distances, while prohibiting ‘double counting’ due to multiple variants in the block intersecting the same dataset. We repeat this procedure 2,000 times, generating a null distribution with stable P-values. The expected intersection distributions are used to calculate Z-scores and P-values for the observed intersection. The final reported P-values are Bonferroni corrected (P_c) for the 1,544 TF datasets tested. We calculate the relative risk by dividing the observed intersection by the mean expected intersection. We also considered a higher-stringency null model that only considers variants located within DNase-seq peaks in any of the 22 available EBV-infected B cell line datasets, which controls for the known association of autoimmune variants and B cell regulatory regions²³.

We validated the RELI procedure as follows. First, we compared the Z-score-based P-values produced by RELI to empirically calculated P-values. We selected 187 ChIP-seq datasets with European SLE GWAS RELI corrected P-values that are evenly distributed between 1 and 10⁻⁷. An upper bound of 10⁻⁷ was chosen due to the amount of time required to run the simulations. Overall, we observe very strong concordance between these 187 empirically-derived P-values and the P-values estimated by RELI (Supplementary Figure 8, Panel A), with a Pearson correlation coefficient of 0.82 (P<10⁻⁴⁵). We also performed 200,000,000 simulations examining the EBNA2 Mutu ChIP-seq vs. European SLE variant relationship. Across these simulations, we observed a maximum of 16 loci intersecting EBNA2 Mutu ChIP-seq peaks (Supplementary Figure 8, Panel B), conservatively setting an empirically-determined P-value lower bound at 5 x 10⁻⁹ and further supporting our estimated P-value of P_c<10⁻²⁴ for the 26 observed locus intersections. To validate our choice of 2,000 simulations, we compared the P-values obtained for the 187 datasets when using 2,000 vs. 5,000,000 simulations. Nearly identical P-values were obtained (Supplementary Figure 8, Panel C).

We also estimated RELI false positive rates. We first generated a “false library” of 1,544 ChIP-seq datasets that match the “real library” by randomly repositioning each peak within the genome. This random “false library” of ChIP-seq results matches the number of datasets, the number of peaks each dataset contains, and the width of those peaks. Upon running the European SLE variants with RELI using 10 different “false libraries”, only one of the 15,440 datasets achieved a P-value less than our P_c<10⁻⁶ threshold (Supplementary Data Set 3). Further, the P-value for this dataset (P_c<10⁻⁸) is much less significant than those for EBNA2, RELA, etc. (P_c ≪ 10⁻²⁰). We thus estimate our overall false positive rate to be ~ 1/15,440 (~0.006%).

Identification of optimal clusters

We identified optimal clusters (red outlines in Figure 1) by comparing the observed number of TF/locus intersections to results from simulations. First, loci (X-axis) and TFs (Y-axis) were sorted in decreasing order of the number of intersections (colored boxes in the heatmap). We then iteratively considered every possible sub-matrix boundary, starting at the upper left corner. In each trial, the total number of intersections is kept fixed, but the locations of the intersecting positions are randomly permuted across loci. A Gaussian null distribution was obtained from 10,000 random trials. P-values were calculated for each sub-matrix by comparing the observed number of intersections within the sub-matrix to the null distribution, using a standard Z-score transformation. The optimal cluster was defined as the sub-matrix with the best P-value.

Cell line genotyping and imputation

We genotyped five EBV-infected B cell lines with available ChIP-seq data (Supplementary Table 2) on Illumina OMNI-5 arrays, as previously described⁷¹. Genotypes were called using the Gentrain2 algorithm within Illumina Genome Studio. Quality control was performed as previously described⁷¹. Quality control data cleaning was performed in the context of a larger batch of non-disease controls to allow for the assessment of data quality. Briefly, all cell lines had call rates >99%; only common variants (minor allele frequency >0.01) were included; and all variants were previously shown to be in Hardy-Weinberg equilibrium in control populations at P>0.0001⁷¹. We performed genome-wide imputation using overlapping 150 kb sections of the genome with IMPUTE2⁷² and a composite imputation reference panel of pre-phased integrated haplotypes from 1,000 Genomes (June 2014). Included imputed genotypes met or exceeded a probability threshold of 0.9, an information measure of 0.5, and the same quality-control criteria described above for the genotyped markers.

Detection of allele-dependent sequencing reads using MARIO

We developed the MARIO (Measurement of Allelic Ratios Informatics Operator) pipeline to identify allele-dependent behavior at heterozygous genetic variants in functional genomics datasets. In brief, the pipeline downloads a set of reads, aligns them to the genome, calls peaks using MACS2 (parameters: --nomodel --extsize 147 -g hs -q 0.01), identifies allele-dependent behavior at heterozygotes within peaks (described below), and annotates the results (Supplementary Figure 9).

To estimate the significance of the degree of allelic imbalance of a given dataset at a given heterozygote, we developed the Allelic Reproducibility Score (ARS), which is based on a combination of two predictive variables: the total number of reads overlapping the variant and the imbalance between the number of reads for each allele. Other variables tested were uninformative (below). The ARS value also accounts for the number of available experimental replicates and the degree to which they agree. ARS values were calibrated using seven TFs with four replicate ChIP-seq experiments available in the same cell line (GM12878 or K562): SPI1 (set 1), SPI1 (set 2), NRSF, REST, RNF2, YY1 and ZBTB33.

ARS values were calculated as follows:

1. Determine the number of reads mapping to each allele of each heterozygous variant in each replicate. We applied our pipeline to each experimental replicate and counted the number of reads for each allele that overlap each heterozygous variant. Insertions and deletions were not considered. All duplicate reads were removed using the “MarkDuplicates” tool from the PICARD software package (https://broadinstitute.github.io/picard/). Before mapping reads using Bowtie2⁷³ (parameters -N 1 --np 0 --n-ceil 10 --no-unal), we masked all common variants in the GrCh37 (hg19) reference genome to N, which removes bias generated by reads carrying non-reference alleles. We designate the allele with the greater number of reads the strong allele, and the other the weak allele (Supplementary Figure 10a).
2. Identify predictive variables of reproducible allele-dependent behavior across replicates. We collected a set of seven datasets, {D}, with each dataset comprised of four experimental replicates, {R} (Supplementary Figure 10b). Each replicate contains a set of variants {V} that are heterozygous in the given cell type. For each of these variants, we calculated the value of four variables {X}: the ratio between the number of weak and strong allele reads, the total number of reads available at the variant, distance to peak center, and peak width.

   We evaluated the performance of each of these variables using a true-positive set of reproducible variants. This set was created by identifying all variants that share the same strong allele across all four replicates (Supplementary Figure 10c). Each variable was assessed based upon its ability to effectively separate reproducible variants from non-reproducible variants (all other variants). The reproducible variants are enriched for allele-dependent behavior, whereas the non-reproducible variants are depleted (Supplementary Figure 10d, left-most panel). Of the four variables tested, two were predictive of reproducible allele-dependent binding: the ratio between the number of weak and strong allele reads (WS_ratio), and the total number of reads available at the variant (num_reads), which we designate the predictive variables.

3. Determine a function mapping the values of the predictive variables to a single ARS value. Our approach accommodates datasets containing any number of experimental replicates and rewards greater agreement between replicates. Within each of the seven datasets in the set {D}, we consider all possible combinations of one, two, or three replicates. Without loss of generality, we describe the procedure for the case of two replicates, which considers the subsets {R₁,R₂}, {R₁,R₃}, {R₁,R₄}, etc. We first identify the set {H} of reproducible variants (as described above) for each subset. We then threshold the WS_ratio into ranges, {(0 – 0.1), (0 – 0.2), (0 – 0.3), … (0 – 1 )}, and for each range, we calculate the fraction of variants that are contained in the reproducible variant set as a function of num_reads (Supplementary Figure 11a). At this stage, this fraction still accounts for all variants, both allele-dependent and non-allele-dependent. We therefore adjust each curve by the normalized cumulative frequency of non-allele-dependent variants within the given range. For example, consider the WS_ratio=0.3 curve (Supplementary Figure 11a). Each point on this curve is divided by a single value representing the normalized cumulative frequency of the non-reproducible variants, which is obtained from the Y-axis at the X=0.3 position in the WS_ratio plot depicted in Supplementary Figure 10d. Before dividing, 1 is added to this value to avoid divide-by-zero errors. Collectively, this approach selectively penalizes non-allele-dependent behavior by accounting for the proportion of non-allele-dependent variants within each curve. These values were averaged across the seven datasets, yielding the final ARS values. This entire procedure is repeated for the cases of one, two, or three available replicates, generating the points shown in Supplementary Figure 11b. Curves were fit to these points using a saturating function:

   ARSw=Aw1+Bw×r-Aw,

   where w is the WS_ratio, r is num_reads, and A_w and B_w are the fitting parameters. The resulting functions yield ARS values for any given heterozygous variant in any dataset, as a function of the number of experimental replicates, the WS_ratio, and num_reads. When multiple replicates are available, we only report an ARS value for a variant if the strong allele is consistent in the majority of cases. A direct interpretation of the ARS values can be seen in the relationship between ARS values and the WS_ratio (Supplementary Figure 11c).

Statistical significance of the number of EBNA2 allele-dependent binding events

We observed a total of 21 cases of allele-dependent EBNA2 binding to an autoimmune risk variant (Table 2). To establish the statistical significance of this observation, we collected the full set of 42 autoimmune-disease associated variants that are (1) located within a ChIP-seq peak in at least one of the three available EBNA2 datasets and (2) heterozygous in the cell type from which that peak was obtained. This set represents all autoimmune variants for which we could have observed allele-dependent EBNA2 binding. We next created a pool of non-autoimmune-associated variants that also meet the above two requirements (resulting in a total of 4,160 variants). For each of the 42 autoimmune variants, we chose a corresponding non-autoimmune variant from this pool, while approximately matching for the total number of EBNA2 ChIP-seq reads in the peak (within 10% of the read count). This procedure thus creates a matched set of 42 non-autoimmune variants that have an equal chance of resulting in allele-dependent EBNA2 behavior. There were a sufficient number of variants to repeat the above procedure 10 times, without replacement. In total, we observed 256 significant EBNA2 allele-dependent binding events across these matched non-autoimmune variant sets, which is not significantly different from the frequency that we observed with the autoimmune variants.

EBV Infection of Ramos cells

All cells were confirmed to be free of mycoplasma infection using PlasmaTest (InvivoGen, San Diego, CA). Wild-type EBV was prepared from supernatants of B95–8 cells cultured in RPMI medium 1640 supplemented with 10% FBS for two weeks. Briefly, the cells were pelleted and the virus suspension was filtered through 0.45 μM Millipore filters. The concentrated virus stocks were aliquoted and stored at -80^oC.

We infected ~2 x 10⁶ Ramos Cells (ATCC CRL-1596) in the presence of growth medium containing 2μg/ml of phytohemagglutinin (PHA) for 4 hours. The infected cells were washed, cultured in growth media, and observed daily for multinuclear giant cell formation and morphological changes characteristic of EBV-infected B cells. After 10 passages, the infection was confirmed by measuring the expression of viral EBNA2 protein levels (Supplementary Figure 6).

RNA-seq

RNA was isolated from Ramos cell lines with and without EBV infection using the mirVANA Isolation Kit (Ambion). RNA sequencing targeting 150 million mappable, 125 base pair reads from paired-end, poly-A enriched libraries was performed at the CCHMC DNA Sequencing and Genotyping Core Facility. Sequencing reads were aligned to the GrCh37 (hg19) build of the human genome using TopHat⁷⁴ and Bowtie2⁷³ with Ensembl⁷⁵ RNA transcript annotations as a guide. In parallel, these data were aligned to the EBV genome (NCBI). As expected, 0 reads mapped in the EBV negative dataset, whereas 7,349 reads mapped in the EBV-infected dataset. 82.8% of the sequence reads aligned specifically to the human transcriptome, with a 2.6% increase in the aligned reads in the EBV negative samples. No abnormal quality control (QC) flags were identified following QC analysis with the software FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). For allelic analysis, sequencing reads were aligned to the GrCh37 (hg19) build of the human genome using Hisat2⁷⁶. Differential expression analysis was performed using Cufflinks⁷⁷.

As additional QC, we further compared our results to a study examining host gene expression changes to EBV infection in primary B cells²⁸. Of the 80 genes whose expression is significantly altered by the presence of EBV in our study, 18 of them are also significantly differentially expressed in this dataset. Further, among the 80 differentially expressed genes we detect, many of them represent classic host genes whose expression is modulated by EBV. Genes whose expression is concordantly activated by EBV include CD44⁷⁸, TNFAIP2⁷⁹, MX1⁸⁰, and IFI44⁸¹; genes with lowered expression include VAV3⁸² and CD99⁸³.

Allelic qPCR

gDNA and RNA were extracted from Ramos cells with and without B95.8 EBV infection using the DNeasy Blood & Tissue Kit (Qiagen) and mirVana miRNA Isolation Kit (Invitrogen), respectively. RNA was treated with DNase using the TURBO DNA-free Kit (Ambion) and converted to cDNA using the High-Capacity RNA-to-cDNA Kit (Applied Biosystems). qPCR was performed with a single set of Taqman genotyping primers (Applied Biosystems) to rs8193 using the ABI 7500 PCR system. Fold change of expression was calculated with 2^−ΔΔCT values, where cDNA was normalized to gDNA.

Statistical analyses

Details on statistical analyses are described in the corresponding sections. For statistical details on RELI and MARIO, see the corresponding sections above. The number of replicates or data points (N) is provided in the Figures and legends. Data are represented as means +/− one standard deviation, unless otherwise noted.

Data availability

RNA-seq data are available in the Gene Expression Omnibus (GEO) database under accession number {"type":"entrez-geo","attrs":{"text":"GSE93709","term_id":"93709"}}GSE93709. Full datasets and results, including disease variants (with alleles) and all RELI and MARIO output, are provided in the Supplementary Material.

We thank: J. Lee, C. Schroeder, Y. Huang, X. Lu, Z. Patel, E. Zoller, and The CCHMC DNA Sequencing and Genotyping Core for experimental support; C. Gunawan, K. Ernst, and T. Hong for analytical support; B. Cobb for administrative support; R. Kopan, C. Karp, W. Miller, J. Whitsett, M. Fisher, A. Strauss, S. Hamlin, L. Muglia, H. Singh, J. Oksenberg, I. Chepelev, S. Waggoner, S. Thompson, and H. Moncrieffe for constructive feedback and guidance; Y. Yuan (University of Penn) and D. Thorley-Lawson (Tufts Institute) for generous donation of cell lines (Mutu and IB4, respectively). Thanks to all of our colleagues who have made their data available to us, without which this project and its results would not have been possible. Funding sources: National Institutes of Health (NIH) R01 NS099068, NIH R21 HG008186, Lupus Research Alliance “Novel Approaches”, CCRF Endowed Scholar, CCHMC CpG Pilot study award, and CCHMC Trustee Awards to M.T.W.; NIH R01 AI024717, NIH U01 HG008666, NIH U01 AI130830, NIH P30 AR070549, NIH R24 HL105333, NIH KL2 TR001426, NIH R01 AI031584, NIH P30 DK078392, Kirkland Scholar Award and US Department of Veterans Affairs I01 BX001834 to J.B.H.; NIH R01 DK107502 to L.C.K; NIH DP2 GM119134 to A.B.