Confidence intervals for Markov chain transition probabilities based on next generation sequencing reads data

Abstract

Background: Markov chains (MC) have been widely used to model molecular sequences. The estimations of MC transition matrix and confidence intervals of the transition probabilities from long sequence data have been intensively studied in the past decades. In next generation sequencing (NGS), a large amount of short reads are generated. These short reads can overlap and some regions of the genome may not be sequenced resulting in a new type of data. Based on NGS data, the transition probabilities of MC can be estimated by moment estimators. However, the classical asymptotic distribution theory for MC transition probability estimators based on long sequences is no longer valid.

Methods: In this study, we present the asymptotic distributions of several statistics related to MC based on NGS data. We show that, after scaling by the effective coverage d defined in a previous study by the authors, these statistics based on NGS data approximate to the same distributions as the corresponding statistics for long sequences.

Results: We apply the asymptotic properties of these statistics for finding the theoretical confidence regions for MC transition probabilities based on NGS short reads data. We validate our theoretical confidence intervals using both simulated data and real data sets, and compare the results with those by the parametric bootstrap method.

Conclusions: We find that the asymptotic distributions of these statistics and the theoretical confidence intervals of transition probabilities based on NGS data given in this study are highly accurate, providing a powerful tool for NGS data analysis.

Graphical abstract

Keywords

Markov chains / next generation sequencing / transition probabilities / confidence intervals

INTRODUCTION

Markov chains (MC) have been widely used to model molecular sequences [¹,²]. They have been used to study the dependencies between the bases [³], the enrichment and depletion of certain word patterns [⁴], prediction of occurrences of long word patterns from short patterns [⁵,⁶], and the detection of signals in introns [⁷]. Narlikar et al. [⁸] studied the effect of the order of MCs on several biological problems including phylogenetic analysis, assignment of sequence fragments to different genomes in meta-genomic studies, motif discovery, and functional classification of promoters. Usually the transition probabilities of the Markov chain are estimated from data, and hence the uncertainty of the estimated transition probabilities has to be taken into account. This is done for example in [²] for word counts in long sequences.

Often DNA samples do not consist of a few long sequences, but of many short sequences, for example based on next generation sequencing (NGS) technology. For short sequences the standard asymptotic results do not apply, see for example [⁹]. In [⁹] it was shown that the word count standardization includes an additional scaling factor, the effective coverage d. While [⁹] considered word count statistics under a Markov chain with estimated transition probabilities, in this study we are interested in finding confidence regions for the transition probabilities. We derive a normal and a chi-square approximation for several statistics related to MC based on many short sequences (also called reads) from a long underlying sequence that is assumed to be a realization of a Markov chain; and again the effective coverage d appears.

Consider an r-th order Markov chain A = A₁A₂… A_n of length n. The maximum likelihood estimate of transition probability from a word w of length r to a letter b is given by

{\widehat{P}}_{\mathrm{w},b}=f_{\mathrm{w}b}/f_{\mathrm{w}},

where fw is the number of occurrences of word w within the sequence. Many investigators have studied the limit distribution of {\widehat{P}}_{w,b} and other related statistics when the sequence length n tends to infinity [¹⁰,¹¹]. In the classical papers [¹⁰,¹¹], Billingsley studied many important problems related to MC. In particular, he showed the following important results.

Suppose the long sequence follows a simple first order Markov chain, with finite state space \mathcal{A} of size s. We call the elements of \mathcal{A} letters, as we will be thinking of a Markov chain on the set of nucleotides or amino acids. Denote the transition probability from letter i to letter j by p_ij , and let

(1){\xi }_{ij}=(f_{ij}-f_i{p}_{ij})/f_i^{\frac{\mathrm{1}}{\mathrm{2}}}

Theorem 3.1 in [¹¹] shows that if the Markov chain is stationary and ergodic, then as the sequence length n → ∞ the distribution of the s²-dimensional random vector x = (x_i,j ) converges to the normal distribution with mean 0 and covariance matrix (l_i,j;k,l), where

(2){\lambda }_{i,j;k,l}={\delta }_{i,k}({\delta }_{j,l}{p}_{ij}-p_{ij}{p}_{kl}),

where du,v is the Kronecker delta that is 1 when u = v, and 0 otherwise. Let z = (z1,…, zs) be the random vector with components

(3){\zeta }_i=(f_i-n{p}_i)/n^{\frac{\mathrm{1}}{\mathrm{2}}},

Theorem 3.3 in [¹¹] also shows that the distribution of the random vector z converges to the normal distribution with covariance matrix {\alpha }_{ij}+O\left({\displaystyle \frac{\mathrm{1}}{n}}\right), where

{\alpha }_{ij}={\delta }_{ij}{p}_i-p_i{p}_j+p_i{\displaystyle \underset{m=\mathrm{1}}{\overset{\infty }{\mathrm{\Sigma }}}}(p_{ij}^{(m)}-p_j)+p_j{\displaystyle \underset{m=\mathrm{1}}{\overset{\infty }{\mathrm{\Sigma }}}}(p_{ji}^{(m)}-p_i).

In [¹¹] it is also shown that the statistic

(4)S={\displaystyle \underset{ij}{\mathrm{\Sigma }}}{\displaystyle \frac{{(f_{ij}-f_i{p}_{ij})}^{\mathrm{2}}}{f_i{p}_{ij}}}.

is asymptotically chi-square in distribution, and the number of degrees of freedom is s(s - 1) (assuming the transition probabilities p_{ijare all positive), where s is the size of the alphabet set \mathcal{A}.}

These results were derived under the assumption of one sample consisting of a long realisation of the Markov chain. With the development of next generation sequencing technologies, instead of one long sequence, typically short fragments from the long genome sequence are sampled, while assembly of these genome fragments as in [¹²,¹³] can be challenging.

