Genome-wide association study identifies multiple susceptibility loci for diffuse large B cell lymphoma.

Diffuse large B cell lymphoma (DLBCL) is the most common lymphoma subtype and is clinically aggressive. To identify genetic susceptibility loci for DLBCL, we conducted a meta-analysis of 3 new genome-wide association studies (GWAS) and 1 previous scan, totaling 3,857 cases and 7,666 controls of European ancestry, with additional genotyping of 9 promising SNPs in 1,359 cases and 4,557 controls. In our multi-stage analysis, five independent SNPs in four loci achieved genome-wide significance marked by rs116446171 at 6p25.3 (EXOC2; P = 2.33 × 10(-21)), rs2523607 at 6p21.33 (HLA-B; P = 2.40 × 10(-10)), rs79480871 at 2p23.3 (NCOA1; P = 4.23 × 10(-8)) and two independent SNPs, rs13255292 and rs4733601, at 8q24.21 (PVT1; P = 9.98 × 10(-13) and 3.63 × 10(-11), respectively). These data provide substantial new evidence for genetic susceptibility to this B cell malignancy and point to pathways involved in immune recognition and immune function in the pathogenesis of DLBCL.

Introductory Paragraph

Diffuse large B-cell lymphoma (DLBCL) is the most common lymphoma subtype and is clinically aggressive. To identify genetic susceptibility loci for DLBCL, we conducted a meta-analysis of three new genome-wide association studies (GWAS) and one prior scan, totaling 3,857 cases and 7,666 controls of European ancestry, with additional genotyping of nine promising SNPs in 1,359 cases and 4,557 controls. In our multi-stage analysis, five independent SNPs in four loci achieved genome-wide significance marked by rs116446171 at 6p25.3 (EXOC2; P=2.33×10−21), rs2523607 at 6p21.33 (HLA-B; 2.40×10−10), rs79480871 at 2p23.3 (NCOA1; P=4.23×10−8), and two independent SNPs, rs13255292 and rs4733601, at 8q24.21 (PVT1; P=9.98×10−13 and P=3.63×10−11, respectively). These data provide substantial new evidence for genetic susceptibility to this B-cell malignancy, and point towards pathways involved in immune recognition and immune function in the pathogenesis of DLBCL.

Diffuse large B-cell lymphoma (DLBCL), the most common subtype of non-Hodgkin lymphoma (NHL)[1], has an aggressive clinical course[2]. The risk of DLBCL is increased in individuals with a family history of NHL (odds ratio (OR)=1.4; 95%CI 1.1–2.0)[3], supporting a genetic contribution. Also, relatives of DLBCL patients are at elevated risk for both DLBCL (RR=9.8, 95%CI 3.1–31) and Hodgkin lymphoma (HL, RR=2.0, 95%CI 1.05–4.0), but not indolent lymphomas[4]. Among candidate gene studies investigating susceptibility to DLBCL, only one locus, the LTA252G/TNF-308A haplotype on chromosome 6p21, reached genome-wide significance (P=2.9×10−8)[5]. In small GWAS of all NHL subtypes combined, no conclusive loci for NHL or DLBCL were identified in individuals of European background[6-9], whereas a recent study conducted in East Asia identified a locus at 3q27[10].

To discover new DLBCL susceptibility loci, in stage 1, we genotyped 2,878 DLBCL cases and 2,854 controls of European ancestry from 22 studies using the Illumina OmniExpress Beadchip (Online Methods; Supplementary Table 1; Supplementary Figure 1). A total of 5,346 (93.3%) samples and 611,844 SNPs successfully passed rigorous quality control criteria (Online Methods; Supplementary Table 2). To augment the number of controls, genotype data from 3,536 cancer-free controls previously analyzed with the Omni2.5 SNP microarray were folded into the analytical build[11], resulting in a total of 2,661 cases and 6,221 controls for the stage 1 GWAS analysis (Supplementary Table 2).

