The dNdScv R package is a suite of maximum-likelihood dN/dS methods designed to quantify selection in cancer and somatic evolution (Martincorena et al., 2017). The package contains functions to quantify dN/dS ratios for missense, nonsense and essential splice mutations, at the level of individual genes, groups of genes or at whole-genome level. The dNdScv method was designed to detect cancer driver genes (i.e. genes under positive selection in cancer) on datasets ranging from a few samples to thousands of samples, in whole-exome/genome or targeted sequencing studies.
The background mutation rate of each gene is estimated by combining local information (synonymous mutations in the gene) and global information (variation of the mutation rate across genes, exploiting epigenomic covariates), and controlling for the sequence composition of the gene and mutational signatures. Unlike traditional implementations of dN/dS, dNdScv uses trinucleotide context-dependent substitution matrices to avoid common mutation biases affecting dN/dS (Greenman et al., 2006).
This vignette shows how to perform driver discovery and selection analyses with dNdScv in cancer sequencing data. The current version of dNdScv only supports human data, but future versions will incorporate functions to allow studies of selection on any species.
Main reference: Martincorena I, et al. (2017) Universal Patterns of Selection in Cancer and Somatic Tissues. Under review.
library("seqinr")
library("Biostrings")
library("MASS")
library("GenomicRanges")
library("dndscv")## Warning: replacing previous import 'Biostrings::translate' by
## 'seqinr::translate' when loading 'dndscv'data("dataset_simbreast", package="dndscv")
dndsout = dndscv(mutations)## Warning in dndscv(mutations): Same mutations observed in different
## sampleIDs. Please verify that these are independent events and remove
## duplicates otherwise.For this example, we have used a simulated dataset of 196 breast cancer exomes provided in the package. The simplest way to run the dndscv function is to provide a table of mutations. Mutations are provided as a data.frame with five columns (sampleID, chromosome, position, reference base and mutant base). It is important that only unique mutations are provided in the file. Multiple instances of the same mutation in related samples (for example, when sequencing multiple samples of the same tumour) should only be listed once.
head(mutations)## sampleID chr pos ref mut
## 1 Sample_1 1 871244 G C
## 2 Sample_1 1 6648841 C G
## 3 Sample_1 1 17557072 G A
## 4 Sample_1 1 22838492 G C
## 5 Sample_1 1 27097733 G A
## 6 Sample_1 1 27333206 G AWith the example dataset provided, the function should take about one minute to run. In this example, the function issues a warning as it detects the same mutation in more than one sample, requesting the user to verify that the input table of mutations only contains independent mutation events. In this case, each sample corresponds to a different patient and so the warning can be safely ignored.
We have run dNdScv with default parameters. This includes removing ultra-hypermutator samples and subsampling mutations when encountering too many mutations per gene in the same sample. These were designed to protect against loss of sensitivity from ultra-hypermutators and from clustered artefacts in the input mutation table, but there are occasions when the user will prefer to relax these (see Example 2).
The output of the dndscv function is a list of objects. For an analysis of exome or genome data, the most relevant output will often be the result of neutrality tests at gene level. P-values for substitutions are obtained by Likelihood-Ratio Tests as described in (Martincorena et al, 2017) and q-values are obtained by Benjamini-Hodgberg’s multiple testing correction. The table also includes information on the number of substitutions of each class observed in each gene, as well as maximum-likelihood estimates (MLEs) of the dN/dS ratios for each gene, for missense (wmis), nonsense (wnon), essential splice site mutations (wspl) and indels (wind).
sel_cv = dndsout$sel_cv
print(head(sel_cv), digits = 3)## gene_name n_syn n_mis n_non n_spl n_ind wmis_cv wnon_cv wspl_cv
## 18056 TP53 1 43 5 4 5 113.68 221.7 221.7
## 12976 PIK3CA 3 34 0 0 3 29.76 0.0 0.0
## 9224 KRAS 1 21 0 0 0 125.99 0.0 0.0
## 18923 VHL 3 9 1 0 4 24.76 38.2 38.2
## 1296 APC 2 8 10 0 6 2.82 31.2 31.2
## 1465 ARID1A 1 7 10 0 3 3.30 53.5 53.5
## wind_cv pmis_cv ptrunc_cv pallsubs_cv qallsubs_cv pind_cv
## 18056 138.8 0.00e+00 1.11e-16 0.00e+00 0.0e+00 1.73e-09
## 12976 30.7 0.00e+00 5.67e-01 0.00e+00 0.0e+00 2.12e-04
## 9224 0.0 0.00e+00 8.43e-01 0.00e+00 0.0e+00 1.00e+00
## 18923 204.5 7.43e-09 2.23e-02 7.95e-09 2.0e-05 1.32e-08
## 1296 23.1 4.96e-02 2.82e-10 2.12e-09 6.1e-06 1.81e-06
## 1465 14.4 3.64e-02 2.96e-12 2.39e-11 1.2e-07 1.84e-03
## pglobal_cv qglobal_cv
## 18056 0.00e+00 0.00e+00
## 12976 0.00e+00 0.00e+00
## 9224 0.00e+00 0.00e+00
## 18923 4.00e-15 2.01e-11
## 1296 1.31e-13 5.28e-10
## 1465 1.40e-12 4.68e-09signif_genes = sel_cv[sel_cv$qglobal_cv<0.1, c("gene_name","qglobal_cv")]
rownames(signif_genes) = NULL
print(signif_genes)## gene_name qglobal_cv
## 1 TP53 0.000000e+00
## 2 PIK3CA 0.000000e+00
## 3 KRAS 0.000000e+00
## 4 VHL 2.007394e-11
## 5 APC 5.281677e-10
## 6 ARID1A 4.679831e-09
## 7 PTEN 2.403743e-08
## 8 IDH1 8.742613e-06
## 9 PBRM1 3.169391e-05
## 10 RB1 9.215911e-04
## 11 PIK3R1 9.568278e-04
## 12 CDKN2A 9.568278e-04
## 13 BRAF 4.821497e-03
## 14 KMT2C 2.652243e-02
## 15 FBXW7 3.287783e-02
## 16 EGFR 5.088073e-02Note in the table that the dN/dS ratios of significant cancer genes are typically extremely high, often >10 or even >100. This contains information about the fraction of mutations observed in a gene that are genuine drivers under positive selection. For example, for a highly significant gene, a dN/dS of 10 indicates that there are 10 times more non-synonymous mutations in the gene than neutrally expected, suggesting that at least around 90% of these mutations are genuine drivers (Greenman et al, 2006; Martincorena et al, 2017).
Another output that can be of interest is a table with the global MLEs for the dN/dS ratios across all genes. dN/dS ratios with associated confidence intervals are calculated for missense, nonsense and essential splice site substitutions separately, as well as for all non-synonymous substitutions (wall) and for all truncating substitutions together (wtru), which include nonsense and essential splice site mutations.
print(dndsout$globaldnds)## name mle cilow cihigh
## wmis wmis 1.0222179 0.9855147 1.060288
## wnon wnon 1.0616564 0.9813078 1.148584
## wspl wspl 0.9707887 0.8656096 1.088748
## wtru wtru 1.0326258 0.9651248 1.104848
## wall wall 1.0241562 0.9879328 1.061708Global dN/dS ratios in somatic evolution of cancer, and seemingly of healthy somatic tissues, appear to show a near-universal pattern of dN/dS~1, with exome-wide dN/dS ratios typically slightly higher than 1 (Martincorena et al., 2017). Only occasionally, I have found datasets with global dN/dS<1, but upon closer examination, this has typically been caused by contamination of the catalogue of somatic mutations with germline SNPs. An exception are melanoma tumours, which show a bias towards slight underestimation of dN/dS due to the signature of ultraviolet-induced mutations extending beyond the trinucleotide model (Martincorena et al., 2017). In my personal experience, datasets of somatic mutations with global dN/dS<<1 have always reflected a problem of SNP contamination or an inadequate substitution model, and so the evaluation of global dN/dS values can help identify problems in certain datasets.
The dndscv function also outputs other results that can be of interest, such as an annotated table of coding mutations (annotmuts), MLEs of mutation rate parameters (mle_submodel), lists of samples and mutations excluded from the analysis and a table with the observed and expected number of mutations per gene (genemuts), among others.
head(dndsout$annotmuts)## sampleID chr pos ref mut geneind gene ref_cod mut_cod strand
## 1 Sample_1 1 871244 G C 15262 SAMD11 G C 1
## 2 Sample_1 1 6648841 C G 19346 ZBTB48 C G 1
## 3 Sample_1 1 17557072 G A 12419 PADI1 G A 1
## 4 Sample_1 1 22838492 G C 19338 ZBTB40 G C 1
## 5 Sample_1 1 27097733 G A 1465 ARID1A G A 1
## 6 Sample_1 1 27333206 G A 6236 FAM46B C T -1
## ref3_cod mut3_cod aachange ntchange impact pid
## 1 AGA ACA R133T G398C Missense ENSP00000342313
## 2 TCA TGA V569V C1707G Synonymous ENSP00000366902
## 3 AGA AAA E353E G1059A Synonymous ENSP00000364620
## 4 GGA GCA D776H G2326C Missense ENSP00000364798
## 5 TGA TAA E1108K G3322A Missense ENSP00000320485
## 6 TCA TTA L169L C507T Synonymous ENSP00000289166dNdScv relies on a negative binomial regression model across genes to refine the estimated background mutation rate for a gene. This assumes that the variation of the mutation rate across genes that remains unexplained by covariates or by the sequence composition of the gene can be modelled as a Gamma distribution. This model typically works well on clean cancer genomic datasets, but not all datasets may be suitable for this model. In particular, very low estimates of \(\theta\) (the overdispersion parameter), particularly \(\theta<1\), may reflect problems with the suitability of the dNdScv model for the dataset.
print(dndsout$nbreg$theta)## [1] 5.654063An additional set of neutrality tests per gene are performed using a more traditional dN/dS model in which the local mutation rate for a gene is estimated exclusively from the synonymous mutations observed in the gene (dNdSloc) (Wong, Martincorena, et al., 2014). This test is typically only powered in very large datasets. For example, in the dataset used in this example, comprising of 196 simulated breast cancer exomes, this model only detects ARID1A as significantly mutated.
signif_genes_localmodel = as.vector(dndsout$sel_loc$gene_name[dndsout$sel_loc$qall_loc<0.1])
print(signif_genes_localmodel)## [1] "ARID1A"The dndscv function can take a list of genes as an input to restrict the analysis of selection. This is strictly required when analysing targeted sequencing data, and might also be used to obtain global dN/dS ratios for a particular group of genes.
To exemplify the use of the dndscv function on targeted data, we can use another example dataset provided with the dNdScv package:
library("seqinr")
library("Biostrings")
library("MASS")
library("GenomicRanges")
library("dndscv")
data("dataset_normalskin", package="dndscv")
data("dataset_normalskin_genes", package="dndscv")
dndsskin = dndscv(m, gene_list=target_genes, max_muts_per_gene_per_sample = Inf, max_coding_muts_per_sample = Inf)## Warning in dndscv(m, gene_list = target_genes, max_muts_per_gene_per_sample
## = Inf, : Mutations observed in contiguous sites within a sample. Please
## annotate or remove dinucleotide or complex substitutions for best results.## Warning in dndscv(m, gene_list = target_genes, max_muts_per_gene_per_sample
## = Inf, : Same mutations observed in different sampleIDs. Please verify that
## these are independent events and remove duplicates otherwise.This dataset comprises of 3,408 unique somatic mutations detected by ultra-deep (~500x) targeted sequencing of 74 cancer genes in 234 small biopsies of normal human skin (epidermis) from four healthy individuals. Note that all of the mutations listed in the input table are genuinely independent events and so, again, we can safely ignore the two warnings issued by dndscv. For more details on this study see:
Martincorena I, et al. (2015) High burden and pervasive positive selection of somatic mutations in normal human skin. Science. 348(6237):880-6. doi: 10.1126/science.aaa6806.
In the paper above, we described a strong evidence of positive selection on somatic mutations occurring in normal human skin throughout life. These mutations are detected as microscopic clones of mutant cells in normal skin. The dNdScv analysis below recapitulates some of the key analyses of selection in this study:
sel_cv = dndsskin$sel_cv
print(head(sel_cv[sel_cv$qglobal_cv<0.1,c(1:10,16,17)]), digits = 3)## gene_name n_syn n_mis n_non n_spl n_ind wmis_cv wnon_cv wspl_cv wind_cv
## 47 NOTCH1 20 224 51 32 61 6.00 25.89 25.89 16.853
## 48 NOTCH2 8 106 21 9 25 5.46 16.14 16.14 7.142
## 26 FAT1 20 68 26 8 18 1.53 10.06 10.06 2.770
## 72 TP53 9 43 11 7 7 3.28 13.29 13.29 12.546
## 49 NOTCH3 15 35 8 2 11 1.13 3.79 3.79 3.345
## 62 RBM10 6 14 8 2 1 1.28 7.23 7.23 0.759
## pglobal_cv qglobal_cv
## 47 0.00e+00 0.00e+00
## 48 0.00e+00 0.00e+00
## 26 1.11e-16 2.74e-15
## 72 2.95e-13 5.46e-12
## 49 6.51e-04 9.64e-03
## 62 1.32e-03 1.62e-02print(dndsskin$globaldnds, digits = 3)## name mle cilow cihigh
## wmis wmis 1.29 1.18 1.42
## wnon wnon 2.00 1.69 2.36
## wspl wspl 2.69 2.13 3.39
## wtru wtru 2.18 1.88 2.52
## wall wall 1.37 1.25 1.50Classic maximum-likelihood implementations of dN/dS use a simple substitution model with a single rate parameter. Mutations are classified as either transitions (C<>T, G<>A) or transversions, and the single rate parameter is a transition/transversion (ts/tv) ratio reflecting the relative frequency of both classes of substitutions (Goldman & Yang, 1994). The dndscv function can take a different substitution model as input. The user can choose from existing substitution models provided in the data directory as part of the package or input a different substitution model as a matrix:
library("dndscv")
# 192 rates (used as default)
data("submod_192r_3w", package="dndscv")
colnames(substmodel) = c("syn","mis","non","spl")
head(substmodel)## syn mis non spl
## [1,] "t*AAA>ACA" "t*AAA>ACA*wmis" "t*AAA>ACA*wnon" "t*AAA>ACA*wspl"
## [2,] "t*AAA>AGA" "t*AAA>AGA*wmis" "t*AAA>AGA*wnon" "t*AAA>AGA*wspl"
## [3,] "t*AAA>ATA" "t*AAA>ATA*wmis" "t*AAA>ATA*wnon" "t*AAA>ATA*wspl"
## [4,] "t*AAC>ACC" "t*AAC>ACC*wmis" "t*AAC>ACC*wnon" "t*AAC>ACC*wspl"
## [5,] "t*AAC>AGC" "t*AAC>AGC*wmis" "t*AAC>AGC*wnon" "t*AAC>AGC*wspl"
## [6,] "t*AAC>ATC" "t*AAC>ATC*wmis" "t*AAC>ATC*wnon" "t*AAC>ATC*wspl"# 12 rates (no context-dependence)
data("submod_12r_3w", package="dndscv")
colnames(substmodel) = c("syn","mis","non","spl")
head(substmodel)## syn mis non spl
## [1,] "t*A>C" "t*A>C*wmis" "t*A>C*wnon" "t*A>C*wspl"
## [2,] "t*A>G" "t*A>G*wmis" "t*A>G*wnon" "t*A>G*wspl"
## [3,] "t*A>T" "t*A>T*wmis" "t*A>T*wnon" "t*A>T*wspl"
## [4,] "t*A>C" "t*A>C*wmis" "t*A>C*wnon" "t*A>C*wspl"
## [5,] "t*A>G" "t*A>G*wmis" "t*A>G*wnon" "t*A>G*wspl"
## [6,] "t*A>T" "t*A>T*wmis" "t*A>T*wnon" "t*A>T*wspl"# 2 rates (classic ts/tv model)
data("submod_2r_3w", package="dndscv")
colnames(substmodel) = c("syn","mis","non","spl")
head(substmodel)## syn mis non spl
## AAA>ACA "t" "t*wmis" "t*wnon" "t*wspl"
## AAA>AGA "t*tstv" "t*tstv*wmis" "t*tstv*wnon" "t*tstv*wspl"
## AAA>ATA "t" "t*wmis" "t*wnon" "t*wspl"
## AAC>ACC "t" "t*wmis" "t*wnon" "t*wspl"
## AAC>AGC "t*tstv" "t*tstv*wmis" "t*tstv*wnon" "t*tstv*wspl"
## AAC>ATC "t" "t*wmis" "t*wnon" "t*wspl"We can fit a traditional ts/tv model to the skin dataset using the code below:
library("seqinr")
library("Biostrings")
library("MASS")
library("GenomicRanges")
library("dndscv")
data("dataset_normalskin", package="dndscv")
data("dataset_normalskin_genes", package="dndscv")
dndsskin_2r = dndscv(m, gene_list=target_genes, max_muts_per_gene_per_sample = Inf, max_coding_muts_per_sample = Inf, sm = "2r_3w")## Warning in dndscv(m, gene_list = target_genes, max_muts_per_gene_per_sample
## = Inf, : Mutations observed in contiguous sites within a sample. Please
## annotate or remove dinucleotide or complex substitutions for best results.## Warning in dndscv(m, gene_list = target_genes, max_muts_per_gene_per_sample
## = Inf, : Same mutations observed in different sampleIDs. Please verify that
## these are independent events and remove duplicates otherwise.print(dndsskin_2r$mle_submodel)## name mle cilow cihigh
## 1 t 0.001138919 0.001032197 0.001256676
## 2 tstv 4.034639837 3.734319602 4.359112328
## 3 wmis 1.249162834 1.142169906 1.366178339
## 4 wnon 2.847342924 2.438120672 3.325250396
## 5 wspl 2.182462999 1.766235543 2.696777764sel_cv = dndsskin_2r$sel_cv
print(head(sel_cv[sel_cv$qglobal_cv<0.1, c(1:10,16,17)]), digits = 3)## gene_name n_syn n_mis n_non n_spl n_ind wmis_cv wnon_cv wspl_cv wind_cv
## 47 NOTCH1 20 224 51 32 61 5.07 25.93 25.93 16.853
## 48 NOTCH2 8 106 21 9 25 5.54 20.15 20.15 7.142
## 26 FAT1 20 68 26 8 18 1.44 13.43 13.43 2.770
## 72 TP53 9 43 11 7 7 3.25 13.69 13.69 12.546
## 49 NOTCH3 15 35 8 2 11 1.09 4.25 4.25 3.345
## 62 RBM10 6 14 8 2 1 1.15 7.14 7.14 0.759
## pglobal_cv qglobal_cv
## 47 0.00e+00 0.00e+00
## 48 0.00e+00 0.00e+00
## 26 0.00e+00 0.00e+00
## 72 2.03e-13 3.75e-12
## 49 2.89e-04 4.27e-03
## 62 9.92e-04 1.22e-02In general, the full trinucleotide model is recommended for cancer genomic datasets as it typically provides the least biased dN/dS estimates. The impact of using simplistic mutation models can be considerable on global dN/dS ratios (see Martincorena et al., 2017), and can lead to false signals of negative or positive selection. In general, the impact of simple substitution models on gene-level inferences of selection tends to be smaller.