In this study, we investigate the asymptotic distributions of similar statistics described above as in [¹⁰,¹¹], but under the new model proposed in [⁹,¹⁴] for the NGS short reads data (see Fig. 1 and Section of “Methods” for the details of the model). We demonstrate that, after scaling by the effective coverage d proposed in [⁹], these statistics converge to the same corresponding distributions as those statistics for long sequences. We further apply the asymptotic properties of these statistics for finding the theoretical confidence regions for MC transition probabilities based on NGS short reads data. We validate these asymptotic distributions and our theoretical confidence intervals using simulated as well as real data sets. In addition, we propose a parametric bootstrap method, and compare the bootstrap confidence intervals with our theoretical confidence intervals. We find that, the asymptotic properties of these statistics and the theoretical confidence intervals given in this study are highly accurate for the NGS short reads data, providing a powerful tool for NGS data analysis.

RESULTS

Our theoretical results are mainly presented in Section of “Methods”: we first introduce the probabilistic model for NGS data based on a MC sequence and random sampling of short reads; we then derive the statistics of interest based on NGS reads and the effective coverage d; we then develop two propositions on the asymptotic distributions of these statistics, and finally we propose the theoretical confidence intervals as well as the parametric bootstrap confidence intervals for MC transition probabilities based on NGS short reads data. Although we cannot rigorously prove the major propositions developed in this study, we validate statistical properties of these NGS related statistics and theoretical confidence intervals of transition probabilities using simulated NGS data sets as well as a real viral genome sequence data in this section.

Validating the theoretical results using simulated NGS data sets

We validate Proposition 1, Proposition 2 and the theoretical confidence intervals for MC transition probabilities of Eq. 16 in Section of “Methods” using simulated NGS short reads data sets.

For each simulated NGS data set, we first simulate the underlying genome sequence with MC using the same set of parameters of transition probabilities P = [p_{ij] as in [⁹], which are listed in Supplementary Table S1. The initial position of the simulated genome is sampled based on the stationary distribution of the MC. We set the genome length G = 100, 000 for all simulated data sets. We then randomly sample M reads of length b bps from the simulated underlying genome. The reads are sampled based on the Poisson process as described in Section of “Methods”. The Poisson process can be homogeneous or inhomogeneous.}

We simulate 3 NGS data sets with reads sampled by homogeneous Poisson processes. The parameters of the 3 homogeneous data sets are (1) M = 500, b = 100 (denoted as H1 thereafter); (2) M = 1000, b = 200 (denoted as H2 thereafter); and (3) M = 1000, b = 300 (denoted as H3 thereafter).

For each of these homogeneous data sets, the rate c of Poisson processes is calculated by

c={\displaystyle \frac{M}{G-\beta }}\approx {\displaystyle \frac{M}{G}}.

Noting Remark 1, the sequencing coverage C = bc = bM/G. We follow Eq. 13 and calculate the corresponding effective coverage d as

d=\mathrm{1}+\beta c=\mathrm{1}+{\displaystyle \frac{\beta M}{G}}.

For real NGS data, the reads are generally generated inhomogeneously from the genome. To test our results for inhomogeneous cases, we also simulate 3 NGS data sets with reads sampled by inhomogeneous Poisson processes. To achieve this, we divide the long underlying genome sequence into 100 consecutive non-overlapped bins b1, … , bt, … , b100 with the same size; the sampling rates along the positions within the same bin bt are equal; and the rates of the different bins are proportional to 100 random variables drawn from the gamma distribution G(1, 20) [¹⁵]. The simulated 3 inhomogeneous data sets of different read numbers and read lengths are (1) M = 500, b = 100 (denoted as IH1 thereafter); (2) M = 1000, b = 200 (denoted as IH2 thereafter); and (3) M = 1000, b = 300 (denoted as IH3 thereafter).

Validating the approximate normal distribution of {\xi }_{ij}^{(R)} in Proposition 1

To validate the approximate normal distribution of {\xi }_{ij}^{(R)} by Proposition 1, we first standardize the {\xi }_{ij}^{(R)} as

\widetilde{{\xi }_{ij}^{(R)}}={\displaystyle \frac{{\xi }_{ij}^{(R)}}{\sqrt{d{p}_{ij}(\mathrm{1}-p_{ij})}}}\to N(\mathrm{0},\mathrm{1}),

where the estimated values of d and p_i,j are used in the standardization of {\xi }_{ij}^{(R)}. For each of the 6 data sets (H1, H2, H3, IH1, IH2, and IH3), we repeat the simulating processes by the same parameters 2,000 times to generate 2,000 samples, and then calculate the \widetilde{{\xi }_{ij}^{(R)}} value for each sample, and check whether 2,000 values of \widetilde{{\xi }_{ij}^{(R)}} have the standard normal distribution by the Kolmogorov-Smirnov (KS) test.

For the homogeneous data sets, for all ij, the \widetilde{{\xi }_{ij}^{(R)}}'s have approximately mean of 0 and variance of 1, with p-values being uniformly distributed in the region of (0, 1) (Fig. 2), indicating that all \widetilde{{\xi }_{ij}^{(R)}}'s follow the standard normal distribution. We obtain similar results for the inhomogeneous data sets (Supplementary Fig. S1). These results show strong validations of Proposition 1.

Validating the approximate distributions of {\zeta }_i^{(R)} and S^{(R)} in Proposition 2

To validate the approximate distributions of {\zeta }_i^{(R)} and S^{(R)} in Proposition 2, we first standardize them as follows:

\widetilde{{\zeta }_i^{(R)}}={\displaystyle \frac{{\zeta }_i^{(R)}}{\sqrt{d{\alpha }_{ii}}}}\to N(\mathrm{0},\mathrm{1}),

and

\widetilde{S^{(R)}}={\displaystyle \frac{S^{(R)}}{d}}\to {\chi }^{\mathrm{2}}(s(s-\mathrm{1})).

We then check on the simulated data sets, whether the new statistic \widetilde{{\zeta }_i^{(R)}} has the standard normal distribution for all i ∈ A, and the new statistic \widetilde{{\zeta }_i^{(R)}} has the chi-square distribution with df = s(s - 1) = 4 × 3= 12, respectively.

