SNVHMM: predicting single nucleotide variants from next generation sequencing

Abstract

Background

The rapid development of next generation sequencing (NGS) technology provides a novel avenue for genomic exploration and research. Single nucleotide variants (SNVs) inferred from next generation sequencing are expected to reveal gene mutations in cancer. However, NGS has lower sequence coverage and poor SNVs detection capability in the regulatory regions of the genome. Post probabilistic based methods are efficient for detection of SNVs in high coverage regions or sequencing data with high depth. However, for data with low sequencing depth, the efficiency of such algorithms remains poor and needs to be improved.

Results

A new tool SNVHMM basing on a discrete hidden Markov model (HMM) was developed to infer the genotype for each position on the genome. We incorporated the mapping quality of each read and the corresponding base quality on the reads into the emission probability of HMM. The context information of the whole observation as well as its confidence were completely utilized to infer the genotype for each position on the genome in study. Therefore, more probability power can be gained over the Bayes based methods, which is very useful for SNVs detection for data with low sequencing depth. Moreover, our model was verified by testing against two sets of lobular breast tumor and Myelodysplastic Syndromes (MDS) data each. Comparing against a recently published SNVs calling algorithm SNVMix2, our model improved the performance of SNVMix2 largely when the sequencing depth is low and also outperformed SNVMix2 when SNVMix2 is well trained by large datasets.

Conclusions

SNVHMM can detect SNVs from NGS cancer data efficiently even if the sequence depth is very low. The training data size can be very small for SNVHMM to work. SNVHMM incorporated the base quality and mapping quality of all observed bases and reads, and also provides the option for users to choose the confidence of the observation for SNVs prediction.

Background

In recent years, the advent of NGS technology has largely propelled the genomic research. NGS can generate millions of reads ranging from 30-350 base pairs (bp) based on the sequencing platform used. Continuous improvement in NGS technology brings the increasing of the throughput to a high extent and also lowers the cost [1]. With abundant reads aligned, many novel inferences can be made including regulatory element identification, mutation detection, gene expression estimation and detection of RNA splicing and fusion transcripts. NGS is expected to be a powerful tool for revealing genetic variations contributing to various complex diseases by providing sequence of a set of candidate genes, the whole exome or the whole genome. For example, whole genome sequencing can help in finding the frequency of tumor-specific point mutations for diseases such as multiple myeloma [2], while whole exome sequencing can be used to discover protein-coding mutation as well as small non-coding RNAs and aberrant transcriptional regulation that may contribute to diseases such as MDS [3].

The SNV calling algorithms can be divided into two categories. The first category includes threshold based commercial software packages such as Roche GSMapper and Lasergene, and the second category entails posterior probability based method including Maq [4], SOAPsnp [5], Varscan [6], Atlas-SNP2 [7] etc. For the threshold based prediction methods, a good threshold setting is difficult to obtain and relies heavily on the user experience [8].

In transcriptome based data, the number of reads representing a given transcript is highly variable across all genes making it difficult to determine a minimum depth. Moreover, the confidence for the prediction of each location is unavailable. Compared to the threshold based methods, posterior probability (Bayes) based methods achieve flexibility by considering the confidence of observation of each position on the genome. For the cancer genome sequencing data, sequencing errors, as well as the altered ploidy and tumor cellularity, are important factors affecting the accuracy of SNV calling. Although tools exist for SNVs discovery from NGS data, few are specifically suited to work with data from tumors. Recently, SNVMix [9] addressed this problem by incorporating the dependency of near-by genotypes and the posterior probability to improve the accuracy of SNVs prediction. However, the performance of SNVMix for data with low sequencing depth is not satisfactory compared to its performance with data having high sequencing depth. It has been observed that NGS provides lower sequence coverage in certain areas of genome including regulatory regions [10]. It is necessary to improve the performance of SNVs detection for tumor data with low sequencing depth. Moreover, SNVMix has achieved a relatively high sensitivity in the Bayesian framework, but the specificity is some low. The performance of specificity is needed to be improved further.