In stage 1, with adjustment for gender, age and four eigenvectors (Online Methods), we observed an enrichment of SNPs with smaller P-values compared to the null distribution in the Q-Q plot with a lambda of 1.016 (Supplementary Figure 2). Two SNPs exceeded the threshold for genome-wide significance (P<5×10−8) whereas 20 SNPs showed highly suggestive associations (P<5×10−7) (Supplementary Figure 3). All but one SNP mapped to the HLA region of chromosome 6 (29.5Mb to 33.2Mb on Human Genome version 19 coordinates).

In stage 2, we included data from two unpublished and previously genotyped GWAS (GELA/EPIC and Mayo) plus one published GWAS (UCSF[7]), totaling 1,196 DLCBL cases and 1,445 controls (Online Methods; Supplementary Tables 1, 3). Because different genotyping platforms were used, we imputed common SNPs for each study based on the 1000 Genomes Project release version 3[12] and IMPUTE2[13] (Supplementary Table 4). In meta-analysis of all genotyped and high-quality imputed SNPs from stages 1 and 2 (N=8,363,971), we identified 19 SNPs at genome-wide significance (P<5×10−8) (Supplementary Table 5) and 134 SNPs at a suggestive level of significance (P<5×10−7) (Supplementary Table 6); 123 of the 153 total SNPs mapped to the HLA region on chromosome 6. Based on these results, we selected and successfully designed TaqMan primers for eight promising SNPs (P<5×10−6) outside the HLA region and one SNP from the HLA region for stage 3 de novo genotyping in an additional 1,359 DLBCL cases and 4,557 controls (Online Methods; Supplementary Tables 1, 3).

In a meta-analysis of all three stages (Supplementary Table 7), we identified four non-HLA SNPs in three novel loci at 6p25.3 (rs116446171, P=2.33×10−21) near EXOC2, 8q24.21 (rs13255292, P=9.98×10−13; rs4733601, P=3.63×10−11) near PVT1 and MYC, and 2p23.3 (rs79480871, P=4.23×10−8) near NCOA1 (Table 1; Figures 1a–c). The two 8q24.21 SNPs displayed minimal linkage disequilibrium (LD, r2=0.03 in 1000 Genomes CEU population). Furthermore, in conditional analysis, both rs13255292 (conditional OR=1.22, P =1.39×10−12) and rs4733601 (conditional OR=1.18, P =2.84×10−10) remained genome-wide significant; together these data support the presence of two independent SNPs associated with DLBCL at 8q24.21. We also observed two suggestive SNPs (P<5×10−7) (Supplementary Table 8), one at 5q31.3 (rs79464052, P=5.57×10−8) in ARAP3 (Supplementary Figure 4), and one at 3q13.33 (rs2681416), although the latter SNP did not replicate in stage 2 or 3.

Table 1

Association of novel loci and new independent SNPs wth risk of diffuse large B-cell lymphoma (DLBCL)

Location	Nearest gene(s)	SNP	Positiona	Risk alleleb	Other allele	RAFc	Stage	No. Cases/No. controls	OR	(95% CI)	P	Phet	I2
6p25.3	EXOC2	rs116446171	484,453	G	C	0.019	Stage 1	2,661/6,220	2.26	(1.82–2.81)	1.48×10−13
						0.018	Stage 2	1,194/1,443	2.70	(1.84–3.96)	3.99×10−7
						0.019	Stage 3	1,351/4,460	1.78	(1.29–2.46)	0.00040
							Combined	5,206/12,123	2.20	(1.87–2.59)	2.33×10−21	0.17	32.82