Our simulation results show that: 1) for all i ∈ A, the \widetilde{{\zeta }_i^{(R)}} approximates to the standard normal distribution quite well for both homogeneous data sets (Supplementary Fig. S2) and inhomogeneous data sets (Supplementary Fig. S3); 2) the \widetilde{S^{(R)}} approximates a chi-square distribution with df = 12 (Mean=12 and Variance=24) on both homogeneous data sets and inhomogeneous data sets as well: the calculated values of (Mean, Variance and p-value) are (12.02, 25.15, 0.402) for H1, (12.00, 24.75, 0.690) for H2, (12.14, 24.82, 0.556) for H3, (11.81, 22.92, 0.058) for IH1, (12.15, 23.65, 0.064) for IH2, and y (11.88, 24.18, 0.453) for IH3.

Validating the theoretical confidence intervals of MC transition probabilities using NGS data

One of the major applications of Proposition 1 is to construct the theoretical confidence intervals of MC transition probabilities using NGS data. We validate the derived theoretical level 1 -a confidence interval Eq.16 by simulation study.

For each of the simulated homogeneous (H1, H2 and H3) and inhomogeneous (IH1, IH2 and IH3) data sets, we

1. Simulate one sample, estimate {\widehat{p}}_{ij}s,

2. Calculate theoretical level 1 -a (a ∈ {0.01, 0.05, 0.1}) confidence interval for p_{ij} by Eq.16;

3. Repeat Steps 1 and 2 for 1,000 times, and calculate the fraction of times that the theoretical level 1 -a confidence intervals cover the true p_{ij}.

The simulation results in Fig. 3 and Supplementary Fig. S4 show that, the fractions of times that the theoretical confidence intervals cover the true p_{ij} are very close to 1 -a (a = 0.01, 0.05, 0.1) for both homogeneous and inhomogeneous data sets.

Calculating theoretical confidence intervals of p_{ij} based the estimated \widehat{d}

In real study, the effective coverage d is generally not available for the NGS data sets. To deal with data sets with unknown d, we adopt the approach proposed in [⁹], and estimate the effective coverage d as follows:

(5)\widehat{d}=median\{{(Z_{\mathbf{w}}^{(R)})}^2,w\in {\mathcal{A}}^k\}/\mathrm{0.}\mathrm{456}

where Z_{\mathbf{w}}^{(R)} is given in Eq. 14.

We find that, by plug-in the estimated effective coverage \widehat{d} to the theoretical confidence intervals for p_{ij} in Eq.16, we can still achieve similar accuracies in calculating the theoretical confidence intervals (Supplementary Figs. S5, S6).

Validating the theoretical confidence interval for p_{ij} by parametric bootstrap

For each of the homogeneous NGS data sets (H1, H2, H3), we calculate the confidence intervals by parametric bootstrap method. We repeat 100 bootstrap calculation for each data set. We find that the fraction of times that the bootstrap confidence intervals covering the true transition probabilities approximate the level 1 -a well (Supplementary Figs. S7S9). Meanwhile, the theoretical confidence intervals by Eq.16 achieve similar results as the parametric bootstrap confidence intervals both in 1) the fraction of time cover the true transition probability (Supplementary Figs. S7S9), and 2) lengths of confidence intervals (Supplementary Fig. S10).

Application to real viral genome sequence data

We apply the theoretical results developed in this study on a real viral genome sequence Bacillus phage SP-beta. The length G of the genomic sequence, which was downloaded from NCBI, is 134, 416 in bp. The frequencies of A, G, C and T are 45,580 (33.91%), 21,078 (15.68%), 25,477 (18.95%) and 42,281 (31.46%), respectively. We model the Bacillus genome sequence with a 2ndorder MC after estimating MC order by method developed in [⁹]. We estimate the transition probability based on the genome sequence as

{\hat{p}}_{ij,k}={\displaystyle \frac{f_{ijk}}{f_{ij}}},

which are shown in Supplementary Table S2. The stationary distribution p of the 2nd-order MC is calculated based on the estimated transition probabilities (Supplementary Table S2). Based on Theorem 3.1 in Billingsley [¹¹], we derive the theoretical confidence for long genome sequence as follows

(6){\mathrm{CI}}_{\mathrm{1}-\alpha }^{\mathrm{Long}}=\left({\hat{p}}_{ij,k}-z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}\sqrt{{\displaystyle \frac{{\hat{p}}_{ij,k}(\mathrm{1}-{\hat{p}}_{ij})}{f_{ij}}}},\hspace{0.17em}{\hat{p}}_{ij,k}+z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}\sqrt{{\displaystyle \frac{{\hat{p}}_{ij,k}(\mathrm{1}-{\hat{p}}_{ij,k})}{f_{ij}}}}\right)

See Table 1 for an example of the confidence intervals for p_AA,A by (6).

We then simulate NGS data sets by generating short reads from the Bacillus phage SP-beta genome sequence based on both homogeneous and inhomogeneous Poisson processes. The read length is b = 200 in bp and we generate NGS data sets using coverage C = bc ∈ {0.5, 1, 2, 5, 10}. We then estimate the MC transition matrix using NGS reads as

{\tilde{p}}_{ij,k}={\displaystyle \frac{f_{ijk}^{(R)}}{f_{ij}^{(R)}}},

and extend the theoretical confidence intervals of Eq.16 for NGS to the 2nd-order MC transition probabilities p_ij,k as follows

(7){\mathrm{CI}}_{\mathrm{1}-\alpha }^{\mathrm{NGS}}=\left({\widehat{p}}_{ij,k}-z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}\sqrt{{\displaystyle \frac{d{\widehat{p}}_{ij,k}(\mathrm{1}-{\widehat{p}}_{ij,k})}{f_{ij}^{(R)}}}},\hspace{0.17em}{\widehat{p}}_{ij,k}+z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}\sqrt{{\displaystyle \frac{d{\widehat{p}}_{ij,k}(\mathrm{1}-{\widehat{p}}_{ij,k})}{f_{ij}^{(R)}}}}\right).