Hidden Markov model (HMM) is widely used in many fields such as speech and handwriting recognition, text classification, as well as DNA and protein classification [11]. Recently, a HMM based program VARiD [12] was developed for SNVs prediction for data from multiple sequencing platforms. VARiD is mainly focused on color space sequence and does not fully consider the mapping and base quality of the aligned reads and corresponding bases on the aligned reads in the considered model. Moreover, this method is time consuming for whole genome analysis and has not been used on RNA-Seq or whole exome sequence analysis from tumor data so far.

In this paper we developed an algorithm SNVHMM, for SNVs prediction of tumor data obtained from NGS basing on a discrete HMM. Since non-SNVs are prevalent and continuous in the genome [13], point mutations in cancer data are relevant to certain genes and are concentrated in the corresponding area [14, 15], the contextual information, especially for the non-SNVs, can be considered and made full use of in addition to the information from the overall distribution of traditional Bayesian framework. So SNVHMM is expected to gain more probability power from the contextual information on the genome compared to traditional Bayesian framework, and obtain better performance for SNVs prediction. Moreover, with the contextual information added to the whole distribution information, SNVHMM is also expected to improve the statistical performance of Bayesian method for tumor data with low sequencing depth.

Implementation

Problem formulation and SNVHMM model specification

We denote the length of the considered genome as L. Given the aligned reads for the sequence in study, we can get the depth L _t of the stated position t on the genome. The quality of the reads covering position t and the quality of corresponding bases on the reads are denoted as {\left\{r_i^t\right\}}_{i=1}^{L_t} and {\left\{q_i^t\right\}}_{i=1}^{L_t} respectively (Figure 1). We consider three genotypes for each stated position as {aa, ab, bb}, where {aa} denotes homozygous for the reference allele, {ab} denote heterozygous and {bb} denote homozygous for the non-reference allele. Our aim is to predict the genotype for each position on the genome, given the aligned reads.

Figure 1

Observation of HMM: the illustration for alignment at position t for the sequence in study. The observation sets {\left\{r_i^t\right\}}_{i=1}^{L_t} and {\left\{q_i^t\right\}}_{i=1}^{L_t} are considered as the observation of o _t  in HMM. For the bases on the covered reads, blue color denotes the base is the same as reference allele, yellow color denotes the base is different from reference allele while purple color denotes the base is undecided on the genome in study.

We denote the number of the hidden states as I. The hidden state and observation for each position are noted as S = {s _t }(t = 1, 2, ⋯, L) ∈ {v _i }(i = 1, 2, ⋯, I) and O = {o _t }(t = 1, 2, ⋯, L) respectively, where {\left\{v_i\right\}}_{i=1}^I are all states considered. The underlying genotypes of the sequenced genome are taken as the hidden states, which are interpreted as follows: (1) homozygous for normal; (2) heterozygous; (3) homozygous for mutation (Figure 2). These states are important in detecting single nucleotide polymorphism or point mutation for normal sample as well as cancer sample. The last two states are taken as SNV in our study. For simplicity, we note state {aa} as state 1, state {ab} as state 2 and state {bb} as state 3 in the following initial state distribution and state transition matrix.

• Initial state distribution:

\mathit{\pi }=\left\{{\pi }_1,{\pi }_2,{\pi }_3\right\},{\pi }_i=\mathrm{P}\left(s_1=v_i|t=1\right)

• State transition matrix:

\mathit{A}=\left(\begin{array}{ccc}\hfill a_{11}\hfill & \hfill a_{12}\hfill & \hfill a_{13}\hfill \\ \hfill a_{21}\hfill & \hfill a_{22}\hfill & \hfill a_{23}\hfill \\ \hfill a_{31}\hfill & \hfill a_{32}\hfill & \hfill a_{33}\hfill \end{array}\right)\triangleq \left({\mathit{A}}_1^T,{\mathit{A}}_2^T,{\mathit{A}}_3^T\right),

a_{\mathit{ij}}=\mathrm{P}\left(s_{t+1}=v_j|s_t=v_i\right)

• Emission probability distribution:

