Main

We used massively parallel sequencing technology to sequence the genomic DNA of tumour and normal skin cells obtained from a patient with a typical presentation of French–American–British (FAB) subtype M1 acute myeloid leukaemia (AML) with normal cytogenetics. For the tumour genome, 32.7-fold ‘haploid’ coverage (98 billion bases) was obtained, and 13.9-fold coverage (41.8 billion bases) was obtained for the normal skin sample. Of the 2,647,695 well-supported single nucleotide variants (SNVs) found in the tumour genome, 2,584,418 (97.6%) were also 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 new, including mutations in protocadherin/cadherin family members (CDH24 and PCLKC (also known as PCDH24)), G-protein-coupled receptors (GPR123 and EBI2 (also known as GPR183)), 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. On the basis of deep readcount data, we determined that all of these mutations (except FLT3) were present in nearly all tumour 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 new cancer-associated mutations.

AML refers to a group of clonal haematopoietic 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 predominantly because of 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 have a role in the pathogenesis2. Furthermore, 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 previously unidentified 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,7,8. Recent studies using array-based comparative genomic hybridization and/or single nucleotide polymorphism (SNP) arrays, although identifying important gene mutations in acute lymphoblastic leukaemia9,10 have revealed very few recurrent submicroscopic somatic copy number variants in AML (M.J.W., manuscript in preparation, and refs 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 because of the high cost of conventional capillary-based approaches and the large numbers of primary tumour 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 Caenorhabditis elegans as a proof-of-principle14. This approach has distinct advantages in reduced cost, a markedly 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 tumour genome of a single AML patient and the matched normal genome (derived from a skin biopsy) of the same patient. After alignment to the human reference genome, sequence variants were discovered in the tumour genome and compared to the patient’s normal sequence, to the dbSNP database, and to variants recently reported for two other human genomes15,16; revealing new single nucleotide and small insertion/deletion (indel) variants genome-wide. 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.

Rationale for using the FAB M1 AML subtype for sequencing

Of the eight 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,18,19. The frequency of well-described progression mutations (for example, activating alleles of FLT3, KIT and RAS) is similar to that of other common FAB subtypes5. We therefore decided to sequence the genome of tumour cells derived from a patient with M1 AML, because so little is known about the molecular pathogenesis of this common subtype. The criteria used to select the sample are outlined in Supplementary Information.

Case presentation of UPN 933124

The case presentation is described in detail in the Supplementary Information. In brief, a previously healthy woman in her mid-50s presented suddenly with fatigue and easy bruisability, and was found to have a peripheral white blood cell count of 105,000 cells per microlitre, with 85% myeloblasts. A bone marrow examination revealed 100% myeloblasts with morphological features and cell surface markers consistent with FAB M1 AML (Supplementary Fig. 1). Cytogenetic analysis of tumour 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. At relapse, the bone marrow had 78% myeloblasts, and contained a new clonal cytogenetic abnormality, t(10; 12) (p12; p13). Informed consent for whole-genome sequencing was subsequently obtained from her next of kin.

A typical M1 AML diploid genome and expression profile

The tumour 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 loss-of-heterozygosity (LOH), indicating that the genome was essentially diploid at this level of resolution (see Supplementary Fig. 2). Further analysis of the 933124-derived tumour and skin samples showed 26 inherited copy number variants (that is, detected in both the tumour and skin samples). All but two of these had been previously reported in the Database of Genomic Variants (see Supplementary Table 1). All of the copy number variants detected in this genome were found in at least one other AML patient (89 other 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 tumour 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 Supplementary Information). The expression profile of patient 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).

Coverage depth of the tumour and skin genomes

Because 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 SNPs derived from two different SNP array platforms, Affymetrix 6.0 and Illumina Infinium 550K. For a SNP to be included in the high-quality 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 tumour genome, 46,494 heterozygous SNPs passed the above criteria and were defined as high-quality SNPs. For the skin samples, 46,572 high-quality SNPs were defined.

We performed 98 full runs on the Illumina Genome Analyser to achieve the targeted level of 90% diploid coverage as determined by coverage of the high-quality SNP set. Maq21 was used to perform alignment, determine consensus, and identify SNVs within the 98 billion bases generated from the tumour genome (see Table 1). Maq predicted a total of 3.81 million SNVs (Maq SNP quality ≥ 15) in the tumour genome, including matching heterozygous genotypes for 91.2% of the 46,494 high-quality SNPs. When we lowered the Maq SNP quality cutoff to 0, 94.06% high-quality SNPs were predicted. Further investigation of Maq alignments revealed coverage for both alleles at a further 5.38% of the high-quality SNPs, but Maq did not predict a SNP or matching heterozygous genotype owing to insufficient depth or quality of coverage. Extra analysis revealed coverage at 46,484 of 46,494 high-quality SNPs for at least one allele (that is, 99.98% haploid coverage for the tumour genome).

Table 1 Tumour and skin genome coverage from patient 933124
Full size table

