Extended local similarity analysis (eLSA) of microbial community and other time series data with replicates

Pearson’s correlation coefficient-based analysis

Suppose that the time series data for factors X and Y with replicates are measured simultaneously. We denote them as X = X_[1:n][1:m] and Y = Y_[1:n][1:m], where n is the number of samples (time points) and m is the number of replicates. Let X _{i [1:m]} and Y _{j [1:m]}, or, in more abbreviated form, X _i and Y _j , be the vectors containing the m replicates from the i-th time spot of X and the j-th time spot of Y, respectively. The application of Pearson’s Correlation Coefficient (PCC) requires taking the profile means, i.e. and . Then the PCC between X and Y is defined as:

(1)

where , , and are the means of X and Y, respectively. The statistical significance of r is tested by the fact that follows a t-distribution (degree of freedom: v = n – 2, mean: 0 and variance v / (v – 2)) when m = 1. For a pair of non-replicated series where m = 1, PCC is a straightforward and powerful method to test and identify linear relationship between two bivariate normally distributed random variables. It is widely adopted in the literature but with limitations. Specifically, when the real relationships are more complex, for example, the association between the two factors only occurs in a subinterval of the region of interest or the change of one factor has a time-delay in response to the change of another factor. Several methods, including the original LSA method, have been proposed to overcome such difficulties [11, 16].

Local similarity analysis with replicates

P_i+1,j+1 = max{0,P_i,j + S _XY [F(X _i ),F(Y _j )]} and

S _sgn (X,Y) = sgn[P _max (X,Y) – N _max (X,Y)].

The S _max (X,Y) obtained is the maximum local similarity score possible for all configurations of m-replicated time series X and Y within time-delay D. In this extended algorithm, the scalars x _i ’s and y _i ’ s from the non-replicated series in Ruan et al.[11] are replaced by vector functions F(X _i )’s and F(Y _j )’s to handle data with replicates. Alternatively, we can also consider F(X _i )’s and F(Y _j )’s as the same input data for the original algorithm in Ruan et al.[11], except that they are F-transformed data. In addition, this extended LSA framework easily accommodates the original version of LSA without replicates using m = 1 as a special case.

Different ways of summarizing the replicate data

Notice that the only additional component we introduced in the eLSA algorithm is the function F. Many reports have suggested different possible forms for F, and several computational methods have been proposed for summarizing the additional information available from replicates, including the simple average method (abbreviated as ‘simple’) and the Standard Deviation (SD)-weighted average method (abbreviated as ‘SD’), and the multivariate correlation coefficient method [17–19]. However, the result of the multivariate correlation coefficient method from Zhu et al.[17] can be shown to be the same as the ‘simple’ method. Therefore, in eLSA, we used the first two methods. We also propose the use of median in place of average and Median Absolute Deviation (MAD) in place of SD when robust statistics are needed to handle outliers [20]. The corresponding methods are named simple median method (abbreviated as ‘Med’) and MAD-weighted median method (abbreviated as ‘MAD’), respectively.

The ‘simple’ method is, in spirit, to take the mean profiles to represent the replicated series. In practice, we take F to be the simple average of repeated measurements: . The ‘SD’ method, on the other hand, takes the standard deviation of the replicates into account. Here we take F to be the replicate average divided by its standard deviation (SD): . Importantly, this method utilizes the variability information available, and, as such, it is claimed to be better than the ‘simple’ method in estimating the true correlation [18]. However, in order for the ‘SD’ method to be effective, a relatively large number of replicates, m, are needed, e.g., m ≥ 5. For a small number of replicates, the ‘SD’ method may not work well since the standard deviation may not be reliably estimated. Further, if we replace average with median and SD with MAD, we obtain the ‘Med’ method: F(X _i ) = Median(X _i ) and the ‘MAD’ method: , where MAD(X _i ) = Median(|X _i – Median(X _i )|). The two transformations have similar properties as their corresponding average and SD versions, but they are more robust.

Bootstrap confidence interval for the LS score