\mathit{B}=\left\{b_{s_t}\left({\mathit{o}}_{\mathit{t}}\right)\right\}

Figure 2

States of HMM: the illustration for three states {aa}, {ab} and {bb} in HMM. The meaning of different colors are defined the same as in Figure 1.

The observation to be considered for each position includes the coverage, the mapping quality of the covering reads and the base quality on the covering reads corresponding to each stated position. The observation for each position t is taken as {\mathit{o}}_{\mathit{t}}={\left\{q_i^t,r_i^t\right\}}_{i=1}^{L_t}. The emission probability b_{v_i}\left({\mathit{o}}_{\mathit{t}}\right) is calculated as a conditional probability, given the hidden state:

b_{s_t}\left({\mathit{o}}_{\mathit{t}}\right)=\mathrm{P}\left({\left\{q_i^t,r_i^t\right\}}_{i=1}^{L_t}|s_t=v_i\right)=f\left({\left\{q_i^t,r_i^t\right\}}_{i=1}^{L_t}\right)

(1)

To make full use of the mapping quality and base quality for each position on the genome, we compute b_{s_t}\left({\mathit{o}}_{\mathit{t}}\right) using the whole probability formula by considering if the covered reads are correctly aligned and if the corresponding bases on these reads are correctly called. We use a formula motivated by (5) in [9] by introducing a generalized Binomial distribution in addition to the conditional computation of the base calling probability and aligning probability.

b_{v_i}\left(o_t\right)=\left(\begin{array}{l}{L}_t\\ P_t\end{array}\right){\displaystyle {\prod }_{j=1}^{L_t}\left\{0.25\left(1-r_j^t\right)+0.5r_j^t\left[q_j^t{u}_i+\left(1-q_j^t\right)\left(1-u_i\right)\right]\right\}}

(2)

where {\left\{u_i\right\}}_{i=1}^I is the Binomial distribution parameter for each position on the genome and P _t is the number of reads having the same base with reference allele at position t. The detailed derivation of (2) is given at the supplementary file. In this study, we only considered two types of nucleotides covering the stated position, which have the largest and second largest number at the stated position. In the case of rare third alleles, these reads are assumed to be errors. In this study, u _i denotes the probability of occurrence for the allele having the largest number at the stated position.

Prior distribution of HMM

We take the initial distribution of π as Dirichlet distribution with hyper-parameter δ = (δ₁, δ₂, δ₃), u = (u₁, u₂, u₃) is taken conjugately according to a Beta distribution with hyper-parameter α = (α₁, α₂, α₃) and β = (β₁, β₂, β₃) as follows:

\mathrm{P}\left(\mathit{\pi }|\mathit{\delta }\right)=\mathit{\text{Dirichlet}}\left(\mathit{\pi }|\mathit{\delta }\right)

(3)

\mathrm{P}\left(u_k|{\alpha }_k,{\beta }_k\right)=\mathit{\text{Beta}}\left(u_k|{\alpha }_k,{\beta }_k\right)

(4)

where we take δ = (1000, 100, 100) by assuming that most positions will be homozygous for the reference allele. We also set α = (1000, 500, 1) and β = (1, 500, 1000) by assuming the probability of state {aa} occurring at the stated position is much larger than that it not occurring, vice versa for state {bb}. We also assume the probability of state {ab} occurring at the stated position is the same as that it not occurring. For the initial distribution of state transition matrix, we take the initial distribution of A _i as follows:

\mathbf{P}\left({\mathit{A}}_{\mathit{i}}|{\mathit{\gamma }}_{\mathit{i}}\right)=\mathit{\text{Dirichlet}}\left({\mathit{A}}_{\mathit{i}}|{\mathit{\gamma }}_{\mathit{i}}\right)

(5)

where we take γ ₁  = (1000, 100, 100), γ ₂  = (100, 1000, 100) and γ ₃  = (100, 100, 1000). Since the sum of elements in A _i should be equal to probability 1, a normalization for {\left\{{\mathit{A}}_{\mathit{i}}\right\}}_{i=1}^I is performed after each iteration of SNVHMM.

