|Home | About | Journals | Submit | Contact Us | Français|
Acute myeloid leukemia is a highly malignant hematopoietic tumor that affects about 13,000 adults yearly in the United States. The treatment of this disease has changed little in the past two decades, since most of the genetic events that initiate the disease remain undiscovered. Whole genome sequencing is now possible at a reasonable cost and timeframe to utilize this approach for unbiased discovery of tumor-specific somatic mutations that alter the protein-coding genes. Here we show the results obtained by sequencing a typical acute myeloid leukemia genome and its matched normal counterpart, obtained from the patient’s skin. We discovered 10 genes with acquired mutations; two were previously described mutations thought to contribute to tumor progression, and 8 were novel mutations present in virtually all tumor cells at presentation and relapse, whose function is not yet known. Our study establishes whole genome sequencing as an unbiased method for discovering initiating mutations in cancer genomes, and for identifying novel genes that may respond to targeted therapies.
We used massively parallel sequencing technology to sequence the genomic DNA of tumor and normal skin cells obtained from a patient with a typical presentation of FAB M1 Acute Myeloid Leukemia (AML) with normal cytogenetics. 32.7-fold ‘haploid’ coverage (98 billion bases) was obtained for the tumor genome, and 13.9-fold coverage (41.8 billion bases) was obtained for the normal sample. Of 2,647,695 well-supported Single Nucleotide Variants (SNVs) found in the tumor genome, 2,588,486 (97.7%) also were detected in the patient’s skin genome, limiting the number of variants that required further study. For the purposes of this initial study, we restricted our downstream analysis to the coding sequences of annotated genes: we found only eight heterozygous, non-synonymous somatic SNVs in the entire genome. All were novel, including mutations in protocadherin/cadherin family members (CDH24 and PCLKC), G-protein coupled receptors (GPR123 and EBI2), a protein phosphatase (PTPRT), a potential guanine nucleotide exchange factor (KNDC1), a peptide/drug transporter (SLC15A1), and a glutamate receptor gene (GRINL1B). We also detected previously described, recurrent somatic insertions in the FLT3 and NPM1 genes. Based on deep readcount data, we determined that all of these mutations (except FLT3) were present in nearly all tumor cells at presentation, and again at relapse 11 months later, suggesting that the patient had a single dominant clone containing all of the mutations. These results demonstrate the power of whole genome sequencing to discover novel cancer-associated mutations.
Acute Myeloid Leukemia refers to a group of clonal hematopoietic malignancies that predominantly affect middle-aged and elderly adults. An estimated 13,000 people will develop AML in the United States in 2008, and 8,800 will die from it1. Although the life expectancy from this disease has increased slowly over the past decade, the improvement is due predominantly to improvements in supportive care, not in the drugs or approaches used to treat patients.
For most patients with a ‘sporadic’ presentation of AML, it is not yet clear whether inherited susceptibility alleles play a role in pathogenesis2. Further, the nature of the initiating or progression mutations is for the most part unknown3. Recent attempts to identify additional progression mutations by extensively re-sequencing tyrosine kinase genes yielded very few new mutations, and most were not recurrent4,5. Expression profiling studies have yielded signatures that correlate with specific cytogenetic subtypes of AML, but have not yet suggested new initiating mutations6–8. Recent studies using array-based comparative genomic hybridization and/or SNP arrays, while identifying important gene mutations in acute lymphoblastic leukemia9,10 have revealed very few recurrent submicroscopic somatic copy number variants in AML (Walter, MJ, manuscript in preparation, and references 11–13). Together, these studies suggest that we have not yet discovered most of the relevant mutations that contribute to the pathogenesis of AML. We therefore believe that unbiased whole genome sequencing will be required to identify most of these mutations. Until recently, this approach has not been feasible, due to the high cost of conventional capillary-based approaches, and the large numbers of primary tumor cells required to yield the necessary genomic DNA. “Next-generation” sequencing approaches, however, have changed this landscape.
Our group has pioneered the use of whole genome re-sequencing and variant discovery approaches using the Illumina/Solexa technology with the genome of the nematode worm C. elegans as a proof-of-principle14. This approach has distinct advantages in reduced cost, dramatically increased data production rate, and a low input requirement of DNA for library construction. In the present study, we used a similar approach to sequence the tumor genome of a single AML patient, and the matched normal genome (derived from a skin biopsy) of the same patient. Following alignment to the human reference genome, sequence variants were discovered in the tumor genome and compared to the patient’s normal sequence, to dbSNP, and to variants recently reported for two additional human genomes15,16; revealing novel single nucleotide and small insertion/deletion (“indel”) variants genome-wide. Novel somatic mutations were detected in genes not previously implicated in AML pathogenesis, demonstrating the need for unbiased whole genome approaches to discover all mutations associated with cancer pathogenesis.
Of the 8 FAB subtypes of AML, M1 AML is one of the most common (~20% of all cases). No specific cytogenetic abnormalities or somatic initiating mutations have been identified for this subtype; in fact, about half of the patients with de novo M1 AML have normal cytogenetics17–19. The frequency of well-described progression mutations (e.g. activating alleles of FLT3, c-Kit, and Ras) is similar to that of other common FAB subtypes5. We therefore decided to sequence the genome of tumor cells derived from a patient with M1 AML, since so little is known about the molecular pathogenesis of this common subtype. The criteria used to select the sample are outlined in Supplementary Materials.
The case presentation is described in detail in the Supplementary Materials. Briefly, a previously healthy woman in her mid-50s presented suddenly with fatigue and easily bruisability, and was found to have a peripheral WBC count of 105,000 cells per microliter, with 85% myeloblasts. A bone marrow examination revealed 100% myeloblasts with morphologic features and cell surface markers consistent with FAB M1 AML (Supplementary Figure 1). Cytogenetic analysis of tumor cells revealed a normal 46 XX karyotype. Although the patient experienced a complete remission with conventional therapies, she relapsed at 11 months and expired 24 months after her initial diagnosis was made. Informed consent for whole genome sequencing was subsequently obtained from her next of kin.
The tumor sample from patient 933124 contained no somatic copy number changes at a resolution of ~5 kb (further confirmed on the NimbleGen 2.1M array platform, data not shown), and no evidence of copy number neutral LOH, indicating that the genome was essentially diploid at this level of resolution (see Supplementary Figure 2). Further analysis of the 933124-derived tumor and skin samples revealed 26 inherited copy number variants (i.e. detected in both the tumor and skin samples). All but two of these previously had been reported in the Database of Genomic Variants (see Supplementary Table 1). All of the CNVs detected in this genome were found in at least one other AML patient (89 additional cases, mostly Caucasian, have been queried using the same SNP array platform), and all but one were found in at least one of the 160 Caucasian HapMap and Coriell samples that were studied on the same array platform (Supplementary Table 1).
To determine whether the tumor cells of 933124 were typical of M1 AML, we compared the expression signatures of 111 de novo AML cases using unsupervised clustering (Ward’s method, see Materials and Methods). 933124 clustered with multiple other M1 (and M2) AML cases with normal cytogenetics, suggesting that the genetic events underlying the pathogenesis of this case are similar to those of other cases exhibiting normal cytogenetics (Supplementary Fig. 3).
Since most of the acquired mutations in cancer genomes have been shown to be heterozygous, the complete sequencing of a cancer genome requires the detection of both alleles at most positions in the genome20. We therefore designed sequence coverage metrics to define the point at which 90% diploid coverage had been reached. To minimize errors associated with any single platform or measurement, diploid coverage for this genome was assessed using a set of High-Quality (HQ) SNPs derived from two different SNP array platforms, Affymetrix 6.0 and Illumina Infinium 550K. For a SNP to be included in the HQ set, the following criteria had to be satisfied: (1) identical genotypes were called from both assays at the same genomic positions and (2) the resulting genotype was heterozygous. For the 933124 tumor genome, 46,494 heterozygous SNPs passed the above criteria and were defined as HQ SNPs. 46,572 HQ SNPs were defined for the skin sample.
We performed 98 full runs on the Illumina Genome Analyzer to achieve the targeted level of 90% diploid coverage as determined by coverage of the HQ SNP set. Maq21 was used to perform alignment, determine consensus, and identify SNVs within the 98 billion bases generated from the tumor genome (see Table 1). Maq predicted a total of 3.81M SNVs (Maq SNP Quality ≥ 15) in the tumor genome, including matching heterozygous genotypes for 91.2% of the 46,494 HQ SNPs. When we lowered the Maq SNP Quality cutoff to 0, 94.06% HQ SNPs were predicted. Further investigation of Maq alignments revealed coverage for both alleles at an additional 5.38% of the HQ SNPs, but Maq did not predict a SNP or matching heterozygous genotype due to insufficient depth or quality of coverage. Additional analysis revealed coverage at 46,484 of 46,494 HQ SNPs for at least one allele (i.e., 99.98% haploid coverage for the tumor genome).
We sequenced the genome of normal skin cells from the same patient to enable the identification of inherited sequence variants in the tumor genome. Our targeted diploid coverage goal for the skin-derived genome was 80%. We achieved this goal with only 34 Solexa runs (41.8 billion bases), utilizing improved reagents and longer read lengths to attain 82.6% diploid and 84.2% haploid coverage (Table 1).
To begin evaluating the quantity and quality of the detected sequence variants in the tumor and skin genomes, we compared the overlap and uniqueness of this genome’s variants with respect to the Watson and Venter genomes, and to dbSNP (v127) (Figure 1). Of the 3.68 M single nucleotide variants (SNVs, Maq SNP Quality ≥ 15, excluding SNVs found on chromosome X) predicted by Maq in the tumor genome, 2.36 M were present in dbSNP, 2.36 M were detected in the skin genome (Fig. 1A), 1.50 M were detected in the Venter genome, and 1.58 M were found in the Watson genome (Fig. 1B). Ultimately, 1.70 M SNVs were unique to the 933124 tumor genome. Upon filtering the 933124 SNVs at different Maq quality values to determine the stability of results, we observed that the proportion of 933124 SNVs that also are in dbSNP increases from 63.9% to 69.48% when the Maq quality threshold score increases from 15 to 30, as expected.
Because the number of sequence variants initially detected by Maq was high, we developed improved filtering tools to effectively separate true variants from false positives. To this end, we generated an experimental dataset by re-sequencing Maq-predicted SNVs, randomly selecting a training subset and a test dataset, whose annotations and features were submitted to Decision Tree C4.522. This approach identified parameters that separated true variants from false positives, revealing that SNV-supporting read counts (unique based on read start position, and base position in supporting reads), base quality, and Maq quality scores are major determinants for identifying false positives. Implementing rules obtained from the Decision Tree analysis resulted in 91.9% sensitivity and 83.5% specificity for validated SNVs.
The patient had 3,813,205 sequence variants in her tumor genome, as defined by Maq scores of >15 (Table 1). Of these, 2,647,695 were supported by the Decision Tree analysis in the tumor genome, of which 2,584,418 (97.7%) also were detected in the skin genome (Figure 2). The detailed algorithm for selecting putative somatic variants is described in Supplemental Materials. Most of the 63,277 tumor-specific variants we detected were either present in dbSNP or were previously described in the Watson or Venter genomes (31,645), or occurred in non-genic regions (20,440). A total of 11,192 variants were located within the boundaries of annotated genes; 216 of these variants were in untranslated regions, and 10,735 were in introns (but not involving splice junctions) and were not explored further in our analysis. Of the coding sequence variants, 60 were synonymous, and not further evaluated. The remaining 181 variants were either non-synonymous, or were predicted to alter splice site function. By sequencing PCR-generated amplicons from the tumor and skin samples (and also from the relapse tumor sample obtained 11 months after the original presentation), we determined that 152 of these variants were false positive (i.e. wild type) calls, 14 were inherited SNPs, and eight were somatic mutations in both the original tumor and the relapse sample (Table 2). Seven variants could not be validated, either because the regions involved were repetitive, or because all attempts to obtain PCR amplicons failed. All of the PCR-amplified exons from the eight genes containing validated somatic mutations were sequenced in 187 additional cases of AML using samples from our discovery and validation sets23; no additional somatic mutations were detected in these genes (data not shown). A description of how we estimated the false negative (12.45%) and false positive (0.06%) rates for SNVs over the entire genome is presented in Supplementary Materials. Using these estimates, we can predict that very few somatic, non-synonymous variants were missed by our analysis of this deeply covered genome.
To better define the percentage of tumor cells that contained each of the discovered somatic mutations, we amplified each mutation-containing locus from non-amplified genomic DNA derived from the de novo and relapse tumor samples, and from the skin biopsy obtained at presentation. The resulting amplicons were sequenced using the Roche/454 FLX platform, and the frequency of reads containing the reference and variant alleles were defined (Figure 3, and Table 3). Control amplicons containing a known heterozygous SNP in BRCA2 (encoding N372H) and a homozygous SNP in TP53 (encoding P72R) were analyzed similarly. The BRCA2 SNP yielded ~50% variant frequencies in the tumor and skin samples, while nearly 100% of the TP53 alleles were variant in all three samples, as expected. Remarkably, all eight somatic SNVs were detected at ~50% frequencies in the primary tumor sample (100% blasts), and at ~40% frequencies in the relapse sample (78% blasts; if the variant frequencies are corrected for blast counts [i.e. multiplied by 1.28], the frequencies at relapse also were ~50%). The NPMc mutation also was detected at a frequency of ~50%, but the FLT3 ITD allele was detected in only 35.1% of the 454 reads at diagnosis, and 31.3% at relapse, suggesting that the mutation was not present in all tumor cells at diagnosis or relapse.
Surprisingly, the variant alleles also were detected at frequencies of ~5–10% in the skin sample (the PCLKC gene was the single outlier, with 23% variant reads in the skin sample for unclear reasons). In retrospect, it is clear that the skin sample contained contaminating leukemic cells, since the patient’s WBC count at presentation was 105,000 per microliter, with 85% blasts. This information was used to inform the Decision Tree analysis described above: we allowed high quality tumor variants to move forward in the discovery pipeline if they were detected at a low frequency (two or fewer reads) in the skin sample (as defined by a binomial test).
To discover small indels (<6 bp) from sequence reads (32–35 bp long), we started with a set of 236 million reads that were not confidently aligned by Maq to the reference genome. We applied Cross_Match and BLAT to identify gapped alignments that are unique in the genome. To detect indels longer than 6 bp, we developed a novel “split reads” algorithm (see Supplementary Materials) that aligns sub-segments of reads independently to the genome, and computes a mapping quality for the derived gapped alignment based on the number of hits and the quality of the bases. These efforts resulted in the identification of 726 putative small indels (1 to 30 bp in size) that occur in coding exons, 393 of which (54.2%) were found in dbSNP. After manual review, we selected a set of 28 putative somatic coding indels for validation using PCR-based dye terminator sequencing. Of these putative indels, 22 were validated, but were found present in both tumor and skin (15 of these were in dbSNP), 2 were false positive calls, 2 had no coverage, and two were previously validated somatic insertions in NPM1 (4 bp) and FLT3 (30 bp).
In this report, we describe the sequencing and analysis of a primary human cancer genome using next-generation sequencing technology. Our patient’s tumor genome was essentially diploid, and contained ten non-synonymous somatic mutations that may be relevant for her disease. These mutations affect genes participating in several well-described pathways that are known to contribute to cancer pathogenesis, but most of these genes would not have been candidates for directed re-sequencing based on our current understanding of cancer. Hence, these results justify the use of next-generation whole genome sequencing approaches to reveal somatic mutations in cancer genomes.
As we demonstrated in our re-sequencing of the genome of the C. elegans N2 Bristol strain14, and again in this study, massively-parallel short-read sequencing provides an effective method for examining single nucleotide and short indel variants by comparison of the aligned reads to a reference genome sequence. By sequencing our patient’s tumor genome to a depth of >30-fold coverage, and gauging our ability to detect known heterozygous positions across the genome, we have produced a sufficient depth and breadth of sequence coverage to comprehensively discover somatic genome variants. A slightly lower coverage of the normal genome from this individual helped to identify nearly 98% of potential variants as being inherited, a critical filter that allowed us to more readily identify the true somatic mutations in this tumor. Our results strongly support the notion that hypothesis-driven (e.g. candidate gene-based) examination of tumor genomes by PCR-directed or capture-based methods is inherently limited, and will miss key mutations. An additional and important consideration is the demand for large amounts of genomic DNA by these techniques; this is a serious limitation when precious clinical samples are being studied. The Illumina/Solexa technology requires only ~1 ug of DNA per library, enabling the study of primary tumor DNA rather than requiring the use of tumor cell lines, which may contain genetic changes and adaptations required for immortalization and maintenance in tissue culture conditions.
A total of 10 non-synonymous somatic mutations were identified in this patient’s tumor genome. Two are well known AML-associated mutations, including an internal tandem duplication of the FLT3 receptor tyrosine kinase gene, which constitutively activates kinase signaling, and portends a poor prognosis5,24,25, and a four base insertion in exon 12 of the NPM1 gene (NPMc)26–28. Both of these mutations are common (25–30%) in AML tumors, and both are thought to contribute to progression of the disease, rather than to cause it directly29. Interestingly, the frequency of the mutant FLT3 allele in the primary and relapse tumor samples (35.08% and 31.30%, respectively) was significantly less than that of the other 9 mutations (p<0.000001 for both the primary and relapse samples). These data suggest that the FLT3 ITD may not have been present in all tumor cells, and further, that it may have been the last mutation acquired.
The other eight somatic mutations that we detected are all single base changes, and none has previously been detected in an AML genome. Four of the genes affected, however, are in gene families that are strongly associated with cancer pathogenesis (including PTPRT, CDH24, PCLKC, and SLC15A1). The other four somatic mutations occurred in genes not previously implicated in cancer pathogenesis, but whose potential functions in metabolic pathways suggest mechanisms by which they could act to promote cancer (including KNDC1, GPR123, EBI2, and GRINL1B). We speculate regarding the roles of these mutations for the pathogenesis of this patient’s disease in Supplementary Materials.
The importance of the eight newly defined somatic mutations for AML pathogenesis is not yet known, and will require functional validation studies in tissue culture cells and mouse models to assess their relevance. Even though we could not detect recurrent mutations in the limited AML sample set we surveyed, several lines of evidence suggest that these mutations may not be random, “passenger” mutations, as follows:
Although we performed whole genome sequencing on this cancer sample, we restricted our initial validation studies to the 1–2% of the genome that encodes genes. This raises the issue of whether sequencing the cDNA transcriptome of this tumor would have been a faster, cheaper, and more efficient way of finding the mutations. While this approach will undoubtedly be an important adjunct to whole genome sequencing, there are several advantages to the approach we used:
The additional non-coding and non-genic somatic variants in this genome (which we currently estimate at 500–1000, based on our assessment of false positive and negative rates for non-synonymous mutations), which will be fully described later, will provide a rich source of potentially relevant sequence changes that will be better understood as more cancer genomes are sequenced.
In summary, we have successfully used a next-generation whole genome sequencing approach to identify new candidate genes that may be relevant for AML pathogenesis. We cannot overemphasize the importance of parallel sequencing of the patient’s normal genome to determine which variants were inherited; the identification of the true somatic mutations in this tumor genome would not have been feasible without this approach. Furthermore, until hundreds (or perhaps thousands) of normal genomes and additional AML tumors are sequenced, the contextual relevance of the mutations found in this genome will be unknown. Regardless, the somatic mutations that we did find were neither predicted by the curation of previously defined cancer genes, nor by the study of this tumor using unbiased, high-resolution array-based genomic approaches. For AML and other types of cancer, whole genome sequencing may therefore be the only effective means for discovering all of the mutations that are relevant for pathogenesis.
Sequence end reads (average length for tumor genome, 32 bp, and for skin, 35 bp) were generated from Illumina/Solexa fragment libraries derived from the tumor or skin cells of patient 933124, using the Illumina Genome Analyzer. The analyzed reads were aligned to the human reference (NCBI Build 36) using Maq21. Coverage of the tumor and normal genomes was ascertained by comparison to the patient’s heterozygous SNPs, established by compiling shared SNP calls monitored on the Affymetrix 6.0 and Illumina Infinium 550K genotyping platforms. We examined the Maq alignments by Decision Tree analysis to discover single nucleotide variants, as well as to identify copy number variants. Nonaligned reads were further analyzed for indel discovery. For all putative variants, we attempted validation using custom PCR and capillary sequencing on the ABI 3730 platform. All validated somatic mutations were further analyzed by Roche/454 sequencing of PCR-generated amplicons made from primary genomic DNA to compare readcounts of wild-type and mutant alleles in the primary tumor, skin, and relapse tumor samples. A complete description of the AML case sequenced, and the Materials and Methods used to generate this dataset, are provided in the Supplementary Materials.
The authors dedicate this work to our AML patients and their families, and to Alvin J. Siteman, whose generous and visionary gift provided the major funding source for this study. We also thank Drs. Gerry Flance, David Kipnis, and Kenneth Polonsky for their support. The authors also thank Drs. Clara Bloomfield, Michael Caligiuri, and James Vardiman from the Cancer and Leukemia Group B for providing important AML samples for validation studies. We also thank the staff of The Genome Center at Washington University for their support of and their many contributions to this project. Additional funding was provided by the National Cancer Institute (TJL), the National Human Genome Research Institute (RKW), and the Barnes-Jewish Hospital Foundation (TJL).
Sequence Variant deposition in dbGaP
High quality sequence variants defined by Decision Tree (2,647,695 variants) will be deposited in the dbGaP database http://www.ncbi.nlm.nih.gov/sites/entrez?Db=gap) for review by qualified investigators.