With replicate data, researchers can study the variation of quantities of interest and to give their confidence intervals. Due to the complexity of calculating the LS score, the probability distribution of the LS score is hard to study theoretically. Thus, we resort to bootstrap to give a bootstrap confidence interval (CI) for the LS score. Bootstrap is a re-sampling method for studying the variation of an estimated quantity based on available sample data [21]. In this study, we use bootstrap to estimate a confidence interval for the LS score. For a given type I error α, the 1 – α confidence interval is the estimated range that covers the true value with probability 1 – α. Thus, for a given number, B, of bootstraps, we construct the bootstrap sample set , where each and are samples with replacement from X _i and Y _j , respectively. The rest of the calculation is the same as that used for the original data, and we obtain . Without the loss of generality, we suppose that these values are sorted in ascending order: . Then, a 1 – α bootstrap CI of S _max can be estimated by , as suggested by Efron et al.[21].

Data normalization

eLSA analyses require the series of factors X and Y to be normally distributed, but this may not be the case in the real dataset. Therefore, through normalization, the normality of the data can be enforced. To accommodate possible nonlinear associations and the variation of scales within the raw data, we apply the following approach [22] to normalize the raw dataset before any LS score calculations. We use F(X _i ) to denote the F-transformed data of the i-th time spot of an variable X _i . First, we take

(2)

Then, we take

(3)

where Φ is the cumulative distribution function of the standard normal distribution. We will take Z = Z_[1:n] obtained through the above procedure as the normalization of X. Therefore, the normalization steps are taken after the F-transformation.

Permutation test to evaluate the statistical significance of LSA association

It is important to evaluate the statistical significance of the LS score measured by the p-value, the probability of observing a LS score no smaller than the observed score when two factors are not associated locally or globally. To achieve this objective, permutation test is used. To perform the test, we fix Y and reshuffle all the columns of X for each permutation. For a fixed number of permutations L, suppose {X⁽¹⁾,X⁽²⁾,…,X^(L)} is the permuted set of X; then the p-value P _L is obtained using

(4)

where I(·) is the indicator function. With large enough number of permutations, we can evaluate the p-value to any desired accuracy.

False discovery rate (FDR) estimation

In most biological studies, a large number of factors need to be considered. If there are T factors, there will be eLSA pairwise calculations, representing its quadratic growth in T. In order to avoid many falsely declared associated pairs of factors, we need to correct for multiple testing. Many methods have been developed to correct for multiple testing and here we use the method by Storey et al.[23] to address this issue. In particular, we report the q-value, Q, for each pair of factors. The q-value for a pair of factors is the proportion of false positives incurred when that particular pair of factors is declared significant.

Computation complexity and implementation

For a single pair of time series, the time complexity for calculating the LS score using the dynamic programming algorithm is O(n), where n is the number of time points. The estimation of the bootstrap confidence interval for the LS score using B bootstraps will need O(Bn) calculations. The estimation of statistical significance for each pair of factors using L permutations will need O(Ln) calculations. Thus, the number of calculations for a full analysis of each pair of factors will be O(BLn). If there are a total of T factors, there are a total of pairs of factors that need to be compared. Thus, the number of calculations for a full analysis of T factors will be in the order of O(T²BLn), which can be computationally intensive.

In summary, the internal support for replicates and the use of CI estimates are the two major methodological enhancements to LSA. The eLSA software, however, also incorporates other new features, such as faster permutation and false discovery rate evaluations and more options to handle missing values. Other implementation details are available from the software documentation.

Simulations and benchmarks

We generated simulated data to show the efficacy of eLSA in capturing time-dependent association patterns, such as time-delayed associations and associations within a subinterval. We also studied the difference between the eLSA inference using the simple average (referred to as ‘simple’) method, the SD-weighted average method (referred to as ‘SD’), the median (referred to as “Med”) method, and the MAD (referred to as ‘MAD’) method.

Time-delayed association

In this case, X and Y are assumed to be positively correlated with a time delay D. For a particular example with D = 3, we assume that (X_j+3,Y _j )’s follows a bivariate normal distribution with mean μ = 0 and covariance matrix , for j = 1,2,…,20, where ρ = 0.8. X _j ’s are assumed to be standard normal for j = 1,2,3. The generated (X _j ,Y _j )’s are further perturbed m times by a measurement disturbance ε _ij : N(0, 0.01) to obtain the m-replicated series. A pair of simulated series is shown in Figure 1a for a typical simulation with m = 5.

We see that the two series closely follows each other if we shift the Y series three units toward right. In this particular example, the PCC is -0.258 (P=0.272) while the LS score using ‘simple’ averaging method is 0.507 with a p-value of 0.006. We did 1000 bootstraps and the 95% bootstrap confidence interval for this particular example is (0.448, 0.549). Therefore, this time-delayed association is only found significant by the eLSA analysis.

The eLSA analysis pipeline

In this subsection, we briefly describe the eLSA analysis pipeline implemented into the eLSA software package, as shown in Figure 2.

Microbial community data analysis

As an immediate application, we applied the eLSA pipeline to a set of real microbial community time series data. This San Pedro Ocean Time Series (SPOTs) dataset, originally reported in Steele et al.[2] and Countway et al.[25], was collected following a biological feature (i.e. the chlorophyll maximum depth) off the coast of Southern California. The bacterial community was analyzed using the ARISA [4] technique and the protistan community was analyzed using the T-RFLP [26] technique. The dataset is composed of monthly sampled data from September 2000 to March 2004, including 40 time points without replicates. We analyzed the dataset with a delay limit of 3 months and 1000 permutations to evaluate the statistical significance of the LSA score. In this dataset, the factor names, including the operational taxonomic units and environmental factors, are previously defined by Steele et al. [2].

If we look at the top five positive and negative absolute highest LS scores from the unique associations (|D| ≤ 1) found by eLSA (Q ≤ 0.05 and P ≤ 0.05, see Table 3), we can see most of them are time-dependent associations, either time-shifted or within a subinterval. The majority of these are, in any case, beyond the capacity of PCC. In addition, eLSA provides more information about its findings. For example, in the table, Bac609 and Bac675 factors are associated with a shift of one and Euk97 and boxy (oxygen) factors are best associated within a time interval of length 21 starting at time point 15 with no delay. This kind of additional information is not easily obtainable from the PCC analysis but very important for further functional analysis. For instance, we construct an association network using all above unique eLSA associations, as shown in Figure 3. The obtained network obviously reveals some interesting dynamics of the microbial community, such as the domination of positive directed associations, the existence of environmental factors as hubs that are associated with many other factors, (e.g. nutrients such as NO₂, PO₄, SiO₃ and oxygen), and the existence of some highly connected clusters formed by certain bacteria or eukaryote groups.

X	Y	LS	Xs	Ys	Len	D	P	PCC	Ppcc	Q	Qpcc
Euk239	Euk269	0.82	1	1	40	0	0	0.09	0.59	0.02	1.00
Bac609	Bac675	0.77	1	2	39	-1	0	0.14	0.41	0.00	1.00
Euk381	Euk462	0.77	1	1	40	0	0	0.44	0.00	0.02	0.11
Euk583	Bac989	0.68	2	1	39	1	0	0.30	0.06	0.02	0.73
Euk229	Euk339	0.57	1	2	39	-1	0	0.05	0.77	0.02	1.00
Euk97	boxy	-0.62	15	15	21	0	0	-0.42	0.01	0.00	0.17
Euk98	boxy	-0.62	15	15	21	0	0	-0.42	0.01	0.00	0.17
Euk109	boxy	-0.62	15	15	21	0	0	-0.42	0.01	0.00	0.17
Euk112	boxy	-0.62	15	15	21	0	0	-0.42	0.01	0.00	0.17
Euk116	boxy	-0.62	15	15	21	0	0	-0.42	0.01	0.00	0.17

The 5 positive and 5 negative highest absolute LS Scores from associations uniquely found by eLSA in the microbial community dataset. The columns in succession are X (first factor), Y (second factor), LS (Local Similarity score), Xs (start of the best alignment in the first sequence), Ys (start of the best alignment in the second sequence), Len (alignment length), D (shift of the second sequence compared to the first sequence, -: X is ahead of Y, + otherwise), P (p-value for the LS score, 0.00 stands for P<0.005), PCC (Pearson’s Correlation Coefficient), Ppcc (P-value for PCC), Q (q-value calculated for P, 0.00 stands for Q<0.005), Qpcc (q-value for Ppcc).