Estimation of HMM parameters

For simplicity, we denote the model parameters of HMM as λ ≜ (π, u, A) and learn the unknown HMM by using EM algorithm and computing the maximum likelihood estimation when the observed data are incomplete [16]. The aim is to find the model parameter λ maximizing the observation probability i.e. L(o, λ) ≜ P(o|λ) or log P(o|λ),where the later one is usually used when the length of the observation is large. We use a special case of EM algorithm, Baum-Welch algorithm [11], to learn the unknown parameters. For the training of HMM, we use the following auxiliary function Q\left(\mathit{\lambda },\overline{\mathit{\lambda }}\right) as the objective function for the optimization of the HMM parameters.

Q\left(\mathit{\lambda },\overline{\mathit{\lambda }}\right)={\Sigma }_{\mathit{s}\in \mathit{S}}\text{log}\mathrm{P}\left(\mathit{O},\mathit{S}|\overline{\mathit{\lambda }}\right)\mathrm{P}\left(\mathit{S}|\mathit{O},\mathit{\lambda }\right)

(6)

It is proved that maximizing the following auxiliary function can lead to the increase of the likelihood P(O|λ), i.e. \frac{\text{max}}{\overline{\mathit{\lambda }}}Q\left(\mathit{\lambda },\overline{\mathit{\lambda }}\right)\to \mathrm{P}\left(\mathit{O}|\overline{\mathit{\lambda }}\right)>\mathrm{P}\left(\mathit{O}|\mathit{\lambda }\right)[11]. Given model parameter set λ, P(O, S|λ) can be calculated as:

\mathrm{P}\left(\mathit{O},\mathit{S}|\mathit{\lambda }\right)={\pi }_{s_0}{\prod }_{t=1}^L{a}_{s_{t-1}{s}_t}{b}_{s_t}\left({\mathit{o}}_{\mathit{t}}\right)

(7)

Replacing term P(O, S|λ) in (6) with (7), (6) can be rewritten as:

\begin{array}{ll}Q\left(\mathit{\lambda },\overline{\mathit{\lambda }}\right)=& {\displaystyle {\sum }_{\mathit{s}\in \mathit{S}}\text{log}{\overline{\mathit{\pi }}}_{s_0}\mathrm{P}\left(\mathit{O},\mathit{S}|\mathit{\lambda }\right)}+{\displaystyle {\sum }_{\mathit{s}\in \mathit{S}}\left({\displaystyle {\sum }_{t=1}^L\text{log}{\overline{a}}_{s_{t-1}{s}_t}}\right)}\mathrm{P}(\mathit{O},\mathit{S}|\mathit{\lambda })\\ +{\displaystyle {\sum }_{\mathit{s}\in \mathit{S}}\left({\displaystyle {\sum }_{t=1}^L\text{log}{\overline{b}}_{s_t}\left({\mathit{o}}_{\mathit{t}}\right)}\right)}\mathrm{P}\left(\mathit{O},\mathit{S}|\mathit{\lambda }\right)\end{array}

(8)

The update of model parameters π _i and a _ij with constraints {\displaystyle {\sum }_{i=1}^N{\pi }_i}=1 and {\displaystyle {\sum }_{j=1}^N{a}_{\mathit{ij}}}=1 can be obtained by maximizing the first and second term of (8) with respect to π _i and a _ij respectively as follows:

{\overline{\pi }}_i=\frac{\mathrm{P}\left(\mathit{O},s_0=i|\mathit{\lambda }\right)}{\mathrm{P}\left(\mathit{O}|\mathit{\lambda }\right)}

(9)

{\overline{a}}_{\mathit{ij}}=\frac{{\displaystyle {\sum }_{t=1}^T\mathrm{P}\left(\mathit{O},s_{t-1}=i,s_t=j|\mathit{\lambda }\right)}}{{\displaystyle {\sum }_{t=1}^T\mathrm{P}\left(\mathit{O},s_{t-1}=i|\mathit{\lambda }\right)}}

(10)

