{'fraction_non_alphanumeric': 0.043481745681881935, 'fraction_numerical': 0.05113613593555885, 'mean_word_length': 4.911869225302061, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.

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 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).
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 (r 2 =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, r 2 =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 (r 2 ≥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 r 2 values of this SNP with both of our SNPs were low (r 2 <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 socalled 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.

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.

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 r 2 >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 r 2 >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 (r 2 =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. Pvalues 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 r 2 ≥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 .

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.

Author Manuscript
Cerhan et al.