Taking a closer look at one of the topmost ranked association: Bac609 and Bac675 (see Table 3), we found that they are closely following each other with a time shift of one month, where Bac609 precedes Bac675. Further inspection suggests a yearly pattern that recurs with near regularity for this association, such that Bac609 blooms in early springtime each year (time spots 6, 18 and 29 are February, January and March, respectively), and Bac675 blooms one month later (see Figure 4a). From the binning definition in Steele et al. [2], Bac609 is a Bacteroidetes group bacterium while Bac675 is an undefined bacterium. Since these microbial groups are uncultured, this association as well as many others uniquely identified by eLSA provides new insight into their ecological role in the ocean surface waters. Notice there is an unexpected abundance jump at time spot 35 of the Bac675 series. The reason for this outlier however is unknown to us. While such prominent time-delayed associations as the Bac609 and Bac675 are easily visible, we must caution that time-dependent associations could also be too subtle to be viewed directly. Thus, statistical significance can provide a much more reliable guideline.

Gene expression data analysis

Although LSA had its roots grounded in microbial community analysis, the technique can be readily applied to other biological time series data, such as replicated gene expression time series data from microarray and RNA-Seq experiments [27–29]. Here we show an example of applying eLSA to the dauer exit gene expression profile time series data of 446 genes from a C elegans study. The result of the original study suggests that the 446 genes under investigation have similar kinetics in both the dauer exit and the L1 starvation time course [30]. Here we use the dauer exit time series data consisting of 12 hourly time spots, each with four replicates. We analyzed the dataset with a delay limit of 3 hours and with 1000 permutations and 100 bootstraps.

The results are summarized in Table 2. Comparing the C. elegans results to those of the microbial community, we see that gene-gene associations in this network are much denser, since a smaller number of genes end up with a much larger, rather than smaller, number of eLSA significant associations (e.g. 2804 versus 57991 for Q ≤ 0.05 and P ≤ 0.05, see Table 2). Also different is that about 93% of these associations are found by PCC analysis as well. The high congruence between PCC and eLSA analysis may be due to the fact that about 90% of the eLSA findings are without delays, which thus are also amenable to PCC analysis.

Because these genes do not change expression level in both dauer exit and L1 starvation conditions, they are considered as common feeding response genes [30]. However, it is not clear whether they are correlated with each other in expression profiles under the dauer exit condition. To study this, we combined all eLSA and PCC significant associations with Q ≤ 0.05 and P ≤ 0.05, and found the average degree of the resulting association network is around 169, while that of previous microbial community data is around 12. Such high average degree for C. elegans genes shows the high similarity of their expression profiles, which also reflects their intimate functional coordination along the process. Therefore, our result suggests those feeding response genes are likely to be co-expressed under the dauer exit condition.

We next analyzed the unique eLSA associations. These associations form a dense association network themselves with a long-tailed degree distribution, as shown in Figure 5. While the degree distribution peaks at five, the most highly connected gene 48941 has a degree of 189. We also looked at the top 5 positive and 5 negative highest absolute LS scores unique associations by eLSA. Because replicates are available for this dataset, we are able to obtain the bootstrap confidence intervals for the LS score and they are given in Table 4. Interestingly, we found most of the top LS associations involve high degree nodes, such as genes 48941(189), 29494(129), 29504(128), 27993(116), 436287(106), 32607(58), and 51986(52) (degree in parenthesis). These high degree nodes could be regulation hubs in the feeding response pathway. Here we show an example of time-delayed association of gene 32607 and gene 51986 in Figure 4b. In the figure, gene 51986 leads gene 32607 in expression profile change.

We also analyzed all the eLSA associations together, including both unique and non-unique eLSA findings. Though most of the genes are still hypothetical protein coding genes, we do find a group of eukaryotic initiation factors: 30080(eIF-3E), 33683(eIF-3K), 21358(eIF-3D), 33525(eIF-4E), 32503(eIF-1A) and 23975(eIF-2B) in the 446 selected genes. This is as expected because both L1 starvation recovery and dauer exit will increase translation activities and result in high expression level of these genes. In addition, in the translation process, these factors work closely together to form different translation related complexes [31], so their expression levels should be highly correlated with each other. Actually, if we check the associations found by eLSA, we do see these factors form a clique together with all edges being positive associations and statistically significant (see Figure 6). The coherence of the eLSA finding and our biological knowledge shows that eLSA associations do reveal true associations within the biological system. However, as the majority of genes are still hypothetical, a thorough examination for true functional discoveries will require biological experiments.