The update for {\left\{u_i\right\}}_{i=1}^I can be obtained by maximizing the third term with respect to {\left\{u_i\right\}}_{i=1}^I, however, the close-form expression is not available due to the complicated structure of the observation term f\left({\left\{q_i^t,r_i^t\right\}}_{i=1}^{L_t}\right). We use a Newton iteration with respect to the first and second derivation of the third term as follows:

u_i^{\mathit{\text{new}}}=u_i^{\mathit{\text{old}}}-\frac{{\displaystyle {\sum }_{t=1}^T\frac{\partial f}{\partial u_i}*\frac{\mathrm{P}\left(\mathit{O},s_t=i|\mathit{\lambda }\right)}{f}}}{\frac{\partial f}{\partial u_i}\left\{{\displaystyle {\sum }_{t=1}^T\frac{\partial f}{\partial u_i}*\frac{\mathrm{P}\left(\mathit{O},s_t=i|\mathit{\lambda }\right)}{f}}\right\}}

(11)

The forward and backward algorithm [11] is used to update π _i and a _ij . In the implementation of Baum-Welch, the update of u _i is used for the update of the emission probability. Finally, we use Viterbi algorithm [11] to infer the hidden states of the sequence in study.

Results

Dataset

Two types of tumor data are used to verify the effectiveness of our model. The first type is the lobular breast tumor data with two different sequencing depths, which includes 497 positions generated using the Illumina GA II platform and was validated by Sanger. These positions were sequenced using Sanger capillary-based technology and were predicted to be non-synonymous protein-coding. 305 of these positions were confirmed as SNV and are taken as positive (TP), while 192 were not confirmed and are taken as true negative (TN). We take these positions as ground truth for the computation of TP, false positive (FP), TN, false negative (FN). These data can be obtained from the supplementary dataset 2A and 2C [9] along with their corresponding ground truth for SNVs in supplementary dataset 2B and 2D. The depths of supplementary dataset 2A and 2C are 10X and 40X respectively. Moreover, we use these datasets to compare between SNVHMM and SNVMix2, which is more efficient than SNVMix1 [9]. For better training of SNVMix2, we also use the supplementary dataset 3A and 3C [9].

The second dataset that came from MDS tumor data comprises of 7 MDS samples including 5 samples from RNA-Seq having depth <20 and 2 samples from whole exome sequencing having depth >150. These data are all from our lab and 4 mutated MDS genes were validated by PCR, along with other 23 common mutated MDS genes [14, 15] were also checked by SNVHMM for point mutation detection on these genes.

Statistical metrics

We take states of ‘ab’ and ‘bb’ as SNV for each location on the genome. Accuracy, sensitivity, specificity and F-score are proposed to evaluate the performance of SNVHMM and SNVMix2, which are defined as follows:

\mathit{\text{Accuracy}}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{TN}+\mathrm{FN}},\mathit{\text{Sensitivity}}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}

\mathit{\text{Specificity}}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}},F=\frac{2\left(\mathrm{Pr}\mathit{\text{ecision}}\times \mathrm{Re}\mathit{\text{call}}\right)}{\mathrm{Pr}\mathit{\text{ecision}}+\mathrm{Re}\mathit{\text{call}}}

where \mathit{\text{Precision}}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}} and \mathit{\text{Recall}}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}, are the proportion of true SNVs being predicted among the total predicted positives and the total true positives respectively.

Statistical performance