We calculate the theoretical confidence intervals by Eq. (7) as well as the bootstrap confidence intervals by Eq. (19) for transition probability p_AA,A on homogeneous NGS data sets (Table 1). Although we do not have true transition probabilities to validate these confidence intervals, we find that the theoretical confidence intervals by Eq. (7) as well as the parametric bootstrap by Eq. (19) are very close to the confidence intervals estimated from the long sequence by Eq. 6 (Table 1).

The theoretical confidence intervals by Eq. (7) for inhomogeneous NGS data sets are list in Supplementary Table S3.

DISCUSSION

In this study, we proposed statistics for estimations of MC transition probabilities and their corresponding confidence intervals based on next generation sequence reads. Although we did not rigorously prove the two propositions (Propositions 1 and 2) on the asymptotic distributions of {\xi }_{ij}^{(R)}, {\zeta }_i^{(R)} and S^(R), our simulation studies showed that the propositions and the derived formulas for the confidence intervals of the transition probabilities perform decently well on simulated NGS datasets with wide regions of sequencing coverage.

In addition, we proposed a parametric bootstrap method for the confidence intervals of MC transition probabilities based on NGS reads data, and showed that both theoretical confidence intervals derived from the two propositions and the bootstrap methods achieve similar results with high accuracies, further supporting the validity of Propositions 1 and 2. However, it is worth noting that the parametric bootstrap based approach is computationally expensive in sampling data sets many times; while the theoretical confidence intervals bypass the need of the time consuming sampling step, thereby making it highly applicable to large-scale NGS datasets.

The developed theory and methods implicitly assume error-free for the NGS short reads data. Since the sequencing error rates of bulk sequencing are lower than 1% for Illumina sequencing machines, the error-free assumption on NGS data will not be a big issue. The rapid advances in single-cell sequencing (SCS) [¹⁶,¹⁷] provide unprecedented insights of cell heterogeneity. The statistical theory and methods developed in this study will be promising for the quantitative characterizations of SCS data. However, due to limited molecular materials within a single cell as well as high dropout events, the SCS data are error-prone. For SCS data analysis, the sequencing errors are non-negligible and will blur the estimators, making their asymptotic distributions biased. We will study the asymptotic distribution theory of statistics based on the error-prone NGS data in our future work.

METHODS

Probabilistic modeling of a MC sequence and random sampling of the reads using NGS

In NGS, a large number of reads are randomly sampled from the genome. Hence two random processes are involved in the generation of the short reads data: i) the generation of the underlying genome sequence with genome length G base pairs (bps), and ii) random sampling of the M reads with length b from the simulated genome in i). See Fig. 1 for an illustration.

In this study, we assume the underlying genome sequence is generated by a Markov Chain. We treat mainly the case of a first-order Markov chain, as the results can be easily adapted to higher order Markov chains. We thus use a first-order stationary and ergodic MC to model the underlying genome sequence with each letter taking values in a finite alphabet set \mathcal{A} of size s. Since our study is based on genomic sequences, \mathcal{A} = {A, C, G, T} and s = 4. We denote the transition matrix of MC as P = [p_ij]4×4.

We model the random sampling of the reads using a Poisson process as in the Lander-Waterman model [¹⁸]. The Poisson process \rho is a process on [0, ∞) and is assumed to be independent of the Markov chain. As in [¹⁴,¹⁸-²¹], we assume that the genome is a sequence of contiguous bases and that the distribution of reads along the genome follows a potentially inhomogeneous Poisson process with rate c(x) = c(x, G) at position x. We suppress the argument G unless we take limits. If c(x) = c for all x, we refer to the sampling of the reads as homogeneous. Allowing for inhomogeneity in the sampling Poisson process reflects observations about NGS data in [²²] that the sampling rates of reads at different positions are not homogeneous.

We assume that all sampled reads have the same length of b bps. This length b does not depend on the length of the underlying sequence, nor on the parameters of the Markov chain or the Poisson process. For NGS data, b is given by the technology which is applied to obtain the reads.

A total of M reads are independently sampled from the genome of length G bps. The M reads are our data; they are a realisation of the above Poisson process and the Markov chain sample A, as follows. Suppose that an event in \rho occurs at t. Let s be the closest integer ≤t. Then the Poisson event results in a read R = (As, As+1, ... , As+b ) from the underlying Markov chain A = (A1, ... , AG).

Statistics related to Markov chains based on NGS data

We extend the statistics {\xi }_{ij}^{}, {\zeta }_i^{} and S for long sequence in [¹¹] to {\xi }_{ij}^{(R)}, {\zeta }_i^{(R)} and S(R) for the NGS short sequence/reads data accordingly as follows.

Let f_{ij}^{(R)}be the number of occurrences of dinucleotide (ij) within the M reads and f_i^{(R)}={\displaystyle {\mathrm{\Sigma }}_j{f}_{ij}^{(R)}}. Define

(8){\xi }_{ij}^{(R)}={\displaystyle \frac{f_{ij}^{(R)}-f_i^{(R)}{p}_{ij}}{\sqrt{f_i^{(R)}}}}.

Note that the moment estimation of p_ij from NGS data is {\widehat{p}}_{ij}^{(R)}={\displaystyle \frac{f_{ij}^{(R)}}{f_i^{(R)}}}, thus {\xi }_{ij}^{(R)} can also be written as

{\xi }_{ij}^{(R)}=\sqrt{f_i^{(R)}}({\widehat{p}}_{ij}^{(R)}-p_{ij}).

Similarly, we define

(9){\zeta }_i^{(R)}={\displaystyle \frac{f_i^{(R)}-n_{}^{(R)}{p}_i}{\sqrt{n_{}^{(R)}}}},

and

