Identifying cancer driver genes in tumor genome sequencing studies

2.1 Definition of P-values for identifying driver genes

For each gene, we test if the number of samples with ‘driver-like’ non-silent mutations is higher than that expected by the background mutation model M₀.

Let

then Y_ij is a Bernoulli random variable.

Define

Since we assume a different mutation rate for each sample, the probability s_ij varies across samples j = 1,…, J. It is calculated from the background mutation model which will be described in the Section 2.2.

We assign a score to every possible non-silent mutation according to its expected impact on the protein function: higher score for mutations with stronger impact. As will be shown, the order between scores rather than the actual scores determines the test statistics. Therefore, one can assign any score to each non-silent mutation to reflect the order of its impact on the protein function. We assign scores so that they comply with the following order: missense < inframe indel < mutation in splice sites < frameshift indel = nonsense. We also assign different scores to different types of missense mutations based on BLOSUM80 matrix, which is a matrix of scores for each of the 190 possible substitutions of the 20 standard amino acids.

Let T_ij be the maximum score of the non-silent mutations that occurred in sample j for gene i. If no mutation occurred, let T_ij = 0. Define F_ij(x) = P(T_ij < x|Y_ij = 1, M₀). We can obtain the distribution F_ij from the background mutation model described in the Section 2.2.

Then,

We use as our test statistic and define P-values by P(Z_i < z_i|M₀), where z_i is the observed value of Z_i. Z_i can be interpreted as a sum of mutated sample indicators weighted by . Since s_ij is larger than 0.5, is less than one and thus, the weights are negative. The larger the s_ij and F_ij(T_ij), the smaller the weight. Therefore, mutations with higher scores (stronger impact on the protein function) occurring in samples with low mutation rates (samples with large s_ij) contribute more in decreasing Z_i and thus, P-values.

We can generate the distribution of Z_i under the background mutation model by simulating Y_ij from Bernoulli(s_ij) and T_ij from F_ij for j = 1,…, J and then calculating the sum of . Then we can approximate the P-value P(Z_i < z_i|M₀) by computing the tail area of this background distribution beyond the observed value z_i for each gene i.

2.2 Background mutation model

The most distinguishing features of our background mutation model are that it does not assume separate mutation rates for non-silent and silent mutations and that it assumes separate mutation rates for different samples. We assume that each passenger mutation is generated from one background mutation rate process and that whether the mutation is non-silent or silent depends on the genetic code.

There are six types of mutations:

The transitions A:T → G:C and G:C → A:T change a purine to another purine or a pyrimidine to another pyrimidine. The transversions change a purine to a pyrimidine or vice versa. Because transitions occur more frequently than transversions, we assume separate mutation rates for transitions and transversions. Also it is generally observed that a mutation occurs more often at C:G than A:T and that a C:G appearing in CpG dinucleotides has a higher mutation rate than a C:G appearing in non-CpG dinucleotides. Therefore, we assume a separate background mutation rate for each combination of base pair types and CpG dinucleotides context.

We also assume that different tumor samples have different mutation rates. To keep the number of parameters manageable, we assume that relative frequencies of different types of mutations are same for each sample. Thus, the mutation rate in sample j for mutation type m is defined as the product of p_m, the ratio of mutation rate of the type m relative to the type 1 (A:T → G:C) and q_j, the mutation rate of the sample j for the mutation type 1 (Table 1).

To estimate the parameters in the background mutation model, we could fit the model in Table 1 to the sequences for which silent mutations were identified. (Most previous projects have evaluated silent mutations for only a subset of the genes.) To estimate the background mutation rate for insertions and deletions (indels), which are non-silent, however, we included in our estimation genes which have at most one non-silent mutation across all tumor samples; these genes are not likely to be related to tumorigenesis and thus the non-silent mutations in these genes are likely to be passenger mutations. However, since we selected these genes based on the total number of mutations occurring in each gene, the estimated background rates for these genes may be biased. Since the selection was based on the total number of mutations, it is unlikely that the relative frequencies of different types of mutations are subject to the bias, but the sample-specific mutation rates may be. Let q_j′ be the mutation rate of the sample j for mutation type 1 in the selected genes. Then we assume q_j′ = r · q_j, where r is the selection bias and q_j is the unbiased sample-specific mutation rate.