We compare the sensitivity, specificity, precision and F-score between SNVHMM and SNVMix2 on one lobular breast tumor data with sequencing depth 10X and 40X respectively. To get better classification results for SNVMix2, we used the supplementary datasets 3A and 3C [9] for training (SNVMix2_TO) which includes 14649 positions on different chromosomes. To test the gain of ability for SNVHMM in utilizing contextual information in addition to the information of whole distribution in SNVMix2, we train SNVHMM only by the test dataset itself, which is much smaller than the datasets 3A and 3C. For comparison purpose, we also tested SNVMix2 by training it on datasets 2A and 2C for 10X and 40X lobular breast tumor data respectively (SNVMix2_TI). The initial parameter for SNVHMM is set up as described in Section Prior distribution of HMM. Although we have specified how to encode base and mapping uncertainty into the emission probability in SNVHMM, obviating the need for taking thresholds for these quantities. However, different threshold setting of mapping quality (MQ) and base quality (BQ) as well as minimum and valid coverage (MVC) for each location enable us to achieve better performance. Here, MVC is the number of the least reads used to support the prediction of SNV. The results of SNVHMM and SNVMix2 for 10X and 40X lobular breast tumor data under different MQ and BQ condition are reported in Tables 1 and 2 respectively. The results are an average estimator basing on 20 independent runs for SNVHMM and SNVMix2 under different threshold settings. For each MQ and BQ condition, we choose the MVC achieving the best precision result for SNVHMM. The precision and F-score of SNVHMM with respect to different MVC for 10X and 40X lobular breast tumor data are plotted in Figure 3a, b and Figure 4a, b respectively. It is observed from Table 1 that SNVHMM performs significantly better than SNVMix2_TI and gains about 10% for precision and 3% for F-score in average. SNVHMM also outperforms SNVMix2_TO for nearly all the MQ and BQ conditions except MQ30_BQ10, although the gain for F-score is not obvious. It is also observed that SNVHMM decreases the false positive rate and increases the true negative rate compared with SNVMix2_TI and SNVMix2_TO, which leads to the improvement of specificity with a gain of 10% in average over SNVMix2_TO and 51% in average over SNVMix2_TI while maintaining a relatively high sensitivity. The same trend can also be observed in Table 2. SNVHMM outperforms both SNVMix2_TO and SNVMix2_TI for nearly all threshold conditions. Comparing SNVHMM with SNVMix2_TI, SNVHMM gains about 3.5% for precision and 14% for specificity. Comparing SNVHMM with SNVMix2_TO, SNVHMM gains about 2.1% for precision and 8.7% for specificity while maintaining a relatively high sensitivity. It is also noted that the performance of SNVHMM attains its peak and also achieves the best gain over SNVMix2_TI and SNVMix2_TO at MQ30&BQ10 condition.

Table 1 Comparison of statistical performance of SNVHMM with SNVMix2 with different mapping quality (MQ) and base quality (BQ) threshold for 10X data

Table 2 Comparison of statistical performance of SNVHMM with SNVMix2 with different mapping quality (MQ) and base quality (BQ) threshold for 40X data

Figure 3

Plot of the (a) accuracy and (b) F-score with the change of minimum and valid value of coverage for 10X lobular breast cancer data.

Figure 4

Plot of the (a) accuracy and (b) F-score with the change of minimum valid value of coverage for 40X lobular breast cancer data.

Comparing Table 1 with Table 2, we can observe the overall improvement for performance of SNVHMM, SNVMix2_TI and SNVMix2_TO. It is not surprising for more information can be used for SNVs prediction with the increase of sequencing depth. It is also noticed that the performance gap between SNVMix2_TI and SNVMix2_TO narrowed with the increase of sequencing depth while the gain of accuracy for SNVHMM over SNVMix2_TI is obvious. From Figures 3 and 4, it is observed that the best MVC value for precision and F-score is between 2-5 for 10X data and between 6-9 for 40X data respectively at different threshold conditions. So a moderate MVC ranging from 15% to 25% in proportion to the sequencing depth is needed for efficient prediction under different threshold settings. It is not surprising as a small MVC value can result in estimating some locations with low confidence and a large MVC value can cause some locations with high confidence being excluded from prediction. The performance of SNVHMM with respect to different MQ and BQ threshold and MVC condition is presented in Additional file 1: Table S1 and Table S2 for 10X and 40X data respectively. We also compute the p-value of accuracy and F-score for SNVHMM against SNVMix2_TI and SNVMix2_TO for different MQ and BQ threshold settings. For SNVMix2_TI and SNVHMM with sequencing depth 10X, the p-values for accuracy and F-score are less than 1.8e-23 and 3.2e-4 respectively (ANOVA, t-test), while the p-values for accuracy and F-score are less than 1.5e-12 and 9.9e-3 respectively (ANOVA, t-test) for SNVMix2_TI and SNVHMM with sequencing depth 40X. For SNVMix2_TO and SNVHMM with sequencing depth 10X, the p-values for accuracy and F-score are less than 1.3e-23 and 8.8e-12 respectively (ANOVA, t-test), while the p-values for accuracy and F-score are less than 1.9e-7 and 4.1e-5 respectively (ANOVA, t-test) for SNVMix2_TI and SNVHMM with sequencing depth 40X. So the gain of SNVHMM over SNVMix2_TI is statistical significant for both 10X and 40X data while SNVHMM is also significantly better than SNVMix2_TO, although the absolute gain of SNVHMM over SNVMix2_TO is limited. The ROC of SNVHMM and SNVMix2_TO for 20 independent runs over different MQ and BQ threshold are plotted in Figure 5, which show an obvious improvement for SNVHMM against SNVMix2_TO. Since the performance of SNVMix2_TO is better than SNVMix2_TI, the plot of ROC for SNVMix2_TI is not presented here.