(10)S_{}^{(R)}={\displaystyle \underset{ij}{\mathrm{\Sigma }}}{\displaystyle \frac{{(f_{ij}^{(R)}-f_i^{(R)}{p}_{ij})}^{\mathrm{2}}}{f_i^{(R)}{p}_{ij}}},

where

n^{(R)}=M\beta ,

with M being the number of reads.

To see the centering in the above equations, we first calculate the expectation of f_i^{(R)} and of f_{ij}^{(R)}. For this purpose we introduce some notations. For a read Rk = (r_k,1, . . . , r_k,b ), we denote its starting position on genome G by S(Rk ), rounded to the lower integer. The left-hand end of the read is a Poisson process which creates the read starts at t such that s≤t<s + 1 for a nonnegative integer s, then we set S(R_k ) = s. Let f_i be the number of occurrences of letter i in the genome, we have that

f_i={\displaystyle \underset{s=\mathrm{1}}{\overset{G}{\mathrm{\Sigma }}}}\mathrm{1}(A_s=i),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbb{E}{f}_i=G{p}_i,

where 1 denotes the indicator function which is 1 when the argument is true, and 0 otherwise.

To calculate the counts of letter i in the M reads, the f_i^{(R)} can be written in the form of S(R_k ) as follows

\begin{array}{l}{f}_i^{(R)}={\displaystyle \underset{s=\mathrm{1}}{\overset{G-\beta }{\mathrm{\Sigma }}}}{\displaystyle \underset{k=\mathrm{1}}{\overset{M}{\mathrm{\Sigma }}}}\mathrm{1}(S(R_k)=s){\displaystyle \underset{\ell \mathrm{=1}}{\overset{\beta }{\mathrm{\Sigma }}}}\mathrm{1}(r_{k,\ell }=i)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle \underset{s=\mathrm{1}}{\overset{G-\beta }{\mathrm{\Sigma }}}}{\displaystyle \underset{k=\mathrm{1}}{\overset{M}{\mathrm{\Sigma }}}}\mathrm{1}(S(R_k)=s){\displaystyle \underset{\ell \mathrm{=1}}{\overset{\beta }{\mathrm{\Sigma }}}}\mathrm{1}(A_{s+\ell -\mathrm{1}}=i).\end{array}

By the independence of the Poisson process and the Markov chain, we have

\begin{array}{l}{{{\displaystyle Ef}}_i}^{(R)}={{\displaystyle p}}_i\beta \underset{s=\mathrm{1}}{\overset{G-\beta }{\Sigma }}E\left(\underset{k=\mathrm{1}}{\overset{M}{\Sigma \mathrm{1}(S({{\displaystyle R}}_k)=s)}}\right)\\ ={{\displaystyle p}}_i\beta \underset{s=\mathrm{1}}{\overset{G-\beta }{\Sigma }}E(\rho (s,s+\mathrm{1}))\\ ={{\displaystyle p}}_i\beta {{\displaystyle \int }}_{\mathrm{0}}^{G-\beta }c(x)dx\\ =(G-\beta )\overline{c}\beta {{\displaystyle p}}_i=M\beta {{\displaystyle p}}_i=n^{(R)}{{\displaystyle p}}_i,\end{array}

where \bar{c} is the averaged coverage, which is calculated as follows

\bar{c}={\displaystyle \frac{\mathrm{1}}{G-\beta }}{\displaystyle {\int }_{\mathrm{0}}^{G-\beta }}c(x)dx=M/(G-\beta ).

This explains the centering in Eq. 9.

For the centering in Eq. 8 we first calculate

\begin{array}{l}{f}_{ij}^{(R)}={\displaystyle \underset{s=\mathrm{1}}{\overset{G-\beta }{\mathrm{\Sigma }}}}{\displaystyle \underset{k=\mathrm{1}}{\overset{M}{\mathrm{\Sigma }}}}\mathrm{1}(S(R_k)=s){\displaystyle \underset{\ell \mathrm{=1}}{\overset{\beta -\mathrm{1}}{\mathrm{\Sigma }}}}\mathrm{1}(r_{k,\ell }=i,\hspace{0.17em}{r}_{k,\ell +\mathrm{1}}=j)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle \underset{s=\mathrm{1}}{\overset{G-\beta }{\mathrm{\Sigma }}}}{\displaystyle \underset{k=\mathrm{1}}{\overset{M}{\mathrm{\Sigma }}}}\mathrm{1}(S(R_k)=s){\displaystyle \underset{\ell \mathrm{=1}}{\overset{\beta -\mathrm{1}}{\mathrm{\Sigma }}}}\mathrm{1}(A_{s+\ell -\mathrm{1}}=i,\hspace{0.17em}{A}_{s+\ell }=j).\end{array}

Note that

\mathbb{E}(\mathrm{1}(A_{s+\ell -\mathrm{1}}=i,\hspace{0.17em}{A}_{s+\ell }=j)|A_{\mathrm{1}},\mathrm{...},A_{s+\ell -\mathrm{1}})=p_{ij}\mathrm{1}(A_{s+\ell -\mathrm{1}}=i),

we then have

\mathbb{E}{f}_{ij}^{(R)}=p_{ij}{\displaystyle \underset{s=\mathrm{1}}{\overset{G-\beta }{\mathrm{\Sigma }}}}{\displaystyle \underset{k=\mathrm{1}}{\overset{M}{\mathrm{\Sigma }}}}\mathrm{1}(S(R_k)=s){\displaystyle \underset{\ell \mathrm{=1}}{\overset{\beta }{\mathrm{\Sigma }}}}\mathrm{1}(A_{s+\ell -\mathrm{1}}=i)=p_{ij}{f}_i^{(R)}.

This conditional expectation calculation makes the centering in Eq. 8 reasonable. Note that we do not center by \mathbb{E}{f}_{ij}^{(R)}=M\beta p_i{p}_{ij}. The centering in Eq. 10 follows from the centering in Eq. 8.

The effective coverage d of NGS data