8q24.21	PVT1	rs13255292	129,076,573	T	C	0.321	Stage 1	2,661/6,221	1.19	(1.11–1.28)	1.25×10−6
						0.315	Stage 2	1,195/1,444	1.30	(1.14–1.47)	4.29×10−5
						0.330	Stage 3	1,322/4,498	1.22	(1.09–1.36)	0.001
							Combined	5,178/12,163	1.22	(1.15–1.29)	9.98×10−13	0.37	8.30
		rs4733601	129,269,466	A	G	0.477	Stage 1	2,661/6,221	1.19	(1.11–1.27)	4.22×10−7
						0.479	Stage 2	1,196/1,445	1.19	(1.05–1.33)	0.004
						0.487	Stage 3	1,337/4,523	1.19	(1.07–1.32)	0.0016
							Combined	5,194/12,189	1.18	(1.11–1.25)	3.63×10−11	0.09	43.85
6p21.33	HLA-B	rs2523607	31,322,790	A	T	0.120	Stage 1	2,661/6,221	1.45	(1.29–1.64)	7.10×10−10
						0.123	Stage 2	1,195/1,444	1.14	(0.96–1.35)	0.14
						0.109	Stage 3d	1,114/1,102	1.25	(1.04–1.51)	0.019
							Combined	4,970/8767	1.32	(1.21–1.44)	2.40×10−10	0.26	21.63
2p23.3	NCOA1	rs79480871	24694472	T	C	0.076	Stage 1	2,660/6,220	1.35	(1.17–1.55)	3.51×10−5
						0.057	Stage 2	1,195/1,443	1.56	(1.22–1.99)	0.00037
						0.063	Stage 3	1,344/4,524	1.19	(0.98–1.46)	0.084
							Combined	5,199/12,187	1.34	(1.21–1.49)	4.23×10−8	0.15	34.59

Position according to human reference NCBI37/hg19;

Allele associated with an increased risk of DLBCL;

Risk allele frequency in controls;

Not genotyped in NCI Replication study.

Figure 1

Association results, recombination hot-spots, and linkage disequilibrium (LD) plots for the regions newly associated with diffuse large B-cell lymphoma (DLBCL)

(a–d) Top, association results of GWAS data from stage 1 DLBCL-GWAS (grey diamonds) and combined data of stages 1–3 (red diamond) are shown in the top panels with −log10(P) values (left y axis). Overlaid are the likelihood ratio statistics (right y axis) to estimate putative recombination hotspots across the region on the basis of 5 unique sets of 100 randomly selected control samples. Bottom, LD heatmap based on r2values from combined control populations for all SNPs included in the GWAS. Shown are results for 6p25.3 (a), 8q24.21 (b), 2p23.36 (c), and p21.33 (d) regions.

Within the HLA region, rs2523607 (P=3.35×10−9) was carried forward for replication in stage 3. This SNP, localized at 6p21 in HLA-B, reached a combined P=2.40×10−10 in a meta-analysis of all three stages (Table 1; Figure 1d). To further evaluate the association of HLA variants with DLBCL risk, we imputed classical HLA alleles at six loci (HLA-A, B, C, DRB1, DQA1, and DQB1) in the four GWAS datasets from stages 1–2 and conducted a meta-analysis (Online Methods). The imputation accuracy of HLA types was high (>95.2%) when compared to HLA sequencing (four-digit resolution) previously performed on a subset of the NCI samples[14] scanned as part of this study in stage 1 (Online Methods). Of all SNPs and classical HLA alleles tested across the MHC, only the SNP rs2523607 (OR=1.34, P=3.3×10−9 in stages 1 and 2) and the classical allele HLA-B*08:01 (OR=1.30, P=3.16×10−8 in stages 1 and 2) reached genome-wide significance (Supplementary Table 9). These markers were in very high LD (r2=0.91), and after adjusting for the effect of HLA-B*08:01, the association of rs2523607 was greatly weakened (P=5.5×10−3).

To gain additional insight into potential biological mechanisms, expression quantitative trait loci (eQTL) analyses were performed in two datasets consisting of lymphoblastoid cell lines (Online Methods). In one of the datasets, significant associations were observed for rs116446171 with HIST1H3F and rs2523607 with HCG27 (Supplementary Table 10), while in the other dataset significant associations (FDR<0.05) were observed for rs2523607 (using rs3130923 as a proxy, r2=0.94) with LY6G6E, FLOT1, and RNF5 (Supplementary Table 11); no associations were observed for the other DLBCL-associated loci.