To estimate the parameters r, q_j and p₁,…, p₈, we first define the position of base pairs across all the sequenced genes. Since we assume that background mutation rates are independent of genes, we do not need to differentiate genes. Therefore, we concatenate all the sequenced genes and determine the position of each base pair from 1 to N, the total number of base pairs that are sequenced. Let K denote the subset of positions of the base pairs belonging to the genes used for silent mutation detection and let L denote the subset of positions of the base pairs belonging to the genes which have at most one non-silent mutation across all samples.

For position k in genes for which silent mutations have been evaluated, the probability that a silent transition of type i (second column in Table 1) occurs in sample j equals q_jc_kp_i where c_k is 1 if a transition at position k results in a silent mutation, otherwise c_k is 0. The probability that a silent transversion of type i occurs at that position equals q_jd_kp_i, where d_k is the number of silent transversions possible at position k (0, 1 or 2). The full set of probabilities definitions are shown in Table 2 based on the indicators X_jk, which indicate the type of mutation occurring at position k in sample j. Since any mutation is either a silent mutation or a non-silent mutation, c_k + e_k = 1 (number of possible transition mutations) and d_k + f_k = 2 (number of possible transversion mutations). When a mutation occurs within splice sites, it is considered to be non-silent, therefore c_k = d_k = 0, e_k = 1, f_k = 2 if k belongs to splice sites.

Table 2.

Definition of probabilities of X_jk

Open in a new tab

a_k = (c_kp_tk + d_kp_vk)I(k ∈ K)+r(e_kp_tk + f_kp_vk + p₇ + p₈)I(k ∈ L);

I(x), indicator function, 1 if x is true and 0 otherwise;

c_k, number of silent transitions possible at position k (0 or 1); d_k, number of silent transversions possible at position k (0, 1 or 2); e_k, number of non-silent transitions possible at position k (0 or 1); f_k, number of non-silent transversions possible at position k (0, 1 or 2); t_k, mutation type ID for the transition at position k (1, 3 or 5); v_k, mutation type ID for the transversion at position k (2, 4 or 6); non, no mutation; sts, silent transition; stv, silent transversion; nts, non-silent transition; ntv, non-silent transversion; iid, inframe indel; fid, frameshift indel.

All of the constants shown in Table 2 can be determined from the gene sequence and genetic code. However, the values are ambiguous in cases where genes have several alternative transcripts, and where some base pairs belong to different codons in alternative transcripts. We describe how to determine the values of the constants in such cases in the Supplementary Material.

2.3 Estimation of parameters

We use the method of moments to estimate r and p_m. The process of obtaining the method of moments estimates and is described in the Supplementary Material.

The estimation of q_j is more complex because the number of base pairs sequenced per sample is not sufficient to estimate the extremely small mutation rate q_j accurately. For example, no mutations were found for many samples in the data from Ding et al. (2008). Therefore, the maximum likelihood estimate of q_j would be zero for those samples, which are problematic point estimates. To improve the accuracy of the estimates, we use empirical Bayes methods to estimate the distribution of q_j. Empirical Bayes methods borrow information from all the samples for estimating each q_j, therefore give more robust estimates of q_j. (Casella, 1985)

We assume the prior distribution f of q_j is uniform on (α, β). As estimates of α, β, we use the values maximizing the marginal likelihood given the estimates, and .

The posterior distribution of q_j given the data X_jk for k ∈ K ∪ L, j = 1,…, J and the estimated parameters is

where the product is over all positions k, the probabilities are computed from the formulas in Table 2 and λ is the normalizing constant.

