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 tumorspecific point mutations for diseases such as multiple myeloma [2], while whole exome sequencing can be used to discover proteincoding mutation as well as small noncoding 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], AtlasSNP2 [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 nearby 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 RNASeq 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 nonSNVs 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 nonSNVs, 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 and 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 nonreference 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 and 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 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:
• State transition matrix:
• Emission probability distribution:
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 . The emission probability is calculated as a conditional probability, given the hidden state:
To make full use of the mapping quality and base quality for each position on the genome, we compute 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.
where 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 hyperparameter δ = (δ_{1}, δ_{2}, δ_{3}), u = (u_{1}, u_{2}, u_{3}) is taken conjugately according to a Beta distribution with hyperparameter α = (α_{1}, α_{2}, α_{3}) and β = (β_{1}, β_{2}, β_{3}) as follows:
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:
where we take γ_{1} = (1000, 100, 100), γ_{2} = (100, 1000, 100) and γ_{3} = (100, 100, 1000). Since the sum of elements in A_{i} should be equal to probability 1, a normalization for 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, BaumWelch algorithm [11], to learn the unknown parameters. For the training of HMM, we use the following auxiliary function as the objective function for the optimization of the HMM parameters.
It is proved that maximizing the following auxiliary function can lead to the increase of the likelihood P(Oλ), i.e. [11]. Given model parameter set λ, P(O, Sλ) can be calculated as:
Replacing term P(O, Sλ) in (6) with (7), (6) can be rewritten as:
The update of model parameters π_{i} and a_{ij} with constraints and can be obtained by maximizing the first and second term of (8) with respect to π_{i} and a_{ij} respectively as follows:
The update for can be obtained by maximizing the third term with respect to , however, the closeform expression is not available due to the complicated structure of the observation term . We use a Newton iteration with respect to the first and second derivation of the third term as follows:
The forward and backward algorithm [11] is used to update π_{i} and a_{ij}. In the implementation of BaumWelch, 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 capillarybased technology and were predicted to be nonsynonymous proteincoding. 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 RNASeq 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 Fscore are proposed to evaluate the performance of SNVHMM and SNVMix2, which are defined as follows:
where and , 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 Fscore 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 Fscore 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 Fscore in average. SNVHMM also outperforms SNVMix2_TO for nearly all the MQ and BQ conditions except MQ30_BQ10, although the gain for Fscore 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) Fscore 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) Fscore 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 Fscore 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 pvalue of accuracy and Fscore 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 pvalues for accuracy and Fscore are less than 1.8e23 and 3.2e4 respectively (ANOVA, ttest), while the pvalues for accuracy and Fscore are less than 1.5e12 and 9.9e3 respectively (ANOVA, ttest) for SNVMix2_TI and SNVHMM with sequencing depth 40X. For SNVMix2_TO and SNVHMM with sequencing depth 10X, the pvalues for accuracy and Fscore are less than 1.3e23 and 8.8e12 respectively (ANOVA, ttest), while the pvalues for accuracy and Fscore are less than 1.9e7 and 4.1e5 respectively (ANOVA, ttest) 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.
Additional file 1: Table S1. Statistical performance of SNVHMM for different minimum and valid coverage (d), as well as for different MQ and BQ value when the sequencing depth of lobular breast cancer data is 10X. Table S2: Statistical performance of SNVHMM for different minimum and valid coverage (d), as well as for different MQ and BQ value when the sequencing depth of lobular breast cancer data is 40X. Table S3: 23 reported mutated genes in Bejar,R. et al. (2011) and Thol,F. et al. (2012) are checked by SNVHMM. 4 new genes that are found in 5 MDS RNASeq sample and 2 MDS whole exome samples are found by SNVHMM and validated by our lab. Table S4: description of 4 MDSrelated mutated genes found by SNVHMM and validated by our lab in 5 RNASeq and 2 whole exome samples.
Format: DOC Size: 332KB Download file
This file can be viewed with: Microsoft Word Viewer
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 Fscore 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/ webcite. 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 largesize 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 RNASeq 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 RNASeq 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 nonsynonymous, which coincides with our validated result.
Table 3. Number of point mutations found in 5 MDS RNASeq 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 nearby 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 lowcoverage 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 nonSNVs 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 nonreference 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 nonreference bases covered or insufficient nonreference 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/ webcite.
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 lowdepth 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: 64bit Linux
Programming language: C
Other requirements: Linux Ubuntu 3.0.0 or higher
License: GNU GPL
Any restrictions to use by nonacademics: license needed.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
JB and XZ performed the design of the study. JB, CL, JX and HW participated in the implementation of the algorithm. CL and JB prepared and analyzed the breast and MDS data. JB drafted the manuscript. PK revised the manuscript. All authors read and approved the final manuscript.
Acknowledgements
We would like to acknowledge the members of Translational Biosystems Lab in Cornell Medical School and Dr. Jing Su for his help with programming.
Funding
This work was supported by Funding: NIH R01LM01018503 (Zhou), NIH U01HL11156001 (Zhou), NIH 1R01DE02267601 (Zhou), U01 CA16688601 (Zhou).
References