In order to state our results we need to introduce some more notations. For a position s ∈{1, 2, . . . , G}, let L(s) denote the (random) number of reads which cover position s. For n = 1, 2, . . . ,

C(n)=\{s\in \{1,2,\mathrm{...},G\}:L(s)=n\}

denote the set of positions that is covered by exactly n reads. As in [⁹], we also call C(n) the region of coverage n. Note that a region does not need to be contiguous; we group all the parts of the genome with coverage n into the n-th region.

We define

(11)d_n={\displaystyle \frac{\mathrm{1}}{G}}\mathbb{E}|C(n)|,

which is the expected percentage of the genome that has coverage n. We then define the effective coveraged as follows

(12)d=d(G)={\displaystyle \frac{{\displaystyle {\mathrm{\Sigma }}_{n=\mathrm{1}}^{\infty }}{n}^{\mathrm{2}}{d}_n}{{\displaystyle {\mathrm{\Sigma }}_{n=\mathrm{1}}^{\infty }}n{d}_n}}.

To calculate d, we note that

\begin{array}{l}E\left|C(n)\right|=\underset{s=\mathrm{1}}{\overset{G}{\Sigma }}P(L(s)=n)\\ =\underset{s=\mathrm{1}}{\overset{G}{\Sigma }}P(\rho (s-\beta ,s]=n).\end{array}

can be with an explicit expression, since \rho (s-\beta ,s] has a Poisson distribution with parameter {\displaystyle {\int }_{s-\beta }^s}c(x)dx. In the homogeneous case, \rho (s-\beta ,s] has Poisson distribution with parameter bc. Then d_n = P(Poi(bc) = n), where Poi(bc) denotes a Poisson random variable with parameter bc. Thus in the homogeneous case, we have

{\displaystyle \underset{n=\mathrm{1}}{\overset{\infty }{\mathrm{\Sigma }}}}n{d}_n=\mathbb{E}Poi(\beta c)=\beta c,

and

{\displaystyle \underset{n=\mathrm{1}}{\overset{\infty }{\mathrm{\Sigma }}}}{n}^{\mathrm{2}}{d}_n=\mathbb{E}Po{i}^{\mathrm{2}}(\beta c)=\beta c+{(\beta c)}^{\mathrm{2}}.

Thus, in the homogeneous case, d has a simple expression as

(13)d=1+\beta c.

Remark 1. The results in [⁹] also include d, but the notation is slightly different. The rate c(x) in this paper translates into the rate c(x)/b in [⁹] so that bc in this paper is c in [⁹]. We argue that in the setup it is more intuitive to let the rate of the Poisson process not be scaled by the length of the reads, and hence we change in notation.

Asymptotic distributions of statistics related to Markov chains based on NGS data

The previous study [⁹] showed that under the same NGS short reads generative model as in Methods section which models the underlying long genome sequence as a MC and uses the LanderWaterman model for random sampling of the reads, the traditional standardization of word count as

(14)Z_w^{(R)}={\displaystyle \frac{f_w^{(R)}-{\mathbb{E}}_w^{(R)}}{{\widehat{\sigma }}_w^{(R)}}},

where {\left({\widehat{\sigma }}_{\mathrm{w}}^{\left(R\right)}\right)}^2={\mathbb{E}}_{\mathrm{w}}^{\left(R\right)}(1-\frac{f_{{-}_{\mathrm{w}}}^{\left(R\right)}}{f_{{-}_{\mathrm{w}}-}^{\left(R\right)}})(1-\frac{f_{{-}_{\mathrm{w}}-}^{\left(R\right)}}{f_{{-}_{\mathrm{w}}}^{\left(R\right)}}), for NGS short reads does not converge to the standard normal distribution as that for a long sequence. Thus [⁹], proposed the concept of effective coverage d for NGS short reads data given in Eq. 12, and proved that {\displaystyle \frac{Z_w^R}{\sqrt{d}}} approximates to the standard normal distribution.

In this study, we show through simulations that the effective coverage d is also critical in the asymptotic distribution theory of MC transition probability estimators based on NGS short reads. We have the following two propositions for the asymptotic distributions of {\xi }_{ij}^{(R)}, {\zeta }_i^{(R)} and S^(R), all of which are related to MC for NGS short reads data.

Proposition 1. Assuming the underlying genomic sequence follows a stationary and ergodic MC with transition probability matrix P={[p_{ij}]}_{4\times 4}. We also assume that a set of reads by NGS are randomly sampled according to the Poisson process, possibly inhomogeneous, with rate function c(x). Then, as the genome length G tends to infinity, {\left({\displaystyle \frac{\mathrm{1}}{\sqrt{d}}}{\xi }_{ij}^{(R)}\right)}_{i,j\in \mathcal{A}} converge in distribution to a mean zero normal vector with covariance S with entries {\mathrm{\Sigma }}_{i,j;k,l}={\lambda }_{i,j;k,l} where d is the scaling factor given in Eq. 12 and {\lambda }_{i,j;k,l} is given in Eq. 2.

Proposition 2. Under the same assumptions of Proposition 1, as G → ∞,

1. The vector {\left({\displaystyle \frac{\mathrm{1}}{\sqrt{d}}}{\zeta }_i^{(R)}\right)}_{i\in \mathcal{A}} converges in distribution to a mean zero normal vector with covariance matrix

{({\alpha }_{ij})}_{i,j\in A}

.

2. The statistic S^(R)/d has an approximate {\chi }^{\mathrm{2}}-distribution with s(s - 1) degrees of freedom.

Remark 2. In this study, we only show the validity of the two propositions by simulation studies. Currently, we cannot provide rigorous mathematical proofs for the two propositions, due to the complex intercorrelations of words within reads. We leave them as our conjectures for future mathematical justifications.

Theoretical confidence intervals for MC transition probabilities based on NGS reads data

We show in Proposition 1 that, {\xi }_{ij}^R has an asymptotic normal distribution