We use the posterior distribution of q_j in calculating s_ij=P(Y_ij=0|M₀) rather than using the point estimates of q_j to take into account the uncertainty in the point estimate. Therefore,

where G_i is the subset of positions of the base pairs belonging to the gene i and . The integration with regard to the posterior distribution of q_j is performed numerically. The resulting values of s_ij are used as described in Section 2.1 for computing statistical significance.

The distributions F_ij(x) = P(T_ij < x|Y_ij = 1, M₀) are also needed for the significance tests used to identify driver genes and are computed from:

where T′_jk is the score of the mutation occurring in position k and sample j. The distribution P(T′_jk < x|X_jk = nts, ntv, iid , or fid ) can be easily calculated from the genetic code and background mutation model. The process of the derivation of F_ij(x) is explained in the Supplementary Material.

3.1 Simulation study

We first performed a simulation study to evaluate our method. For the comparison with the method of Ding et al. (2008), we did not include the mutation score T_ij in the test statistics, that is, we use the test statistic instead of .

We generate simulated data based on the data of Ding et al. (2008). We first generate passenger mutations by shuffling the locations of all observed non-silent and silent mutations across the genes sequenced. There are 1013 non-silent mutations observed in 623 genes and 108 silent mutations observed in 250 genes. For these mutations, we change the base pair positions in which the mutation occurred as follows: we randomly sample the base pair positions from the base pair positions within the sequence of all genes, which correspond to the same base pair types as the mutations. If a mutation occurred in the base pair A, we sample its new base pair position from all the base pair positions within the sequence of all genes corresponding to a base pair A. If the base pair is G or C, we also restrict the sampling by the CpG dinucletoide context. We then determine which of these mutations are non-silent or silent according to the genetic code. Since we randomly sample the base pair positions of all the mutations, they become evenly spread across all genes.

To see the effect of variation of mutation rates across samples, we change the sample ID in which mutations occurred by sampling a new sample ID under two different distributions namely, moderate sample variation and high sample variation.

The first distribution, moderate sample variation, is estimated from the background mutations of the data from Ding et al. (2008). We sample each sample ID with the probability proportional to the number of passenger mutations (silent mutations and non-silent mutations observed in genes with at most one non-silent mutations) that occurred in the sample. For the second distribution, high sample variation, we increase the mutation rates of the 10 samples with highest mutation rate by a factor of 10.

Finally, we make 20 driver genes by adding five non-silent mutations to 20 selected genes.

In our simulations, we have used the true expected ratio of non-silent to silent mutations (R) in applying the method of Ding et al. (2008) because we did not have their software for estimating R. This may somewhat overestimate the accuracy of their method.

Each simulation was repeated for 200 replications. The average number of true and false positive driver genes claimed based on P-value cutoffs of 0.005 and 0.01, respectively, are shown in Table 3. Our method finds more true positives and fewer false positives than the method of Ding et al. (2008). We did Wilcoxon signed rank test of the null hypothesis that the distribution of number of true or false positives from both methods are same and presented one-sided P-values in the last column of Table 3. For moderate sample variation, the P-values for false positives are 0.0001 and 0.0008, and the P-values for true positives are less than 10⁻¹⁶. For high sample variation, all the P-values are less than 10⁻¹⁶. This shows that the difference in the number of true positives or false positives between two methods is significant and it gets more significant as the sample variation grows larger.