To explore plausible mechanisms for the non-coding variants identified in our GWAS, the sentinel SNPs and those in high linkage disequilibrium (r2≥0.8) in Europeans in the 1000 Genomes Project were analyzed using HaploReg v2[15] (Online Methods; Supplementary Table 12). In addition, B-cell specific chromatin dynamics were assessed in a lymphoblastoid cell line (GM12878) using ChroMoS[16], which utilizes the pre-computed chromatin state data for 9 cell lines (including GM12878)[17]. Of the 173 SNPs queried, 61 had information for GM12878 (Supplementary Figure 5), and 3 SNPs were identified as active or weak promoters only in GM12878, while 22 SNPs were identified as strong or weak enhancers in GM12878. In the other 8 cell lines, these regions were mostly defined as neutral, weakly transcribed or polycomb repressed. These results suggest that some of our SNPs are within regions of active chromatin state predominantly within B cells and have a role in the B-cell cis-regulatory network. These results are consistent with growing evidence that disease variants from GWAS are more likely to map to active chromatin sites than neutral sites, as was shown recently for systemic lupus erythematosus[17]. HaploReg showed that the majority of DLBCL-related SNPs were observed in regions of DNAse hypersensitivity common across multiple cell lines (e.g., rs116446171, rs2523607, rs13255292, rs4733601 near EXOC2, HLA-B, PVT1 or 7SK) whereas rs147193201 was specific to B-cells. The preponderance of DNAse hypersensitivity points to the existence of motifs, such as enhancers, silencers, promoters, insulators and other control elements of gene regulation. The proteins bound at these sites are known transcription factors such as NF-κB, c-MYC, GATA2 or genes that regulate transcription such as POL24H8, USF1 or POL2. These suggested mechanisms of action will require laboratory follow-up.

The susceptibility locus at 6p25.3 (rs116446171) maps near a plausible DLBCL candidate gene, EXOC2 (exocyst complex component 2), which is part of a large multiprotein complex responsible for vesicle trafficking and maintenance and intercellular transfer of viral proteins and virions[18]. EXOC2 functions at the interface between host defense and cell death regulation[19]. EXOC2 interacts with Ral proteins, and the Ral-exocyst regulatory node has a crucial role in the maintenance of epithelial cell polarity, cell motility and cytokinesis[20,21], and in proliferation and metastasis[20,22]. It is notable that IRF4 is centromeric to EXOC2 and genetic variation in this region has been linked with chronic lymphocytic leukemia (CLL) risk[23,24], and nominally to DLBCL risk[25]. However, rs116446171 was not in LD with the IRF4 CLL GWAS SNP rs872071[23].

Two 8q24.21 variants (Figure 1b), rs13255292 and rs4736601 positioned at chr8:129.07Mb and chr8:129.26Mb, respectively, are approximately 1Mb telomeric to the 8q24 region linked with multiple cancers[26], including CLL[27]. Both variants are in close proximity to PVT1, which is a non-coding RNA implicated in the MYC activation. Notably, a variant at 8q24.21 (rs2019960) has been linked to HL[28], but the pair-wise r2 values of this SNP with both of our SNPs were low (r2<0.02). The close proximity of PVT1 and the MYC oncogene, which is known to be deregulated in Burkitt lymphoma[29,30] and some DLBCLs [31,32], suggests that germline variation in this region could also contribute to DLBCL risk.

The susceptibility locus at 2p23.3 (rs79480871) maps near NCOA1, nuclear receptor coactivator 1 and ITSN2, intersectin 2. The former gene acts as a transcriptional coactivator for steroid and nuclear hormone receptors and is a member of the p160/steroid receptor coactivator (SRC) family[33], while the latter gene encodes a protein that is a member of a family of proteins involved in clathrin-mediated endocytosis[34] and may also augment the induction of T-cell receptor endocytosis[35]. However, our bioinformatics analysis did not identify a clear link to genes in this region, supporting the need to refine this signal in future work.