Figure 5

ROC curve of SNVHMM and SNVMix2_TO for lobular breast cancer data.

Implementation and robust analysis

The proposed algorithm is implemented in C and supports both Maq [4] and SAMtools [17] pileup format. Running SNVHMM on the lobular breast cancer data with sequencing depth 40X takes 1 ~ 2 seconds and it needs ~20 seconds for the lobular breast cancer data including 14649 locations with sequencing depth 40X on 64 bit Linux Ubuntu 3.0.0. SNVHMM is robust under different MQ and BQ threshold settings. The standard deviations of accuracy and F-score are between 0.001 and 0.003 respectively for both 10X and 40X lobular breast cancer data.

The software is available online at https://sites.google.com/site/snvhmm4/. The initial setting and trained parameters are also available for the lobular breast cancer data and MDS data.

Performance of SNVHMM on MDS sample

To test the effectiveness of SNVHMM on large-size tumor data, we use SNVHMM on two groups of MDS samples to explore some common mutated genes for MDS. The two groups of data include 5 RNA-Seq samples and 2 whole exome samples. We take MQ = 50, BQ = 20 and MVC = 4 for RNA_Seq data and take MQ = 50, BQ = 20 and MVC = 7 for the whole exome data. We use ANNOVAR [18] for the annotation of point mutations. The number of point mutations detected by SNVHMM is reported in Table 3 and the corresponding annotated genes are reported in Additional file 1: Table S3. 18 common MDS mutated genes from [14] and 5 common MDS mutated genes from [15] are checked. Moreover, 4 MDS genes are validated by our lab. The validated information is presented in Additional file 1: Table S4. It can be seen that majority of the 27 genes in Additional file 1: Table S3 were detected in either RNA-Seq data or whole exome data. Some MDS mutated genes are detected in only few samples or not detected such as IDH1, IDH2 and PTPN11, which are proved to be rare in MDS [19]. Finally, 4 new mutated MDS genes: MLL3, IQGAP2, DIDO1 and EIF4G2 are all detected by SNVHMM to be non-synonymous, which coincides with our validated result.

Table 3 Number of point mutations found in 5 MDS RNA-Seq data and 2 whole exome data for SNVHMM

Discussion