{\xi }_{ij}^R={\displaystyle \frac{f_{ij}^R-f_i^R{p}_{ij}}{\sqrt{f_i^{(R)}}}}\sim N(\mathrm{0},d{p}_{ij}(1-p_{ij})).

We applied this property to obtain the theoretical confidence interval for p_ij using NGS data. By dividing the numerator and denominator by f_i^{(R)}, we have

{\xi }_{ij}^R={\displaystyle \frac{{\widehat{p}}_{ij}-p_{ij}}{\sqrt{\mathrm{1}/f_i^{(R)}}}}\sim N(\mathrm{0},d{p}_{ij}(1-p_{ij})),

where, {\widehat{p}}_{ij}=f_{ij}^R/f_i^R, and hence,

{\displaystyle \frac{{\widehat{p}}_{ij}-p_{ij}}{\sqrt{d{p}_{ij}(\mathrm{1}-p_{ij})/f_i^{(R)}}}}~N(\mathrm{0},\mathrm{1}).

For 0 <a< 1, let z_1-a be the number such that the area under the standard normal density function to the left of z_1-a is 1 -a. Thus, we have

(15)P\left(-z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}\le {\displaystyle \frac{{\widehat{p}}_{ij}-p_{ij}}{\sqrt{d{p}_{ij}(\mathrm{1}-p_{ij})/f_i^{(R)}}}}\le z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}\right)=\mathrm{1}-\alpha .

Manipulation of the inequalities gives

P\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\gamma }}\left({\widehat{p}}_{ij}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\gamma -\sqrt{\gamma {\widehat{p}}_{ij}(\mathrm{1}-{\widehat{p}}_{ij})+{\gamma }^{\mathrm{2}}}\right)\le p_{ij}\le {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\gamma }}\left({\widehat{p}}_{ij}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\gamma +\sqrt{\gamma {\widehat{p}}_{ij}(\mathrm{1}-{\widehat{p}}_{ij})+{\gamma }^{\mathrm{2}}}\right)\right)=\mathrm{1}-\alpha ,

wher \gamma ={\displaystyle \frac{z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}^{\mathrm{2}}d}{f_i^{(R)}}}. Therefore, we obtain the level 1 -a confidence interval for p_ij as

{\mathrm{CI}}_{\mathrm{1}-\alpha }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\gamma }}\left({\widehat{p}}_{ij}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\gamma -\sqrt{\gamma {\widehat{p}}_{ij}(\mathrm{1}-{\widehat{p}}_{ij})+{\gamma }^{\mathrm{2}}},\hspace{0.17em}{\widehat{p}}_{ij}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\gamma +\sqrt{\gamma {\widehat{p}}_{ij}(\mathrm{1}-{\widehat{p}}_{ij})+{\gamma }^{\mathrm{2}}}\right).

In practice, f_i^{(R)} can be very large, making f_i^{(R)}\gg z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}^{\mathrm{2}}d and \gamma \approx \mathrm{0}. We thus discard terms of resulting in a simplified confidence interval for p_ij as

(16){\mathrm{CI}}_{\mathrm{1}-\alpha }^{\mathrm{NGS}}=\left({\widehat{p}}_{ij}-z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}\sqrt{{\displaystyle \frac{d{\widehat{p}}_{ij}(\mathrm{1}-{\widehat{p}}_{ij})}{f_i^{(R)}}}},\hspace{0.17em}{\widehat{p}}_{ij}+z_{\mathrm{1}-\frac{\alpha }{\mathrm{2}}}\sqrt{{\displaystyle \frac{d{\widehat{p}}_{ij}(\mathrm{1}-{\widehat{p}}_{ij})}{f_i^{(R)}}}}\right).

Parametric bootstrap for obtaining the confidence intervals for transition probabilities and the stationary distribution of MC

Let X be a stationary and ergodic Markov chain with transition matrix P= [p_{ij}]_s×s. Assume that a set of reads are randomly sampled from X according to the Poisson process, denoted by {R₁, R₂, … , R_M}. We can estimate the transition matrix P and stationary distribution \pi by \widehat{P}=\left[{\widehat{p}}_{ij}\right] and \widehat{\pi }=\left[{\widehat{\pi }}_i\right], where,

(17){\widehat{p}}_{ij}=\frac{f_{ij}^{(R)}}{f_i^{(R)}},

(18){\widehat{\pi }}_i=\frac{{\mathrm{\Sigma }}_{j\in S}{f}_{ij}^{(R)}}{{\mathrm{\Sigma }}_{m=1}^M\left(\left|R_m\right|-1\right)},

where |R_m| is the length of read R_m, m = 1, 2, …, M . Moreover,

\widehat{\pi }\widehat{P}=\widehat{\pi },

\underset{i\in S}{\mathrm{\Sigma }}{\widehat{\pi }}_i=1.

When G → ∞,

\underset{G\to \infty }{\lim }\widehat{P}=P,\underset{G\to \infty }{\lim }{\widehat{\pi }}_i={\pi }_i\mathrm{.}

We provide a parametric bootstrap method for finding the confidence intervals for p_ij and p_{ias follows:}

B1 Simulate the NGS data set. Simulate a genome sequence of length G by the stationary and ergodic Markov chain P , and the sequence starts from the stationary distribution p; generate NGS data set by sampling M reads \{R_1,R_2,\cdots ,R_M\} from the simulated genome by the Poisson process. Note that, for a real data set, we skip Step B1 and go directly to B2.

B2 Given a the set of reads of the NGS data set (either a simulated one from Step B1 or a real one), estimate parameters of transition matrix \widehat{P} and stationary distribution \widehat{\pi } by Eqs. 17 and 18.

B3 Generate bootstrap data set. Simulate a new Markov genome sequence of length G based on the estimated parameters \widehat{P} and \widehat{\pi } in B2; and generate a bootstrap data set with M reads \{R_1^{*},R_2^{*},\cdots ,R_M^{*}\} by the Poisson process.