We sequenced the genome of normal skin cells from the same patient to enable the identification of inherited sequence variants in the tumour 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), using 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 tumour and skin genomes, we compared the overlap and uniqueness of this genome’s variants with respect to the James D. Watson and J. Craig Venter genomes, and to dbSNP (v127; Fig. 1). Of the 3.68 million single nucleotide variants (SNVs; Maq SNP quality ≥15, excluding SNVs found on chromosome X) predicted by Maq in the tumour genome, 2.36 million were present in dbSNP, 2.36 million were detected in the skin genome (Fig. 1a), 1.50 million were detected in the Venter genome, and 1.58 million were found in the Watson genome (Fig. 1b). Ultimately, 1.70 million SNVs were unique to the 933124 tumour genome. On 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.

Figure 1: Overlap of SNPs detected in 933124 and other genomes.
figure 1

a, Venn diagram of the overlap between SNPs detected in the 933124 tumour genome and the genomes of J. D. Watson and J. C. Venter. b, Venn Diagram of the overlap among the 933124 tumour genome, the skin genome and dbSNP (ver. 127). SNVs were defined with a Maq SNP quality ≥15.

PowerPoint slide

Full size image

Refining the detection of potential somatic mutations

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 data set by re-sequencing Maq-predicted SNVs, randomly selecting a training subset and a test data set, whose annotations and features were submitted to Decision Tree C4.5 (ref. 22). This approach identified parameters that separated true variants from false positives, revealing that SNV-supporting read counts (unique on the basis of read start position and base position in supporting reads), base quality and Maq quality scores are chief 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.

Identification of somatic mutations in coding sequences

The patient had 3,813,205 sequence variants in her tumour genome, as defined by Maq scores of >15 (Table 1). Of these, 2,647,695 were supported by the Decision Tree analysis in the tumour genome, of which 2,584,418 (97.6%) were also detected in the skin genome (Fig. 2). The detailed algorithm for selecting putative somatic variants is described in Supplementary Information. Most of the 63,277 tumour-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 polymerase chain reaction (PCR)-generated amplicons from the tumour and skin samples (and also from the relapse tumour sample obtained 11 months after the original presentation), we determined that 152 of these variants were false positive (that is, wild type) calls, 14 were inherited SNPs, and eight were somatic mutations in both the original tumour 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 further cases of AML using samples from our discovery and validation sets23; no further 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 Information. Using these estimates, we can predict that very few somatic, non-synonymous variants were missed by our analysis of this deeply covered genome.

Figure 2: Filters used to identify somatic point mutations in the tumour genome.
figure 2

See text for details. UTR, untranslated regions.

PowerPoint slide

Full size image
Table 2 Non-synonymous somatic mutations detected in the AML sample
Full size table

Defining mutation frequencies in the tumour sample

To better define the percentage of tumour 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 tumour 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 (Fig. 3 and Table 3). Control amplicons containing a known heterozygous SNP in BRCA2 (encoding N372H) and a homozygous SNP in TP53 (encoding P72R) were analysed similarly. The BRCA2 SNP yielded ∼50% variant frequencies in the tumour and skin samples, whereas 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 tumour sample (100% blasts), and at ∼40% frequencies in the relapse sample (78% blasts; if the variant frequencies are corrected for blast counts—that is, multiplied by 1.28—the frequencies at relapse also were ∼50%). The NPMc (cytoplasmic nucleophosmin) mutation was also detected at a frequency of ∼50%, but the FLT3 internal tandem duplication (ITD) allele was only detected in 35.1% of the 454 reads at diagnosis and 31.3% at relapse, suggesting that the mutation was not present in all tumour cells at diagnosis or relapse.

Figure 3: Summary of Roche/454 FLX readcount data obtained for ten somatic mutations and two validated SNPs in the primary tumour, relapse tumour and skin specimens.
figure 3

The readcount data for the variant alleles in the primary tumour sample and relapse tumour sample are statistically different from that of the skin sample for all mutations (P < 0.000001 for all mutations, Fisher’s exact test, denoted by a single asterisk in all cases). Note that the normal skin sample was contaminated with leukaemic cells containing the somatic mutations. The patient’s white blood cell count was 105,000 (85% blasts) when the skin punch biopsy was obtained.

PowerPoint slide

Full size image
Table 3 454 Readcount data for somatic mutations and known SNPs
Full size table

Notably, the variant alleles also were detected at frequencies of ∼5–13% in the skin sample. In retrospect, it is clear that the skin sample contained contaminating leukaemic cells, because the patient’s white blood cell count at presentation was 105,000 per microlitre, with 85% blasts. This information was used to inform the Decision Tree analysis described above: we allowed high-quality tumour 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.

Detecting insertions and deletions (indels)

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 ‘split reads’ algorithm (see Supplementary Information) that aligns sub-segments of reads independently to the genome, and computes a mapping quality for the derived gapped alignment on the basis of 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 tumour and skin (15 of these were in dbSNP), two were false positive calls, two had no coverage, and two were previously validated somatic insertions in NPM1 (4 bp) and FLT3 (30 bp).

Discussion

Here we describe the sequencing and analysis of a primary human cancer genome using next-generation sequencing technology. Our patient’s tumour 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 on the basis of 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 tumour 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 tumour. Our results strongly support the notion that hypothesis-driven (for example, candidate gene-based) examination of tumour genomes by PCR-directed or capture-based methods is inherently limited, and will miss key mutations. A further 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 μg of DNA per library, enabling the study of primary tumour DNA rather than requiring the use of tumour cell lines, which may contain genetic changes and adaptations required for immortalization and maintenance in tissue culture conditions.

A total of ten non-synonymous somatic mutations were identified in this patient’s tumour genome. Two are well-known AML-associated mutations, including an internal tandem duplication of the FLT3 receptor tyrosine kinase gene, which constitutively activates kinase signalling, and portends a poor prognosis5,24,25, and a four-base insertion in exon 12 of the NPM1 gene (NPMc)26,27,28. Both of these mutations are common (25–30%) in AML tumours, and are thought to contribute to progression of the disease rather than to cause it directly29. Notably, the frequency of the mutant FLT3 allele in the primary and relapse tumour samples (35.08% and 31.30%, respectively) was significantly less than that of the other nine 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 tumour 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 about the roles of these mutations for the pathogenesis of this patient’s disease in Supplementary Information.

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 that we surveyed, several lines of evidence suggest that these mutations may not be random, ‘passenger’ mutations. First, somatic mutations in this genome are extremely rare. The rarity of somatic variants, and the normal diploid structure of the tumour genome, argues strongly against genetic instability or DNA repair defects in this tumour. Conceptually, this result is further supported by the very small number of somatic mutations discovered in the expressed tyrosine kinases of AML samples4,5; genetic instability does not seem to be a general feature of AML genomes.

Second, on the basis of the equivalent frequencies of the variant and wild-type alleles for the mutations in the tumour genome (except for FLT3 ITD), it is highly probable that all the mutations are heterozygous, and are present in virtually all of the tumour cells (Fig. 3). The latter suggests that these mutations may have all been selected for and retained because they are important for disease pathogenesis in this patient. Alternatively, all may have occurred simultaneously in the same leukaemia-initiating cell, but only a subset of the mutations (or an as-yet undetected mutation) is truly important for pathogenesis (that is, disease ‘drivers’ versus passengers). Although we suggest that the latter hypothesis is very unlikely on the basis of our current understanding of tumour progression, many more AML genomes will need to be sequenced to resolve this issue.

Third, the same mutations were detected in tumour cells in the relapse sample at approximately the same frequencies as in the primary sample. All of these mutations were therefore present in the resistant tumour cells that contributed to the patient’s relapse, further suggesting that a single clone contains all ten mutations. Fourth, seven of the ten genes containing somatic mutations were detectably expressed in the tumour sample. FLT3 and NPM1 messenger RNAs were highly expressed in this tumour sample, as they are in virtually all AML samples. We detected mRNA from the CDH24, SLC15A1 and EBI2 genes on the Affymetrix expression array, whereas expression of GRINL1B and PCLKC were detected by PCR with reverse transcription (RT–PCR; data not shown). Expression of KNDC1, PTPRT and GPR123 was not detected by either approach, but we cannot rule out expression of these genes in a small subset of tumour cells (for example, leukaemia-initiating cells). Furthermore, for the five point mutations where data are available, the mutated base is highly conserved across multiple species (Table 2).

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 complementary DNA transcriptome of this tumour would have been a faster, cheaper and more efficient way of finding the mutations. Although this approach will undoubtedly be an important adjunct to whole-genome sequencing, there are several advantages to the approach we used: (1) coverage models for whole-genome libraries are at present better understood than for cDNA libraries, where transcript abundance can vary over many orders of magnitude; (2) even if the transcriptome had been sequenced, extensive characterization of the normal genome would have been required to distinguish inherited variants from somatic mutations; and (3) relevant non-synonymous mutations could be missed by cDNA sequencing, including mutations that result in RNA instability (splice variants, nonsense mutations), and/or mutations in genes expressed at low levels, or in only a small subset of tumour cells.

The additional non-coding and non-genic somatic variants in this genome (which we presently estimate at 500–1,000 on the basis of our calculated false positive and negative rates for non-synonymous mutations), 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 tumour genome would not have been feasible without this approach. Furthermore, until hundreds (or perhaps thousands) of normal genomes and other AML tumours are sequenced, the contextual relevance of the mutations found in this genome will be unknown. Nevertheless, the somatic mutations that we did find were neither predicted by the curation of previously defined cancer genes, nor by the study of this tumour 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.

Methods Summary

Sequence end reads (average length for tumour genome, 32 bp, and for skin, 35 bp) were generated from Illumina/Solexa fragment libraries derived from the tumour or skin cells of patient 933124, using the Illumina Genome Analyser. The analysed reads were aligned to the human reference genome (NCBI Build 36) using Maq21. Coverage of the tumour 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 SNVs, as well as to identify copy number variants. Non-aligned reads were further analysed 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 analysed 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 tumour, skin and relapse tumour samples. A complete description of the AML case sequenced, and the materials and methods used to generate this data set are provided in the Supplementary Information.

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 approved investigators.