We introduced a new algorithm SNVHMM for SNVs prediction of tumor data from next generation sequencing, which generally yield data with low sequencing depth due to sequencing errors, as well as the altered ploidy and tumor cellularity. SNVHMM was conceived to circumvent the shortcomings of existing algorithms that cannot efficiently predict SNVs for data with low sequencing coverage. In this algorithm, we considered three genotypes concerned as the hidden states of HMM, and incorporated the confidence of the observation into the emission probability in HMM. The performance of SNVHMM was compared with a recently published method SNVMix2. Compared to SNVMix2, SNVHMM considered the relation of state from near-by locations in addition to their distribution. Moreover, SNVHMM predicted the hidden states by maximizing the posterior probability in condition of the whole observation while SNVMix2 predicted the genotype basing on maximizing the posterior probability in condition of the observation from single location. So SNVHMM gained more probability power for prediction from the same dataset. It was shown by experiment from the lobular breast cancer data sequenced with lower depth that SNVHMM improved the performance of SNVMix2 by only using much smaller size of training data. It was also observed that SNVHMM even exceeded the performance of SNVMix2 trained by much larger datasets. If looking into the performance of SNVHMM and SNVMix2_TI for lobular cancer data, we found that SNVHMM corrected 42% ~ 75% and 26% ~ 33% false positives to true negatives for sequencing depth 10X and 40X respectively, with more than 85% of them to have coverage less than 20. For SNVHMM and SNVMix2_TO, SNVHMM also corrected 10% ~ 17% and 2% ~ 8% false positives to true negatives for sequencing depth 10X and 40X respectively, with more than 80% of them to have coverage less than 20. So SNVHMM improved the performance of SNVMix2 for low-coverage sequencing data or at the low depth area of genome by improving the true negative rate largely. This verified the effectiveness of SNVHMM in utilizing the contextual information of non-SNVs by improving the specificity largely while remaining a relatively high sensitivity.

For experiments on MDS samples, SNVHMM could detect the point mutations efficiently. More than 95% of the point mutations detected by both SNVHMM and SNVMix2 are obvious mutation as most of the covered reads have the same non-reference base. From the common mutated 27 genes list of MDS, most of the mutated genes can be found in majority of the samples by SNVHMM. We also examined the region of the genes not detected to find no non-reference bases covered or insufficient non-reference bases covered with low quality.

For the training of the HMM parameters, we give the parameters before and after training for the lobular breast cancer data and MDS data in Table 4. The threshold of mapping quality and base quality are taken as 50 and 20 respectively and the initial setting of the parameters are taken from the distribution as defined in the “Prior distribution of HMM” Section. For the lobular breast cancer data, we both considered the 10X and 40X conditions. It is observed that majority of the parameters changed significantly while only a few did not change much. The trained π differs largely from the initial value for the lobular breast cancer data. It is not surprising as π indicates the proportion of three kinds of bases on the genome, and it is observed from the ground truth of the lobular breast cancer data that majority of the bases have the state “ab”. It is also observed that majority parameters of the trained u and A changed significantly. The initial setting of u and A seems to be some close to the true parameters of these data. For the MDS data, the trained π and A changed largely compared with the initial setting while μ also changed but not as much as π and A. We check the raw data to find that the SNVs are sparse and state “aa” dominates the whole distribution, so it is reasonable that the first value in π increased largely while the second and the third value in π decreased. It is also not surprising that all the states have a large transition probability to state “aa”. All the trained parameters for different thresholds of lobular breast cancer data and MDS data are provided at https://sites.google.com/site/snvhmm4/.

Table 4 Comparison of the parameters of SNVHMM before and after training on lobular breast cancer (LBC) data and MDS RNA_seq data for threshold of mapping quality 50 and base quality 20

Conclusions

We have proposed a new SNVs prediction tool SNVHMM for cancer data from NGS. SNVHMM can gain more probability power from the transition probability in additional to the posterior probability computation for the genotype distribution of whole observation. So SNVHMM is very efficient when the depth of NGS data is very low. Since NGS has lower sequence coverage and poor SNV detection capability in the regulatory regions of the genome, it is very helpful for SNV prediction for the low-depth area on the genome. SNVHMM outperformed an existing SNV prediction tool SNVMix by reducing its false positives and increasing its true negative. Moreover, SNVHMM needs much less data for training while obtaining a better performance than SNVMix. Finally, two types of MDS data with different coverage are tested, which shows the effectiveness of SNVHMM.

Availability and requirements

 Project name: SNVHMM: predicting single nucleo tide variants from next generation sequencing

Project home page: https://sites.google.com/site/snvhmm4/

Operating system: 64-bit Linux

Programming language: C

Other requirements: Linux Ubuntu 3.0.0 or higher

License: GNU GPL

Any restrictions to use by non-academics: license  needed.