Age-related cancer mutations associated with clonal hematopoietic expansion

Cancer types and sample characteristics

We searched for variants present in the blood normal controls across 2,728 cancer patients (Supplementary Table 1a) from 11 diverse cancer types: breast adenocarcinoma (BRCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), brain low grade glioma (LGG), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), ovarian carcinoma (OV), prostate adenocarcinoma (PRAD), stomach adenocarcinoma (STAD), and uterine corpus endometrioid carcinoma (UCEC). The numbers of cases from each tumor type range from 57 (KIRC) to 673 (BRCA) and are listed in Supplementary Table 1b. Patients were diagnosed between 10–90 years (mean 59.5 ± 13.1 years) and 22.1% were deceased at the time of TCGA sample procurement (Supplementary Table 1b). TCGA collects clinical data regarding diagnosis and prior treatment of neoplasms during the sample submission process. To ensure that our dataset was comprised only of individuals having first-time primary cancers and having had no treatment with radiation and/or chemotherapy, we excluded those having reported histories of these events as identified at https://tcga-data.nci.nih.gov/annotations/ and all clinical data (July 30^th, 2014). However, five patients with synchronous tumors not associated with blood were included, since these synchronous tumors would be unlikely to affect variant analysis in corresponding blood samples.

Variant calling and filtering strategies

Variants in the 2,728 blood normal controls were identified with VarScan (single nucleotide variant and indel), GATK (single nucleotide variant and indel), and Pindel (indel) (see Methods). False-positive filters were subsequently applied prior to downstream analysis and interpretation (see Methods). Out of the 49,317,027 variants (previously reported OV counts¹⁰ were not included here) that passed false positive filters, 1,622,485 with minor allele frequency of <1% in the 1000 Genomes reference and in each cancer cohort were retained for further analysis; this consists of 1,025,632 missense, 529,505 synonymous, 19,663 nonsense, 10,976 splice site, 926 nonstop/readthrough, 20,275 frameshift indels, and 15,508 in frame indels (Supplementary Table 1c). We used a stringent filtering strategy described previously¹¹ for standardizing specificity across the Pan-Cancer somatic variant calls for available matched tumor samples (Supplementary Table 2).

Variants contributing to hematopoietic clonal expansion