3.2 Results for the data of Ding et al. (2008

We identified 28 genes as driver genes with the false discovery rate (FDR) controlled at 5% using the Benjamini and Hochberg method. These include EGFR, CDKN2A, KRAS, STK11, TP53, EPHA3, NF1, ATM, RB1, APC, INHBA, ERBB4, PTPRD, FGFR4, PTEN, EPHA5, NTRK3, NTRK1, KDR, LRP1B, PAK3, NRAS, LTK, ZMYND10, EPHA7, MYO3B, NTRK2 and TFDP1. The P-values of the selected genes are given in Table 4.

By the method of Ding et al. (2008), 22 genes were found to be significant at the 5% FDR level. These genes include two genes that we do not find significant. Our method finds eight genes, which they do not find significant. We drew a map of the genes selected by each method versus tumor samples in Figure 2. The genes are ordered by the P-value obtained by our method and samples are ordered according to the total number of genes with non-silent mutation (among all 623 genes). The genes indicated by the red bar are those which Ding et al. (2008) find significant but we do not. The genes indicated by the yellow bar are those which we find significant but they do not. The genes indicated by the blue bar are those which both methods find significant.

For most of the well known oncogene and tumor suppressor genes in the dataset we analyzed (EGFR, CDKN2A, KRAS, STK11, TP53, NF1, RB1 PTEN, NRAS), the mutations occurring in them have very high mutation scores. For example, in STK11, most mutations are frameshift indels, nonsense or mutations in splice sites. Even the missense mutations represent poorly conserved amino acid changes. In RB1, all seven mutations that occurred in the gene are either frameshift indel, nonsense or mutations in splice sites. This is consistent with our scoring system that mutations in driver genes will tend to have strong impact on protein functions. By incorporating mutation scores in calculating P-values, driver genes have smaller P-values and thus are better identified. For the well- known driver genes EGFR, CDKN2A, KRAS, STK11, TP53 which already have computed P-values of zero due to frequent mutations, the effect of scoring makes little difference. But for genes like NF1, RB1, APC, INHBA, ERBB4, FGFR4, PTEN and NRAS, their P-values are about one-third on average of the P-values calculated without incorporating the mutation scores. Most of those genes are well-known cancer driver genes, and incorporating mutation scores helps in their identification.

Not having mutations with high scores does not preclude a gene being a potential driver gene, but these cases tend to be infrequently mutated genes that occur in samples with large mutation rates as was the case for CDH11 and PDGFRA. These two genes are selected to be significant by Ding's method, but not by our method. Figure 2 shows that the mutations in these genes are clustered in the highly mutated samples except one in PDGFRA. However, most mutations in PDGFRA are missense mutations with low mutation scores, offsetting the effect of low mutation rate to the test statistics of PDGFRA. Also, some of the mutations in both genes occur in the same sample, therefore, our method which is based on the number of samples with mutations rather than the total number of mutations assigned larger P-values to them than the method of Ding et al. (2008).

There are eight genes beside the yellow bar. The gene with the smallest P-value is PTEN, a well-known tumor suppressor gene. Ding et al. (2008) did not find it significant because the total number of mutations is so small (four). However, since each of the four mutations occurred in different samples with low mutation rates and the score of each mutation is high (one nonsense mutation and three missense mutation with high score), our method could find it significant.

The gene with the second smallest P-value is NRAS, a well- known oncogene. Although there were only three total mutations in this gene, all of them are the same missense mutation changing glutamine to leucine, which has a high score. Also, one of the mutation occurred in a sample with low mutation rate, thus we could find it significant.

The gene ZMYND10 is a candidate tumor suppressor gene whose association with carcinomas is suggested by Agathanggelou et al. (2003); Cho (2007); Lerman and Minna (2000); Marsit et al. (2005); Qiu et al. (2004).

The gene EPHA7 is a member of the ephrin receptor family and is known to be related to oncogenesis (Kiyokawa et al., 1994). The other two members EPHA3 and EPHA5 are also selected to be significant by both methods, implying that EPHA7 is potentially involved in oncogenesis.

NTRK2 is a member of the neurotrophic tyrosine receptor kinase (NTRK) family, which phosphorylates members of the MAPK pathway. It is known to be potentially implicated in oncogenesis (Marchetti et al., 2008) and also the other two members of the NTRK receptor family, NTRK1 and NTRK3 are selected to be significant by both methods, supporting the implication of NTRK2 in oncogenesis.

TFDP1 is a transcription factor and its overexpression or amplification is known to be associated with carcinomas (Melchor et al., 2009; Yasui et al., 2003). The role of LTK and MYO3B in oncogenesis is not well known.