Through imputation with SNP2HLA,[36] our strongest associations in the HLA region were with the HLA-B SNP rs2523607 and HLA-B*08:01, which are in very high LD, and based on our available sample size we cannot definitively rule out an orthogonal effect of rs2523607 in favor of HLA-B*08:01. HLA-B encodes the HLA class I heavy chain paralogue, which heterodimerizes with a light chain (β2 microglobulin) to play a central role in presenting intracellularly processed self or foreign antigens to CD8+ cytotoxic T lymphocytes. Class I molecules have been linked to a variety of immune-mediated diseases and cancers including HL, follicular lymphoma, DLBCL[7,14,37,38], and more recently marginal zone lymphoma (Vijai, submitted). Our results strongly suggest HLA-B*08:01 as the primary MHC association with DLBCL risk. This classical allele is carried by the so-called ancestral 8.1 haplotype associated with other complex diseases (e.g., type I diabetes).[39] Classical alleles of other HLA loci may also be involved (including those on the 8.1 haplotype), but larger sample sizes will be required to evaluate this question.

Our study represents the largest DLBCL GWAS in individuals of European descent. We did not observe a notable signal for a locus previously reported for DLBCL on 3q27 in East Asia[10], rs6773854 (reported as OR=1.47, P=1.14×10−11), which was based on a discovery set of 253 B-cell NHL cases (148 DLBCLs). Although our current study had a similar MAF of 0.22 among controls, we observed an OR=1.06 and a P-value of 0.81 for this SNP (Supplementary Table 13), suggesting that the reported marker may not be correlated with the functional susceptibility allele in Europeans. Of the two suggestive loci (P<5×10−7) reported in the literature[8,40], we did not observe an association for rs751837 with DLBCL (OR=0.97, P=0.46), identified in a small Japanese GWAS (OR=3.51, P=3.3×10−7)[40], but we did observe a consistent albeit attenuated association for rs10484561 (OR=1.18, P=1.5×10−4) which was initially reported on a subset of the studies in stage 1 (OR=1.36, P=1.46×10−7)[8]. Previously, an InterLymph study of ~1,800 DLBCLs and ~6,500 controls reported a strong signal for a dinucleotide haplotype in the LTA/TNF locus (LTA 252A>G/TNF–308G>A) at 6p21.3 (OR=1.31, P=2.9 x10−8)[5]. Although nearly all of the cases from the previous publication were included in our current GWAS, the signal we observed overall was weaker (OR=1.15, P=8.5×10−4). The attenuation was not explained by study design (case-control, cohort) or adjustment for population substructure (data not shown), but could be due to population sampling differences, heterogeneity, or chance.

To explore the heritability of DLBCL, we estimated the contribution of all common SNPs to the variance explained by fitting all genotyped autosomal SNPs simultaneously using the method proposed by Yang et al[41] in the Stage 1 dataset. We estimated that common SNPs, including but not limited to the loci discovered in this study, explain approximately 16% of the variance for DLBCL overall.

In summary, our findings represent an important step in defining the contribution of common genetic variants to risk for DLBCL. Our findings are notable because we have newly defined associations of several regions with susceptibility to DLBCL, and these regions harbor plausible candidate genes for further investigation. Further studies are required to discover additional common susceptibility loci as well as functional analyses that can explain the biological underpinnings of these new susceptibility loci.

ONLINE METHODS

Stage 1: DLBCL-GWAS

As part of a larger initiative, we conducted a genome-wide association study (GWAS) of diffuse large B-cell lymphoma (DLBCL) using cases and controls of European descent from 22 studies of non-Hodgkin lymphoma (NHL) (Supplementary Table 1), including nine prospective cohort studies, eight population-based case-control studies, and five clinic or hospital-based case-control studies. All studies were approved by their respective Institutional Review Boards, and informed consent was obtained for all participants. Cases were ascertained from cancer registries, clinics or hospitals, or through self-report verified by medical and pathology reports. To determine NHL subtype, phenotype data for all NHL cases were harmonized to the hierarchical classification proposed by the InterLymph Pathology Working Group[42,43] based on the World Health Organization (WHO) classification[44].