B4 Estimate parameters of transition matrix \widehat{P}{}^{*} and stationary distribution {\widehat{\pi }}^{*} from the simulated bootstrap data set in B3 by Eqs. 17 and 18.

B5 Repeat Steps B3–B4B times to get B estimators {\widehat{p}}_{ij}^{*} s and {\widehat{\pi }}_i^{*}s; For each dinucleotide, order{\widehat{p}}_{ij}^{*} from the smallest to largest values as {\widehat{p}}_{ij}^{*} (1), {\widehat{p}}_{ij}^{*} (2), • • • ,{\widehat{p}}_{ij}^{*}(B), and the level 1 -a confidence interval for p_ij is given by

(19){\mathrm{CI}}_{1-\alpha }^B=\left({\widehat{p}}_{ij}^{*}\left(\lfloor B\times \frac{\alpha }{2}\rfloor \right),\hspace{0.17em}{\widehat{p}}_{ij}^{*}\left(\lceil B\times \left(1-\frac{\alpha }{2}\right)\rceil \right)\right);

For each i, order {\widehat{\pi }}_i^{*} from the smallest to largest values as {\widehat{\pi }}_i^{*} (1), {\widehat{\pi }}_i^{*} (2), … , {\widehat{\pi }}_i^{*} (B), and the level 1 -a confidence interval for {\pi }_i is given by

(20){\mathrm{CI}}_{1-\alpha }^{\mathrm{b}}=\left({\widehat{\pi }}_i^{*}\left(\lfloor B\times \frac{\alpha }{2}\rfloor \right),\hspace{0.17em}{\widehat{\pi }}_i^{*}\left(\lceil B\times \left(1-\frac{\alpha }{2}\right)\rceil \right)\right)\mathrm{.}

References

[1] Almagor, H. (1983) A Markov analysis of DNA sequences. J. Theor. Biol., 104, 633–645

[2] Reinert, G., Schbath, S. and Waterman, M. S. (2005) Statistics on words with applications to biological sequences. In: Applied Combinatorics on Words, Lothaire, M. ed., ch. 6, pp. 268–352 New York: Cambridge University Press

[3] Blaisdell, B. E. (1985) Markov chain analysis finds a significant influence of neighboring bases on the occurrence of a base in eucaryotic nuclear DNA sequences both protein-coding and noncoding. J. Mol. Evol., 21, 278–288

[4] Pevzner, P. A., Borodovsky, M.Y., and MironovA. A. (1989) Linguistics of nucleotide sequences. I: The significance of deviations from mean statistical characteristics and prediction of the frequencies of occurrence of words. J. Biomol. Struct. Dyn., 6, 1013–1026

[5] Hong, J. (1990) Prediction of oligonucleotide frequencies based upon dinucleotide frequencies obtained from the nearest neighbor analysis. Nucleic Acids Res., 18, 1625–1628

[6] Arnold, J., Cuticchia, A. J., Newsome, D. A., Jennings, IIIW. W. 3rd and Ivarie, R. (1988) Mono- through hexanucleotide composition of the sense strand of yeast DNA: a Markov chain analysis. Nucleic Acids Res., 16, 7145–7158

[7] Avery, P. J. (1987) The analysis of intron data and their use in the detection of short signals. J. Mol. Evol., 26, 335–340

[8] Narlikar, L., Mehta, N., Galande, S. and Arjunwadkar, M. (2013) One size does not fit all: on how Markov model order dictates performance of genomic sequence analyses. Nucleic Acids Res., 41, 1416–1424

[9] Ren, J., Song, K., Deng, M., Reinert, G., Cannon, C. H. and Sun, F. (2016) Inference of Markovian properties of molecular sequences from NGS data and applications to comparative genomics. Bioinformatics, 32, 993–1000

[10] Billingsley, P. (1961)Statistical Inference for Markov Processes, vol. 2. Chicago: University of Chicago Press Chicago

[11] Billingsley, P. (1961) Statistical methods in Markov chains. Ann. Math. Stat., 32, 12–40.

[12] Pevzner, P. A., Tang, H. and Waterman, M. S. (2001) An Eulerian path approach to DNA fragment assembly. Proc. Natl. Acad. Sci. USA, 98, 9748–9753

[13] Zerbino, D. R. and Birney, E. (2008) Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res., 18, 821–829

[14] Zhai, Z., Reinert, G., Song, K., Waterman, M. S., Luan, Y. and Sun, F. (2012) Normal and compound Poisson approximations for pattern occurrences in NGS reads. J. Comput. Biol., 19, 839–854

[15] Song, K., Ren, J., Zhai, Z., Liu, X., Deng, M. and Sun, F. (2013) Alignment-free sequence comparison based on next-generation sequencing reads. J. Comput. Biol., 20, 64–79

[16] Sun, F., Arnheim, N. and Waterman, M. S. (1995) Whole genome amplification of single cells: mathematical analysis of PEP and tagged PCR. Nucleic Acids Res., 23, 3034–3040

[17] Daley, T.and Smith, A. D. (2014) Modeling genome coverage in single-cell sequencing. Bioinformatics, 30, 22, 3159–3165

[18] Lander, E. S. and Waterman, M. S. (1988) Genomic mapping by fingerprinting random clones: a mathematical analysis. Genomics, 2, 231–239

[19] Zhang, Z. D., Rozowsky, J., Snyder, M., Chang, J. and Gerstein, M. (2008) Modeling ChIP sequencing in silico with applications. PLOS Comput. Biol., 4, e1000158

[20] Daley, T. and Smith, A. D. (2013) Predicting the molecular complexity of sequencing libraries. Nat. Methods, 10, 325–327

[21] Simpson, J. T. (2014) Exploring genome characteristics and sequence quality without a reference. Bioinformatics, 30, 1228–1235

[22] Schwartz, S., Oren, R. and Ast, G. (2011) Detection and removal of biases in the analysis of next-generation sequencing reads. PLoS One, 6, e16685

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature