Computational methods and resources for the interpretation of genomic variants in cancer

Cancer mutation data

The Cancer Genome Atlas (TCGA) consortium, which began in United States in 2006, is a comprehensive and coordinated effort to understand the molecular basis of cancer using several genomic analysis techniques. The data generated by the experiments are made available through the TCGA Data Portal and the Cancer Genomics Hub (CGHub) [34]. After signing the certification data agreement it is possible to access data about 36 cancer types. The files containing the binary version of the short DNA sequence read alignments (referred as BAM files) can be downloaded using the GeneTorrent application available at the CGHub website (https://cghub.ucsc.edu/).

Since the downloading of all the BAM files remains costly in term of time and storage, it is advisable to get all the pre-processed variant files (VCFs) generated by different institutions within the TCGA consortium. Unfortunately, so far the organization of the data in the TCGA repository is not optimized and available data have to be manually selected. In particular, variant data from different platforms, ranging from SNP arrays to SOLID and Illumina sequencers, are reported in different file formats. Out of these files, the most informative ones are the VCF files, which contains the genetic variants for both normal and tumor samples. Currently, VCF file are provided only for a subset of the tumor types (~67%).

To study the occurrence of somatic mutations across different donors, we performed a recurrence analysis calculating the Fraction of Somatic Mutations and the Fraction of Donors corresponding to different subsets of somatic mutations (Section 2 in Additional file 1).

The analysis of whole set of somatic mutations (PanCancer) revealed that a large fraction of them are occurring in a single patient and only ~1.7% are recurring more than once (see Table 1). This percentage decreases to 1.2% if the cancer types are considered separately. The plot in Figure 3A shows the Fraction of Somatic Mutations at different levels of Mutation Recurrence.

In addition, we calculated the Fraction of Donors covered by subsets of recurrent mutations (Figure 3B). Considering only the recurrent mutations, (observed at least in two donors), they are held by 96% of the individuals in all cancer types (Table 1). This percentage decreases to 82% when the cancer types are considered separately. In Figure 3C we plot the Fraction of Somatic Mutations and the Fraction of Donors of affected at different Mutation Recurrence thresholds. This curve allows us to estimate the fraction of donors affected by a subset of recurrent mutations that are more likely to have a functional impact. Although it is well accepted that each cancer sample is different and only a small fraction of variants are recurrent, we show that with mutations recurring more than 30 times (53 variants) explains a comparatively large fraction of patients (~40%). However, this data could be biased toward cancer types with higher number of samples and with higher of mutations detected.

We performed the recurrence analysis on each cancer type fitting the points calculated at different Mutation Recurrence threshold. The results in Figure 4 show different trend in 27 tumor types for which the regression curves can be calculated. Although in some cases the fitting is less accurate (THCA and PBCA), we can compare different tumor types estimating the fraction of somatic mutations covering 95% of the donors (Supplementary Table 2, Additional file 1). The smallest value corresponds to the Esophageal Adenocarcinoma (ESAD) for which we estimated somatic mutation rate of ~0.006%. The highest value is reported for the Pediatric Brain Cancer (PBCA) for which 75% the mutations are needed to cover 95% of the samples.

In general, it is expected that heterogeneous cancer types show a large variety of recurrent somatic mutations, which is inversely proportional to the rate of decay of the complementary cumulative distribution (CCD). Thus, the analysis of recurrent mutations based on the calculation of the CCD can provide an estimation of the mutational heterogeneity of each cancer type (see Figure 4). Furthermore, the CCDs in Figures 3C and 4 provide an estimation of the maximum fraction of patients that can be recovered using a decreasing pool of recurrent variants.

Focusing on the exonic regions for each cancer type, we observed different number of somatic mutations per individual (Figure 5). The statistic shows a difference of ~2.5 orders in magnitude between the number of somatic mutations in Chronic Myeloid Disorders (CMDI) and Skin Cutaneous Melanoma (SKCM) for which their median values per individual in the exonic regions is 1 and ~400 respectively (Table 1). This result is consistent with the difference in the rate of decay of the CCDs observed for the two cancer types in Figure 4.

In addition, we can use the number recurrent mutations in the exonic regions to estimate the tumor type similarities. Using the subset of mutations repeated more than once in the whole dataset (PanCancer) we built a vector, in which each element represents the number of donors affected by a mutated gene. Only genes harboring the aforementioned recurrent mutations are considered. The gene-based vectors of each cancer type are then used to calculate the cosine similarity (Section 3 in Additional file 1). In Figure 6 we report the tumor dendrogram obtained from 33 cancer types. Broadly, based just on the similarity between the vectors of recurrent mutated genes we can cluster the tumors in two major subgroups of cancer. The figure shows on the first subgroup on the top-left side Colon Adenocarcinoma (COAD), and Rectum Adenocarcinoma (READ), in which high fractions of donors are mutated in APC, TP53 and KRAS, in the same cluster. In the second subgroup we observe 4 different liver cancers (LINC, LIRI, LICA and LIAD), with recurrent mutation in CTNBB1, in the same cluster.

These results confirm the hypothesis that recurrent somatic mutations, which represent only 1.7% of the whole PanCancer mutations, still retain meaningful information about the mutation profile of each cancer type that allows us to cluster them according to their similarity. Finally, it is noteworthy that filtering out non-recurrent somatic variants we significantly reduce the number putative driver mutations to a relatively small set of variants (~170,000), which are held by ~96% of the patients (see Table 1).

Although the analysis of somatic mutation data from multiple projects provides interesting insights about the genetic cause of tumor, we need to remember that data can have biases due to the use of different sequencing protocols and variant-calling procedures. These problems can be amplified with bias in sampling, significantly different size of the cohorts.

Curated databases and cancer gene lists

One problem that has been creating a constant roadblock in developing better tools and methods for cancer genome analysis is the lack of reference benchmark datasets of known cancer variants. The datasets discussed in this section are curated known collections of cancer variants that are used frequently as a substitute of a reference dataset.

The most popular collection of the somatic mutations in cancer is the Catalogue Of Somatic Mutation In Cancer (COSMIC). The last release of COSMIC (version 71, Sep 2014) contains ~2.1 million unique somatic variants detected in ~1 million tumor samples. An important manually curated resource made available through the COSMIC web page is the Cancer Gene Census list. Currently it contains 547 genes for which mutations have been causally implicated in cancer. Out of 547 cancer genes, 85% of them harbor only somatic mutations, 7% harbor only germline variants and the remaining 8% harbor both types.

Beside the Cancer Gene Census there are other lists of putative cancer genes identified using computational or expert-based approaches. One is a list of 125 cancer genes, reported in a publication by Vogelstein and collaborators [9]. The list comprises 54 oncogenes (43%) and 71 tumor suppressor genes (57%). A larger list of 2,125 genes collected by Bushman Lab [45] has been obtained as a union of 8 different datasets. More recently, an analysis of 21 cancer types was published. The authors prioritized the cancer causing genes according to their observed mutation frequency across multiple samples [46]. The final list thus obtained, combining all the cohorts, contains 146 genes, 64 of which are highly significantly mutated (44%), 49 significantly mutated (34%) and 32 near significant (22%). A better description of the available lists of putative cancer genes obtained from the analysis large cancer studies is provided in recent publications [47–49].

Other useful databases for the annotation of cancer genome are NCG [50] and Cancer3D [51]. The NCG 4.0 database contains 2,000 cancer genes, 537 of which have been experimentally validated and 1,463 statistically inferred from the analysis of candidate cancer genes in 3460 human exomes and genomes from 23 different cancer types. For each cancer gene, NCG provides information about duplicability, evolutionary origin, expression, functional annotation, interaction network with human proteins and microRNAs. Finally, protein structure is also a valuable piece of information for predicting the impact genetic variants [52] and effectiveness of drugs [53, 54]. Cancer3D database collects ~275,000 somatic mutations mapped to ~15,000 proteins that correspond to more than 24,300 structure from the Protein Data Bank [55]. The idea is to study the impact of missense somatic mutations on protein three-dimensional structure. The Cancer3D database includes predictions from e-Drug and e-Driver, two recently developed tools for predicting sensitivity to 24 anticancer compounds and cancer driver proteins [56].

Resources for the visualization and analysis of cancer genome data

In the last few years, several web portals have been developed for the visualization and analysis of all kinds of data in cancer genome studies. Among the most used is the cBioPortal maintained by the Memorial Sloan Kettering Cancer Center [57, 58]. The cBioPortal can be queried either using the web interface or through programmable application programming interface (API). The cancer data displayed by the cBioPortal are extracted from 18,469 samples from 80 cancer genome studies. To facilitate the submission of queries, cBioPortal implements the Onco Query Language (OQL) for selecting and defining genetic alterations for any subset of data. The user interface allows searching for somatic mutations or copy number alteration across different samples and cancer studies selecting variations occurring in user-defined gene list or in 22 predefined groups cancer-specific genes. When available, expression level from RNA-Seq experiments, methylation and protein phosphorylation data can also be retrieved.

A possible alternative to cBioPortal is the Cancer Genome Browser [59]. The Cancer Genome Browser allows extraction and visualization of information from 574 datasets. Among these two resources, cBioPortal provides better integration with pathways data and the Cancer Genome Browser allows better visualization of clinical data.

Other three interesting resources for the analysis of cancer genome data are the Cancer Gene annotation system for Cancer Genomics (CaGe) [60], the DriverDB [61] and the Integrative Oncogenomics Cancer Browser (IntOGen) [62]. CaGe is a cancer genome annotation system for the classification of candidate genes from cancer genome studies, using either previously reported or novel categories of cancer genes, providing insights about the underlying carcinogenic mechanisms through pathway analysis. DriverDB incorporates data from 6,079 exome-sequencing experiments from 33 cancer studies. It integrates annotation databases and eight bioinformatics algorithms for detecting driver genes and mutations. Finally, Integrative OncoGenomics (IntOGen) is a web available resource integrating data from 6,792 genomes/exome sequencing experiments from 28 cancer types (release 2014.12) [63]. It collects and analyses somatic mutations in a large set of tumor samples to identify putative cancer driver genes. The prediction of putative cancer driver genes is performed by OncodriveFM [64], which has been developed by member of the same group (see next section). The web interface allows browsing data by gene name, cancer type, sites and projects, reporting the frequencies of mutation at gene and genomic location InOGen also allows the study of genomic alterations in cancer in the contest of pathways [65].

The URLs of the cancer genome data repositories, the cancer gene lists and resources for cancer genome analysis are summarized in Table 2.

Detection of cancer driver genes

In general, the insurgence of cancer is accompanied with an accumulation of somatic mutations. However, this does not imply that all mutations are of equal importance in the oncogenesis and cancer progression. Rather, driver somatic mutations exert a selective advantage to cancer cells. In contrast, passenger mutations are abundant but do not confer any selective advantage. Driver variants are present in a small fraction that may vary in different cancer types [66, 67]. Thus, an obvious approach for the selection of cancer driver genes involves the analysis of recurrent somatic mutations. Following this idea, several methods - such as MuSiC [68], MutSigCV [69] DrGaP [70] and Simon's tool [71] - have been developed to prioritize cancer driver genes using different statistical models for the estimation of a background mutation rate.

The Mutational Significance in Cancer (MuSiC) uses sequence-based input to establish correlations among mutation sites, affected genes and pathways, and to discriminate abundant passenger mutations from significant mutation events. This method aims to identify significantly mutated genes with respect to a background mutation rate. MuSiC has been tested on a set of 316 ovarian cancer samples and it was able to detect 12 significantly mutated genes [68]. MutSigCV is the newest version of the "Mutation Significance" algorithm that uses gene specific background mutation rates including mutation events in gene covariates. The information from known co-varying genes is important for estimating the background mutation rates in the genomic regions where few mutation events are observed. MutSigCV takes an input list of mutation from different samples and builds a model gene-dependent background mutation rate estimated using clusters of genes. MutSigCV has been used to analyze 3,083 tumor-normal pairs, and it was able to discover strong differences in mutation frequency and spectrum across 27 cancer types, providing interesting insights about the etiology of the disease [69]. DrGap [70] is a statistical framework for identifying driver genes and signaling pathways in cancer genome-sequencing studies. This approach integrates biological knowledge of the mutational process in tumors and uses a heuristic strategy to optimize the mixture proportion of chi-square distribution of likelihood ratio test (LRT) statistics. This approach improves the accuracy and sensitivity of the prediction of driver genes avoiding zero estimation of the driver mutation rates due to the small probability of observing any mutation in the available samples. Simon's tool [71] calculates the background mutation rate by accounting for the functional impact of mutations on proteins, variation in background mutation rate among tumors, and the redundancy of the genetic code. Using this algorithm, the authors reanalyzed 623 candidate genes in 188 non-small cell lung tumors and identified 28 driver genes, 6 of which were novel [72].

Alternative approaches have been implemented in ActiveDriver [73], OncodriveFM [64], OncodriveCLUST [74] and ContrastRank [75]. ActiveDriver is able to detect significantly mutated functional sites in cancer genome providing an analysis of the somatic mutations associate to phosphorylation signaling. Indeed, the predictive model has been calculated considering the frequency of somatic variants from ~800 cancer genomes in proximity of ~74,000 phosphorylation sites and 469 kinase domains [73]. The method was able to identify candidate genes, protein complexes and kinases involved in cancer. OncodriveFM [64] is a method for the detection of putative cancer driver genes or gene modules. It computes the functional impact of variants using three established computational approaches (SIFT [76], PolyPhen2 [77] and MutationAssessor [78]). OncodriveFM prioritizes putative cancer driver genes calculating the distribution of the predicted functional impact scores for several variants across tumor samples and its deviation from the null model. The bias towards the accumulation of variants with high functional impact is used detect candidate driver genes. Similar approach has been used to prioritize cancer-associated pathways. OncodriveCLUST has been developed to identify significant bias towards somatic mutations clustered within the same protein. The background model for prioritizing genes is calculated using the rate on silent mutations, which are not assumed to be under positive selection. OncodriveCLUST cluster scores have been calculated using ~240,000 mutations in ~9,500 genes from COSMIC database and has been tested of on four TCGA datasets. This approach provides variable background mutation rate for each gene and it detect recessive cancer genes not identified using OncodriveFM [74].

Finally, we also highlight ContrastRank [75], which is a new method for the prioritization of putative impaired genes in cancer. With respect to previously developed methods, ContrastRank evaluates the background mutation rate using the maximum value between the mutation rates of each gene in 1000 Genomes [79] and normal TCGA samples. For the calculation of the score the methods only considers mutation with allele frequency in 1000 Genomes lower than 0.5% and uses the complement of the cumulative binomial distribution to rank cancer associated-genes. The method has been tested on TCGA from three types of adenocarcinomas. In addition, the Cancer Census, Vogelstein and Bushman cancer gene lists (see Table 2) have been used to assess the quality of the prioritization method. ContrastRank was also used for calculating an exome-based score to discriminate between TCGA normal and tumor TCGA samples resulting in high level of accuracy.

Predicting the impact of non-synonymous variants in cancer

During the last decades several methods have been developed to predict the impact of non-synonymous single nucleotide variants (nsSNVs) at structural [80] and functional levels [2]. In particular, the algorithms for predicting the functional effect of missense variants estimate the probability that a mutation is disease-associated or functionally deleterious. Although the relationships between molecular state and disease are complex and are far from being completely understood, the pathologic effect resulting from amino acid substitution is commonly estimated by predicting its functional impact. Most of the available algorithms are based on the evidence that functionally important protein sites are under purifying selection [81]. Therefore, the majority of disease-causing variants should occur in conserved regions that can be detected by evolutionary analysis [82]. Using this approach many classifiers have been implemented to predict whether a nsSNV has any functional impact. The most famous methods for estimating the impact of genetic variants are SIFT [76] and PolyPhen [77] whose predictions are already embedded in many variant annotation pipelines. Recently, more sophisticated methods, which exploit additional structural and functional information, have been developed [19, 83–85]. Recent advances in the field are represented by the implementation of consensus-based algorithms [86, 87] and a general approach that is able to predict the impact of variants in non-coding regions [81, 88].

Although most variant effect predictor reached an adequate level of accuracy, their predictions do not provide information about the possible phenotypic effect. This problem has been partially addressed with the development of a new class of disease-specific predictors trained on a subset of mutations with defined phenotypic effect. In particular, several methods have been developed for discriminating between passenger and cancer driver mutations [78, 89–93]. Among them, CanPredict [91] has been the first algorithm for predicting cancer-causing mutations. The method uses a conservation measures from PFAM [94] Hidden Markov Models and functional information encoded through Gene Ontology terms. Similar strategy has been implemented in support vector-based approach (Dr. Cancer) that uses sequence profile and cancer-specific functional terms [89]. Dr. Cancer has been tested on 3,163 cancer-drive mutations from 74 proteins. The results show that cancer-related proteins are enriched for particular Gene Ontology terms that can be used to discriminate between cancer and other phenotypes. The same dataset was previously used to train CHASM a random forest-based algorithm [90]. CHASM was tested on ~600 missense mutations in glioblastoma multiforme achieving better performance than PolyPhen and previously developed cancer specific approaches. More recently, new methods such as MutationAssessor [78], TransFIC [93] and FATHMM [92] have been optimized to detect cancer driver variants. MutatationAssessor uses evolutionary information patterns calculating an entropy-based functional impact score from homolog and paralog proteins. The method trained on ~19,000 disease associated variants has been tested on 10,000 mutations in COSMIC database prioritized according to their recurrence and multiplicity. The authors estimated that ~5% of the cancer-relevant mutation involves a change in the protein function rather than standard loss and gain of function mechanism events. In TransFIC approach the prediction of cancer drivers mutation has been tuned considering subset of variants associated to the same Gene Ontology term and selecting a variable threshold for discriminating driver from passenger mutations. The performance of the optimized functional impact score has been tested on subset of the COSMIC mutations classified according to their recurrence. The results show that groups of protein with different function, posses distinct baseline tolerance to deleterious mutations. Finally, FATHMM is a Hidden Markov Model based algorithm using protein domain information that has been optimized for predicting cancer causing mutations. The method, tested on previously collected datasets result in improved performances with respect to previously developed methods. Although the algorithms based on functional information (Gene Ontology terms/ protein domains) achieve better results than standard conservation-based approaches, a fair testing procedure is more difficult because the predictions can be biased toward more abundant functional classes.

Analysis of cancer gene pathways and networks

The accumulation of somatic mutations during the lifespan is the main cause of cancer. Several identified somatic mutations occur in genes involved in many signaling, regulatory and metabolic pathways. Indeed mutated genes such as TP53 and PI3KCA are hubs in pathways and interaction networks which control cell proliferation, growth and apoptosis [95, 96]. In addition, recent sequencing studies [4, 38] revealed that cancer driver genes tend to cluster within a limited number of essential pathways, and rarely mutations on genes in the same pathway co-occur in the same patient. This mutually exclusive genomic events have been observed, for example, in the case of BRAF and KRAS (involved in RAS/RAF pathway) in colorectal cancer [97], APC and CTNNB1 (involved in beta catenin pathway) [98], EGFR and KRAS in lung adenocarcinomas [99], TP53 and MDM2 in many different cancer types [100].

From these observations, it is evident that the analysis of genomic variations across gene pathways and networks is important to understand the combinatorial effect of the mutations and explain the disease mechanism. Pathway and network analysis can be performed using previously annotated gene pathways or testing alternative routes from gene interaction networks. Reference databases collecting information about gene pathways are the Kyoto Encyclopedia of Genes and Genomes Kyoto Encyclopedia of Genes and Genomes (KEGG) [65] and the Molecular Signatures Database (MSigDB) [101]. Biological interaction networks mainly consist of protein-protein interaction data, which are collected in databases such as IntAct [102] and iRefIndex [103]. Another important database is the Reactome, which aggregates data on protein-protein interaction networks and gene pathways [104].

Recently, several methods for the analysis of cancer gene pathway have been developed [105–115]. Few examples are PathScan [114], which is an annotation-based approach, and HotNet2 [115], MDPFinder [107], MEMo [105], and Dendrix [113] which are able to identify driver pathways. PathScan [114] is a probabilistic model that takes into account the length of the genes and differences in their mutation rates under the null hypothesis. The method combines single-sample p-values using an integral approximation that estimates a pathway-specific overall p-value. PathScan represents an alternative to previous approaches, which reduce a gene set into a unique gene simply combining the total mutations and mutable positions of each gene. PathScan was applied to the analysis of data from lung adenocarcinoma sequencing project to find significantly mutated 129 KEGG pathways. The results show that PathScan was able to identify significant pathways detected in previous studies. In addition PathScan found the focal adhesion pathway to be significant in agreement with previous expression studies on prostate and ovarian cancers [114]. HotNet2 [115], has been recently developed for identifying significantly mutated groups of interacting genes from large cancer sequencing studies. It uses a heat diffusion model that encodes both the topology of the network and the significance of the observed frequencies of each mutated gene. HotNet2 has been tested on 3,281 samples from 12 cancer types in TCGA studies. In this analysis it was able to identify 16 significantly mutated subnetworks that comprise well-known cancer signaling pathways. Among them, the well-known TP53, phosphoinositide 3-kinase (PI3K), NOTCH and receptor tyrosine kinases (RTK) signaling pathways.

Although annotation-based methods are able to prioritize important cancer pathways, they present some limitations. The main limitations are the incompleteness of pathway databases and the presence of multiple overlapping gene sets across pathways. Thus, more general methods for detecting significantly mutated gene sets use information from biological interaction networks. The Mutual Exclusivity Modules in Cancer (MEMo) algorithm [105] integrates copy number alteration and mutation data and maps them into biological networks. The method uses correlation analysis and statistical tests to identify network modules of genes recurrently altered across a set of tumor samples, participate in the same biological process, and alteration events are mutually exclusive. MEMo has been tested on a set of ~400 samples from glioblastoma and ovarian cancer. In the first cancer study the method was able to identify important signaling modules such as p53 and PI(3)K. In ovarian cancer, it detected mutually exclusive variation events between BRCA and genes in the Rb module.

The De novo Driver Exclusivity (Dendrix) algorithm [113] is a tool for discovery of mutated driver pathways in cancer using only mutation information from the cancer samples. The method introduces the concepts of coverage and exclusivity to distinguish group of genes with driver mutations from set of genes with passenger mutations. The Dendrix algorithm, which has been applied to the analysis of different cancer types in TCGA study, was able to identify 8 mutation groups mutated in 94% of patients from 17 cancer types and groups of mutually exclusive genes in lung adenocarcinoma and glioblastoma multiforme. Similarly, MDPFinder [107] identifies de novo mutated driver pathways from mutation data solving the maximum weight submatrix problem. The authors implemented an exact method based on binary linear programming and genetic algorithm to combine mutation and expression data. MDPFinder has been tested on a set of ~500 samples from head and neck squamous cell carcinomas glioblastoma and ovarian carcinoma. The results show that the integrative model based on mutation and expression data was able to identify biologically relevant gene sets detected in previous studies [107].

Finally, PARADIGM[116] estimates patient-specific genetic activities incorporating curated gene interactions from the NCI pathway interaction database (PID) [117]. The method, that can use many types of omics data, represents a gene by a factor graph with associated known activity and expression. PARADIGM, tested on breast cancer and glioblastoma samples, was able to identify altered activities in cancer-related pathways with less false-positives than a previous developed method.

An interesting approach, which applies network analysis to the study of the cancer genome, is the Network-Based Stratification (NBS) method [118]. NBS is an algorithm for the classification of cancer subtypes that clusters patients with somatic mutations in similar gene network regions. More in details, the NBS represents each patient with a profile of somatic mutated genes mapped on the human gene interaction network. After a smoothing procedure, the patients are clustered using a non-negative matrix factorization-based approach. The method performs a second clustering procedure to find subgroups of patients frequently co-stratified after random sampling. NBS has been tested on a set of ~1,000 samples from ovarian, uterine endometrial and lung tumors and used to determine relationship between cancer subtypes and patient survival.

Classification of cancer subclonal variants

Tumor progression is an evolutionary process which starts from a single cell and results in the selection of more aggressive subclones [119]. The presence of different cell populations affects the accurate detection of somatic mutations in cancer sample. Although all subclones in a sample have the same origin, they are present in different proportions and have heterogeneous patterns of somatic mutations. Thus, somatic mutations in cancer cells present in small proportion are more difficult to detect because correspond to variations supported by a low number of reads, further reduced by impurity in the tumor sample. In a recent work [13], Ding and colleagues estimated that approximately 340X coverage is needed for detecting (with 99% of chance) at least 3 reads when a heterozygous variant is present in a 5% subclonal population.

A possible solution to this issue lies in single-cell sequencing for revealing diversities in the pattern of somatic mutations within a tumor [120, 121]. Although this technique provides a better characterization of the mutations occurring in each subclonal population, the approach is affected by errors introduced during the amplification process. Furthermore, single-cell experiments need to be performed on multiple tumor cells to have a general overview of the mutation pattern. The presence of subclones in the cell population can be detected analyzing the distribution of reads or variant allele fraction for the somatic mutations. If the distribution is multimodal, the presence of subclonal populations is expected. Figure 7 shows the distribution of variant allele fraction (VAF) in a lung adenocarcinoma cancer sample from TCGA. Although the distribution of variant allele fraction provides an indication about the clonal architecture, the analysis of this data is affected by sample impurity and copy number alterations. In the last few years many computational tools have been developed to estimate the presence of subclonal populations [122–127], some of them are discussed below.

ABSOLUTE [122] uses data from copy number variations for optimizing models of recurrent cancer karyotypes and expected allelic fraction values for somatic SNVs. These models are then used to re-extract information about the absolute cellular copy number of local DNA segments and the number of mutated alleles for somatic SNVs. The ABSOLUTE algorithm has been applied to the analysis of copy-number profiles from 3,155 cancer samples, identifying recurrent genome doubling events that influence tumor progression. ABSOLUTE was also used for the analysis of exome sequencing data from 214 ovarian carcinoma tumor-normal pairs. The method was able to identify large subclone populations with predominant somatic mutations and small subset of subclones with heterozygotes mutations in tumor suppressor genes TP53 and NF1 and in a candidate tumor suppressor gene CDK12.

Alternative methods [124, 125] use phylogenetic approaches to study the evolution of tumor cells. For example, PyClone [125] assumes that clonal population follows a perfect and persistent phylogeny. Under this assumption, for which each somatic mutation only happen once in the evolutionary history and not reverse mutation are allowed, subclonal cells can be identified and their prevalence estimated. PyClone uses a Bayesian clustering algorithm to group sets of somatic mutations belonging to the same cluster of subclonal cell accounting for allelic imbalances introduced by copy-number alterations and normal-cell impurity. The method has been tested on a simulated dataset produced from mixtures of DNA extracted from four 1000 Genomes Project sample and four spatially separated samples from a primary, untreated high-grade ovarian tumor. The results show that PyClone outperforms two genotype-naïve methods using binomial and infinite beta-binomial mixture models. Similarly, PhyloSub [124] infers the phylogeny and genotype of the major lineages in the clonal population calculating the Bayesian prior over the trees clustering the SNVs. A sampling procedure is use to find the optimized the join posterior distribution with higher probability to generate the observed frequencies of somatic mutations. PhyloSub has been tested using deep exome sequencing data from acute myeloid leukemia and chronic lymphocytic leukemia. The results show that PhyloSub is able to identify both linear and branching subclonal lineages.

A recent publication [123] presents a novel combinatorial approach based on binary partition tree (BPT) to model the mechanism of clonal expansion in tumor and estimate the subpopulations of tumor using the variant allele frequencies of somatic mutations. The authors demonstrated that finding a BPT is a NP-complete problem, and derived an approximation algorithm for an optimized version of the problem. Finally they implemented a recursive approach that finds the solution of the optimized BPT problem in a polynomial time. The developed algorithm can detect errors in the estimation of the variant allele frequencies of somatic mutations, which cannot be correctly estimated because of the admixture of normal cells in the tumor sample. The performance of BPT algorithm has been tested on simulated and real cancer data showing it generates more consistent results, and it is faster than previously developed methods.

The methods for the detection of genetic heterogeneity in cancer can be also used to detect subclonal mutation conferring drug resistance. This idea has been investigated using cloneHD [126], a new subclonal reconstruction algorithm optimized using both information about somatic mutations and correlated changes generated by copy-number changes. The method has been applied to the analysis of sequencing data from time-resolved samples from breast cancer and of chronic lymphocytic leukemia. The results demonstrate that cloneHD can be a valuable tool for tracking cancer development and monitoring the response of a patient to therapy regimens.