All DLBCL cases with sufficient DNA (n=2,878) and a subset of controls (n=2,854) frequency matched by age, sex, and study to the entire group of NHL cases, along with 4% quality control duplicates, were genotyped on the Illumina OmniExpress at the NCI Cancer Genomic Research Laboratory (CGR). Genotypes were called using Illumina GenomeStudio software, and quality control duplicates showed >99% concordance. Monomorphic SNPs and SNPs with a call rate of <95% were excluded. Samples with a call rate of ≤93%, mean heterozygosity <0.25 or >0.33 based on the autosomal SNPs, or gender discordance (>5% heterozygosity on X chromosome for males and <20% heterozygosity on the X chromosome for females) were excluded. Furthermore, unexpected duplicates (>99.9% concordance) and first-degree relatives based on identity by descent (IBD) sharing with Pi-hat>0.40 were excluded. Ancestry was assessed using the GLU struct.admix module based on the method by Pritchard et al.[45] and participants with <80% European ancestry were excluded (Supplementary Figure 6). After exclusions, 2,661 (92.5%) cases and 2,685 (94.1%) controls remained (Supplementary Table 2). Genotype data previously generated on the Illumina Omni2.5 from an additional 3,536 controls from three of the studies (ATBC, CPSII, and PLCO) were also included[11], resulting in a total of 2,661 cases and 6,221 controls for the stage 1 analysis. Of these additional 3,536 controls, 703 (~235 from each study) were selected to be representative of their cohort and cancer-free[11], while the remainder were cancer-free controls from an unpublished study of prostate cancer in the PLCO. SNPs with call rate <95%, with Hardy-Weinberg equilibrium P-value<1×10−6, or with a minor allele frequency <1% were excluded from analysis, leaving 611,844 SNPs for analysis (Supplementary Table 4). To evaluate population substructure, a principal components analysis (PCA) was performed using the Genotyping Library and Utilities (GLU), version 1.0, struct.pca module, which is similar to EIGENSTRAT[46]. Plots of the first five principal components are shown in Supplementary Figure 7. Association testing was conducted assuming a log-additive genetic model, adjusting for age, sex, and four significant principal components. All data analysis and management was conducted using GLU.

Stage 2: In Silico Analysis of Three Independent DLBCL GWAS

Three independent DLBCL GWAS provided genotyping data for a meta-analysis, (Supplementary Table 1), which included data generated with the following commercial, SNP microarrays: Illumina HumanHap 660W for Mayo (393 DLBCL and 172 controls), HumanCNV370-Duo for UCSF[7] (254 DLBCLs and 748 controls), and HumanHap 610K for GELA (549 cases). In all studies, subjects with a genotyping call rate <95%, duplicates, related individuals, and SNPs with a call rate <95% were removed prior to imputation (Supplementary Table 4). The GELA study was conducted on cases only; controls were drawn from a pool of 928 individuals from the French component of the EPIC cohort, who were previously scanned on Illumina HumanHap 660W or 610K[47,48]. We subsequently chose a subset of 525 individuals with matched ancestry as determined from the principal components analysis. In total, there were 1,196 cases and 1,445 controls in stage 2.

Imputation was conducted separately for each study in stages 1 and 2 using IMPUTE2[13] and the 1000 Genomes Project version 3[12]. The imputation analysis was restricted to common SNPs (cut-off MAF>0.01 with imputation accuracy INFO score >0.3).