Shendure J, Ji H: Nextgeneration DNA sequencing.
Nat Biotechnol 2008, 26:11351145. PubMed Abstract  Publisher Full Text

Chapman MA, et al.: Initial genome sequencing and analysis of multiple myeloma.
Nature 2011, 471:467472. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Beck D, et al.: Integrative analysis of next generation sequencing for small noncoding RNAs and transcriptional regulation in myelodysplastic syndromes.
BMC Med Genomics 2011, 4:419. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Li H, et al.: Mapping short DNA sequencing reads and calling variants using mapping quality scores.
Genome Res 2008, 18:18511858. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Li RQ, et al.: SNP detection for massively parallel wholegenome resequencing.
Genome Res 2009, 19:11241132. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Koboldt DC, et al.: VarScan: variant detection in massively parallel sequencing of individual and pooled samples.
Bioinformatics 2009, 25:22832285. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Shen Y, et al.: A SNP discovery method to assess variant allele probability from nextgeneration resequencing data.
Genome Res 2010, 20:273280. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Martin ER, et al.: SeqEM: an adaptive genotypecalling approach for next generation sequencing studies.
Bioinformatics 2010, 26:28032810. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Goya R, et al.: SNVMix: predicting single nucleotide variants from next generation sequencing of tumors.
Bioinformatics 2010, 26:730736. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wang W, et al.: Next generation sequencing has lower sequence coverage and poorer SNPdetection capability in the regulatory regions.

Rabiner LR: A tutorial on hidden Markov models and selected applications in speech recognition.
Proc IEEE 1989, 77:257286. Publisher Full Text

Adrian VD, et al.: VARiD: a variation detection framework for colorspace and letterspace platforms.
Bioinformatics 2010, 26:i343i349. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

The International SNP Map Working Group: A map of human genome sequence variation containing 1.42 Million single nucleotide polymorphisms.
Nature 2001, 409:928933. PubMed Abstract  Publisher Full Text

Bejar R, et al.: Clinical effect of point mutations in myelodysplastic syndromes.
N Engl J Med 2011, 364:24962506. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Thol F, et al.: Frequency and prognostic impact of mutations in SRSF2, U2AF1, and ZRSR2 in patients with myelodysplastic syndromes.
Blood 2012, 119:35783584. PubMed Abstract  Publisher Full Text

Dempster AP, et al.: Maximum likelihood from incomplete data via the EM algorithm.

Li H, et al.: The sequence alignment/Map format and SAMtools.
Bioinformatics 2009, 25:20782079. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wang K, Li M, Hakonarson H: ANNOVAR: functional annotation of genetic variants from nextgeneration sequencing data.
Nucleic Acids Res 2010, 38:e164. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bejar R, et al.: Unraveling the molecular pathophysiology of myelodysplastic syndromes.