The collection of both tumor and matched blood normal exome data by TCGA provides a unique comparative resource for identifying those somatic variants in blood that contribute to clonal expansion. We set out to identify both rare truncation variants (RTV), i.e. those having <1% MAF in both the 1000 Genomes collection and the cohort data, and variants overlapping with recurrent somatic mutations (also called Known Hotspot Variants: KHV) found in the analysis of 12 TCGA cancer types (see Methods and Supplementary Table 3a and 3b). Subsequently, 1,598 RTVs in 556 cancer-associated genes based on several recent studies^11–15 (see Methods and Supplementary Table 4) were manually reviewed to remove false positive calls and 136 KHVs in the same set of cancer genes were identified. The resulting list was further filtered to remove polymorphisms present in 1000 Genomes (see Methods) and having greater than 0.1% MAF as reported in the current Exome Variant Server data release (ESP6500SI-V2, http://evs.gs.washington.edu/EVS/) from 6,503 samples drawn from multiple NHLBI Exome Sequencing Project (ESP) cohorts. We focused on those mutations found in blood normal samples, but not present or present only at very low levels in either the tumor samples or tumor adjacent normal samples, as this pattern is highly suggestive of somatic mutation in HSPCs introduced by the clonal expansion process. Inflammatory lymphocytes/macrophages/neutrophils will infiltrate different tumors to different extents. Therefore, the hematopoietic mutations do not have to be completely absent in the tumor sample.

Our analysis of RTVs and KHVs in 556 selected cancer-associated genes identified 70 blood-specific mutations in 58 individuals. We further performed comparative analysis of blood versus tumor samples for these 58 individuals, with the goal of detecting all blood-specific nonsynonymous mutations in the 556 cancer-associated genes; this analysis identified seven additional events (those that were likely loss of heterozygosity related to copy number alterations in the tumor were not included), yielding a final list of 77 blood-specific mutations in 58 cases (Table 1 and Supplementary Table 5a). For five of those 58 cases that also had adjacent normal tissue analyzed, blood variants were absent in the adjacent normal tissue. Interestingly, among the 31 genes harboring these events, 19 have already been linked to hematological malignancies (Fig. 1a). More strikingly, 64 out of the total 77 events (83%) were in these 19 genes, examples as follows: DNMT3A¹⁶ (18 cases), TET2¹⁷ (9 cases), JAK2¹⁸ (8 cases), ASXL1^19,20 (6 cases), TP53²¹ (4 cases), SF3B1²² (2 cases), BCORL1²³ (2 case), ASXL2²⁴ (1 case), and SH2B3²⁵ (1 case) (Table 1 and Fig. 1a). The overall frequency of blood-specific mutations increased with age (logistic regression analysis, P = 2.38e–08), for example, 0.9% of the cases were in their 40s, 1.0% in their 50s, 1.8% in their 60s, 5.3% in their 70s, and 6.1% in their 80s (Fig. 1b). The blood-specific mutations were found in all 11 cancer types (Fig. 1b). Frequencies of individual TET2, DNMT3A, ASXL1, and SF3B1 mutations also show association with age, (all FDR values <0.034, logistic regression). Interestingly, TET2, ASXL1, and SF3B1 mutations were predominantly found in the oldest age groups (70s and 80s) and DNMT3A mutations in 60s to 80s. In contrast, JAK2 mutations, which did not achieve significance, trended in both younger (40s and 50s) and older age groups (70s) (Fig. 1c).

The average age of the 54 individuals with blood-specific mutations in the 19 leukemia/lymphoma-associated genes was 70.0 ± 9.9 years, significantly higher than that of the larger TCGA cohort used in this study (difference of means test, P = 3.4e–10) (Fig. 2a). Notably, the six individuals having two blood-specific mutations in the nine recurrently mutated genes are relatively older, with ages of 64 (DNMT3A: R882H; 36% VAF; TET2: H863fs; 12% VAF), 72 (JAK2: V617F; 73% VAF and TET2: T229fs; 19% VAF), 75 (DNMT3A: Y584fs; 38% VAF and TET2: Q764fs; 33% VAF), 76 (JAK2: V617F; 42% VAF and ASXL1: R548fs; 35% VAF), 83 years (TET2: F381fs; 50% VAF and TET2: Q888*; 20% VAF), and unknown age (BCORL1: G883E; 17% VAF and TP53: Q136*; 18% VAF), respectively (Table 1 and Supplementary Table 5a). We also compared the distribution of variant allele fractions for these 64 events versus inherited sites identified in the same sample set; we observed a clear shift towards lower VAFs in the blood-specific sites (Fig. 2b), suggesting the majority are present in only a fraction of the blood cells.

Although GNAS mutations have been found in leukemia²⁶, activating gain-of-function mutations in GNAS are best known for their involvement in polyostotic fibrous dysplasia and McCune-Albright syndrome²⁷. Interestingly, previous studies showed that activating mosaic GNAS mutations could affect various tissue types and the non-mosaic state for activating GNAS mutations may be lethal for the embryo^28–30. This is consistent with our finding of mosaic GNAS R202H in transcript ENST00000354359 (also known as R201H in transcript ENSP00000360126) in the three blood samples of TCGA cases (11.5%, 14.4%, 21.4% VAFs, respectively). It is also worth noting that two blood-specific truncation mutations were detected in PPM1D, a gene recently found to be associated with breast and ovarian predisposition with mosaic signatures³¹, but not with hematological malignancies.

Due to the fact that the blood-specific variants were present in only a fraction of the blood cells, we postulated that certain low level variants associated with clonal expansion were not captured by the variant detection tools. Therefore, we performed read count-based analysis for a set of known hotspot variants, including R882 in DNMT3A, R132 in IDH1, R172/R140 in IDH2, V617 in JAK2, K700 in SF3B1, and S34 in U2AF1 for the entire TCGA cohort. We compared the distributions of VAFs between tumor and blood normal samples and found that the strongest differences were at R882 in DNMT3A and V617 in JAK2 (Fig. 3 and Table 1). We devised a statistical procedure (see Methods) to identify additional “low VAF” sites significantly above the background error rate. Our analysis identified 14 blood samples having low-level hot spot variants (12 of them are not part of the 58 cases identified), including eight in DNMT3A, four in JAK2, one in both SF3B1 and IDH2, but none in IDH1 or U2AF1. The average age for these 14 cases is 64.0 ± 14.9 years, again higher than the entire TCGA cohort studied (Supplementary Table 5b). We performed deep sequencing for five selected low-level hot spot variants (two R882C and one R882H sample sites in DNMT3A, one V617F in JAK2, and one K700E in SF3B1) to evaluate our detection approach. We achieved more than 450,000× average coverage for each site and all of the five variants were validated as bona fide low-level blood specific mutations (Supplementary Table 5c and see Methods). By including low VAF events, the overall frequencies of blood-specific mutations in different age groups are: 1.2%, 1.3%, 2.2%, 6.1%, and 6.8% in their 40s, 50s, 60s, 70s, and 80s, respectively (Supplementary Fig. 1). Collectively, we estimate that approximately 2% and 3.5% of individuals over the ages of 40 and 60, respectively, and without overt hematological malignancies carry blood-specific mutations that are associated with hematological malignancies.

Known Hotspot Variants in NHLBI exome sequencing project control cohort

To confirm that these mutations are also present in an independent set of normal blood samples, we examined the NHLBI exome sequencing project (ESP) control cohort. We searched for RTVs and KHVs in 6,503 ESP samples, focusing on four AML-associated genes discovered in the TCGA set (DNMT3A, TET2, JAK2, and ASXL1). When we applied a 0.1% MAF threshold in ESP (to focus on rare variants and to prevent potential false positives), we identified 8 RTVs and 13 KHVs in DNMT3A, 13 RTVs in TET2, and 7 RTVs in ASXL1, and 3 KHVs in JAK2 (Supplementary Table 6). For a subset of ESP samples (n = 557 WHISP), we re-aligned reads and performed variant analysis using the same process applied for the TCGA cohort. This allowed us to confirm RTVs and KHVs detected by the pipeline, and also to detect low VAF KHVs potentially missed by variant detection tools. Our pipeline detected one RTV each in DNMT3A, TET2, and ASXL1, four DNMT3A R882H mutations, and two JAK2 V617F mutations. Careful analysis of recurrent hotspot sites (described in the previous section) using the 557 WHISP samples identified an additional three DNMT3A R882, one JAK2 V617, one GNAS R202 and two SF3B1 K700 variants with VAFs ranging from 2% to 10% (Fig. 3 and Supplementary Table 7). Based on these analyses, over 2% of the 557 WHISP cases contain mutations in these five selected genes. However, sequencing of matched non-blood samples from these cases would be required to prove that they are truly blood-specific mutations.

Mutations in TCGA blood samples and patients with hematological malignancies

We next compared the mutation frequencies in 25 genes frequently mutated in at least one of the five following cohorts: TCGA 58 blood samples, 151 myeloproliferative neoplasm (MPN) cases reported by Nangalia et al.³², 150 myelodysplastic syndrome (MDS) cases reported by Walter et al.³³, 160 chronic lymphocytic leukemia (CLL) by Landau et al.³⁴, and 200 TCGA AML cases²¹. DNMT3A, TET2, ASXL1, TP53, and SF3B1 were found to be consistently mutated in at least four groups, while JAK2 was more specifically mutated in 58 TCGA blood samples and MPN patients (Fig. 4a). No mutations were found in TCGA blood samples in genes such as IDH1, NRAS, RUNX1, and PHF6, which are significantly mutated in AML and frequently mutated in MPN and/or MDS. Several genes showed cohort specificity: GNAS was mutated only in TCGA blood samples and CEBPA, WT1, PTPN11, KIT, SMC1A, and SMC3 were preferably mutated in the AML cohort. We reasoned that common mutations among the cohorts (e.g., in DNMT3A, TET2, ASXL1, SF3B1, JAK2, and TP53) may be relevant for initiating HSPC clonal expansion, and are also likely to be early, initiation events for hematological malignancies, such as MPN, MDS, CLL, and AML. On the other hand, genes specific for MPN, MDS, and/or AML (e.g., NRAS, RUNX1, NPM1, and FLT3) are more likely to be subsequent, cooperating mutations that are involved in the progression of these diseases. These observations also show both distinct and common connections among these five cohort groups, suggesting that the TCGA cohort consists of a combination of precursor mutations that may sometimes evolve to MPN, MDS, CLL, and/or AML, although subsequent, collaborating mutations are clearly required (Supplementary Fig. 2). Finally, we compared the average number of mutations among these 4 cohorts (not including MDS) in 556 selected cancer-associated genes; this showed TCGA blood samples having the fewest mutations, MPN and CLL cohorts having intermediate numbers of mutations, and AML patients harboring the highest number of mutations (Fig. 4b).

Variant calling and annotation strategies

Exome sequencing data were aligned to NCBI Build 37 of the human reference using BWA v0.5.9 and de-duplicated using Picard 1.29. Single nucleotide variants were identified by Varscan (version 2.2.6: –min-var-freq 0.1 –p-value 0.1 –min-coverage 8 –map-quality 10), and GATK (revision 5336: –T UnifiedGenotyper –R GRCh37-lite –et NO_ET –l INFO –U ALL –validation_strictness SILENT). Indels were identified using Varscan (version 2.2.9: –min-coverage 3 –min-var-freq 0.2 –p-value 0.1 –strand-filter 1 –map-quality 10), GATK (revision5336: –T IndelGenotyperV2 –R GRCh37-lite –window_size 300 –et NO_ET –U ALL), as well as Pindel (version 0.2.4×, May 8, 2013: –window-size 1). For Pindel analysis, we preset the insertion size to 500 if this information is not provided in the BAM header. SNVs were based on the union of GATK and VarScan. They were subsequently processed through our in-house false-positive filter (–min-homopolymer 10). We required that indels were called by at least 2 out of 3 callers (GATK, Varscan, Pindel). In addition we also included Pindel unique calls (at least 30× coverage and 20% VAF). All combined indels were then processed through our false-positive-filter (–min-homopolymer 10 –min-var-freq 0.2 –min-var-count 6). We then applied additional annotation and minor allele frequency filters as previously reported¹⁰.

Variant transcript annotation is based on all human transcripts obtained from Ensembl Release 70_37_v5. The reference alleles and positions were derived from the sequence and coordinates of GRCh37. All transcripts were annotated and a single representative was selected for each somatic mutation based on the significance of the predicted functional effect of each mutation, ordered from most significant to least significant as follows: nonsense, frameshift, splice site, in frame, missense, no stop (nonstop/readthrough), synonymous, and RNA. Splice site mutations were restricted to substitutions, deletions, or insertions overlapping the 2bp intronic sequence defined as the canonical splice donor or splice acceptor. RNA mutations were restricted solely to transcripts without an annotated open reading frame. Mutations affecting 3′UTR, 5′UTR, intronic sequence, and intergenic sequence were discarded for the purposes of downstream analysis.

Recurrent somatic mutations in 12 cancer types

We collected extensively filtered somatic variants in 3,355 TCGA samples from 12 cancer types (Supplementary Table 3a), and selected recurrent mutations appearing more than once at a given genomic position (Supplementary Table 3b).

Compiling cancer-associated gene list

A total of 556 candidate cancer-associated genes were compiled using nine sources, including recently published large-scale cancer studies, publicly available screening panels, and unpublished preliminary analysis of publicly available data sources. The 204 genes shared across at least two of the nine sources were retained and a literature search was conducted to identify evidence supporting inclusion of any remaining unique genes. A subset of 518 genes originated from recent publications, including 294 genes from Frampton et al.¹⁵, 125 genes from Kandoth et al.¹¹, 212 genes from Lawrence et al.¹³, 194 genes from Pritchard et al.¹⁴, and 124 genes from Vogelstein et al.¹² Thirty-nine additional genes were included based on the analysis of driver mutations in 20 TCGA cancer types (unpublished), recommendations in accordance with the standards and guidelines of the American College of Genetics and Genomics⁴², and 18 novel cancer driver genes identified in recently published large-scale studies (Supplementary Table 4).

Readcount analysis and statistical approaches for identifying significant low level variants

Read counts for variants were determined using our internally developed tool bam-readcount (https://github.com/genome/bam-readcount). For sites to be included in the downstream statistical test, we require greater than 30× coverage for both blood normal and matched tumor samples. We assessed the false-discovery rate (FDR) using a two-stage process. The first step is a purpose-specific pre-filter to eliminate candidates that can be confidently identified as having originated from non-mosaic samples, as these identically fail the inclusion criterion for significance testing. It targets deeply covered sites whose apparent variant reads actually represent base-calling errors. Take the sample space for blood classification as consisting of two mutually exclusive, collectively exhaustive (MECE) statuses, “normal”, S=N, and “mosaic”, S=M, and let the candidate blood site data, D, consist of T spanning reads, of which V and T-V report variant and reference counts, respectively. If we presume that the rate of mosaicism (the fraction of altered cells, assumed as roughly 2%), ρ, is much larger than the Phred-determined likelihood of base-calling error for any single read, ε, i.e. ρ/ε>>1, then the conditional probabilities are binomial, P(D|S)=C^T,V p^V(1−p)^T−V, where $C_{T,V}$ is the number of combinations of T objects chosen V at a time and p=ε for S=N and p=ρ for S=M. Given a prior estimate of P(N)=0.999, we can formulate P(N|D) directly from Bayes’ Rule, from which additional algebraic manipulation and suitable asymptotic approximation show

Candidates are culled as non-mosaic if this probability exceeds 95%, though in practice most cases actually removed have P(N|D)>99.5%. The remaining set of events is subsequently passed to the second step, which is a traditional Fisher exact (table) test for association between sample type and variant allele fraction followed by standard Benjamini-Hochberg FDR assessment. We report events having ≤20% FDR with at least 3 supporting reads and greater than 2% variant allele fraction.

Analysis of NHLBI ESP data

NHLBI variants calls for 6,503 samples were downloaded from the NHLBI Exome Variant Server http://evs.gs.washington.edu/EVS/. All variants were processed using the same tools as for the TCGA cohort. For comparative analysis, all ESP variants are filtered for < 0.1% total MAF to minimize false positives. The Women’s Health Initiative Exome Sequencing Project (WHISP) is part of the National Heart, Lung, and Blood Institute’s (NHLBI), Grand Opportunity Exome Sequencing Project (https://www.fhcrc.org/en/labs/phs/projects/cancer-prevention/projects/whisp.html). WHISP data for 614 samples were downloaded from dbGaP, verified for file integrity, and then imported as BAM files into our data warehouse. Alignment to the reference genome GRCh37-lite was carried out using BWA v0.5.9 with parameters −t 4 −q 5 and marking of duplicates by Picard v1.46. Variant calling was performed as described in the "Variant calling and annotation strategies". For quality control purposes, we included WHISP samples with read mapping rates greater than 80%, duplication rates less than 40%, and at least 10,000 SNVs detected in the target region. The 557 Caucasian WHISP samples selected for this study, on average, had mapping rates of ~95%, duplication rates of ~9%, and ~18,000 SNVs called in the target region.

Deep sequencing validation

Five candidate low-level variants (two R882C and one R882H sample sites in DNMT3A, one V617F in JAK2, and one K700E in SF3B1), 5 positive controls, and 4 negative controls were selected for validation using deep sequencing (Supplementary Table 5c). Primer pairs tailed with sample-specific indexes were designed for individual target sites and further used for PCR amplifications. Indexed libraries were made for tumor and blood pools respectively. We then generated sequencing data using 1 lane of MiSeq run with read length of 2×250. Custom references were created by including specific primer and index sequences. MiSeq reads were aligned to the custom references using BWA (bwa aln −t 8; bwa sampe). Allelic counts for the variants were obtained using in house tool bam-readcount (bam-readcount −b 0 −q 30).