Association testing was conducted for each study using SNPTEST version 2, adjusting for age, sex, and any significant principal components. We evaluated the top 10 eigenvectors for the GELA, Mayo and UCSF studies, respectively, in each baseline risk model adjusting for both age and gender. Based on the significance level (P<0.05) of the regression coefficient for eigenvectors, we chose to adjust for three eigenvectors (EV1, EV7 and EV8) for GELA in the final association model, while no eigenvectors met criteria for adjustment of either the Mayo or UCSF studies.

All meta-analyses were performed using the fixed effects inverse variance method based on the beta estimates and standard errors from each study.

Stage 3: Replication studies and technical validation

In stage 3, eight SNPs in the most promising loci outside of the HLA region and one SNP from the HLA region (Supplementary Table 7) were taken forward for de novo replication in an additional 1359 cases and 4557 controls from four studies (Supplementary Table 1), except for rs2523607, which was not genotyped in one of the studies (NCI replication). Genotyping was conducted using custom TaqMan genotyping assays (Applied Biosystems) at the NCI Cancer Genomics Research Laboratory. Each assay was optimized and validated with 270 HapMap samples and additional CEPH samples (SNP500Cancer), and these samples were used as genotyping controls for clustering and reproducibility. All validated assays had 99% or higher concordance with HapMap and completion with control DNA was >97%. Blind duplicates from stage 3 samples (64 pairs; ~3%) yielded 100% concordance.

In technical validation, we observed a high correlation of genotyping calls from the OmniExpress microarray with confirmatory TaqMan assays in 455 stage 1 duplicate samples for two genotyped (rs13255292, r2=1.00; rs4733601, r2=1.00) and four imputed (rs116446171, r2=0.92; rs2523607, r2=0.99; rs2681416, r2=1.00; rs79480871, r2=0.94) SNPs. We also observed a high correlation of genotyping calls from the Illumina HumanHap 660W microarray with confirmatory TaqMan assays in stage 2 duplicate samples from the Mayo study (N=165) for two genotyped (rs13255292, r2=1.00; rs4733601, r2=1.00) and four imputed (rs116446171, r2=1.00; rs2523607, r2=1.00; rs79480871, r2=0.85; rs79464052, r2=0.95) SNPs.

HLA imputations and analysis

We imputed dense SNPs as well as classical HLA alleles (A, B, C, DRB1, DQA1, DQB1) and coding variants across the HLA region (chr6:29.5–33.2Mb, hg19) in the stage 1 (NCI) and stage 2 (MAYO, USCF2 and GELA/EPIC) studies using SNP2HLA[36]. The imputation was based on a reference panel from the Type 1 Diabetes Genetics Consortium (T1DGC), and consisted of genotypes from 5,225 individuals of European descent who were typed for HLA-A, B, C, DQA1, DQB1, DRB1, DPA1, DPB1 4 digit alleles. To assess imputation accuracy, we compared the imputed HLA alleles to HLA sequencing data (to 4 digits) available on a subset of samples from the NCI GWAS[14], and found high concordance rates for HLA-A (97.3%), B (98.5%), C (98.1%) and DRB1 (97.5%). Due to the limited number of SNPs (N=7,253) in the T1DGC reference set, imputation of HLA SNPs was conducted with IMPUTE2 and the 1000 Genomes reference set as described above. A total of 68,488 SNPs, 201 classical HLA alleles (two- and four-digit resolution) and 1,038 AA markers including 103 AA positions that were ‘multi-allelic’ with three to six different residues present at each position, were successfully imputed (info score >0.3 for SNPs or r2>0.3 for alleles and AAs) and available for analysis. Multi-allelic markers were analyzed as binary markers (e.g., allele present or absent) and using a global test, and a meta-analysis was conducted where we tested SNPs, HLA alleles and AAs across the HLA region for association with DLBCL using PLINK[49] or SNPTEST as described above.

Expression quantitative trait loci (eQTL) analysis

To evaluate the effect of our top loci (and SNPs in LD based on r2>0.8 in HapMap-CEU release 28) on gene expression, we conducted an eQTL analysis on lymphoblastoid cell lines using two independent datasets: Childhood asthma[50] and HapMap[51]. For the childhood asthma dataset [50], peripheral blood lymphocytes were transformed into lymphoblastoid cell lines for 830 parents and offspring from 206 families of European ancestry. Using extracted RNA, gene expression was assessed with the Affymetrix HG-U133 Plus 2.0 chip. Genotyping was conducted using the Illumina Human1M Beadchip and Illumina HumanHap300K Beadchip, and imputation was performed using data from the 1kGP. All SNPs selected for replication were tested for cis associations (defined as gene transcripts within 1 Mb), assuming an additive genetic model, adjusting for non-genetic effects in the gene expression value. To gain insight into the relative importance of associations with our SNPs compared to other SNPs in the region, we also conducted conditional analyses, in which both the DLBCL SNP and the most significant SNP for the particular gene transcript (i.e., peak SNP) were included in the same model. Only cis associations that reached P<6.8×10−5, which corresponds to a false-discovery rate (FDR) of 1%, are reported (Supplementary Table 10).

The HapMap dataset consisted of a publicly available RNAseq dataset[51] from transformed lymphoblastoid cell lines from 41 CEPH Utah residents with ancestry from northern and western Europe (HapMap-CEU), samples available from the Gene Expression Omnibus (GEO) repository (http://www.ncbi.nlm.nih.gov/geo) under accession number GSE16921. Genotyping data for the same HapMap-CEU individuals were directly downloaded from HapMap (www.hapmap.org). Since rs2523607, rs79480871 and rs116446171 were not genotyped in HapMap, we selected rs3130923, rs6746301 and rs7762424 as respective proxies, as they were the strongest linked SNPs available in HapMap (r2=0.94, 0.69 and 0.54 in 1kGP-CEU, respectively). Correlation between expression and genotype for each SNP-probe pair was tested using the Spearman’s rank correlation test with t-distribution approximation and were estimated with respect to the minor allele in HapMap-CEU. P-values were adjusted using the Benjamini-Hochberg false-discovery rate (FDR) correction and eQTLs were considered significant at an FDR≤0.05 (Supplementary Table 11).

Bioinformatics: ENCODE and Chromatin State Dynamics

Using 1000 Genomes data, we identified SNPs with r2≥0.8 with our sentinel SNP that were reported to be non-synonymous or nonsense variants. We utilized HaploReg v2[15], which is a tool for exploring non-coding functional annotation using ENCODE data, to evaluate the genome surrounding our SNPs (Supplementary Table 12). To assess chromatin state dynamics, we used Chromos[16], which has pre-computed data from ENCODE on 9 cell types based on Chip-Seq analyses[17]. These pre-computed data have genome-segmentation performed using multivariate hidden Markov-model to reduce the combinatorial space to a set of interpretable chromatin states. The output from Chromos separates data into 15 chromatin states corresponding to repressed, poised and active promoters, strong and weak enhancers, putative insulators, transcribed regions, and large-scale repressed and inactive domains (Supplementary Figure 5).

Heritability analyses

To estimate the contribution of all common SNPs to the variance explained, we used the method proposed by Yang et al[41], which was extended to dichotomous traits[52] and implemented in the Genome-wide Complex Trait Analysis (GCTA) software[53]. The genetic similarity matrix was estimated from our stage 1 data using all genotyped autosomal SNPs with a minor allele frequency >0.01. We used restricted maximum likelihood (REML), the default option for GCTA, to fit the appropriate variance components model that included the top 10 eigenvectors as covariates. The final estimate of heritability on the underlying liability scale assumed that the lifetime risk of DLBCL was 0.0074[54].

Estimate of recombination hotspots

To identify recombination hotspots in the region we used SequenceLDhot[55], a program that uses the approximate marginal likelihood method[56] and calculates likelihood ratio statistics at a set of possible hotspots. We tested five unique sets of 100 control samples. PHASE v2.1 program was used to calculate background recombination rates[57,58] and LD heatmap was visualized in r2 using snp.plotter program[59].