Identifying local associations in biological time series: algorithms, statistical significance, and applications

Abstract

Local associations refer to spatial–temporal correlations that emerge from the biological realm, such as time-dependent gene co-expression or seasonal interactions between microbes. One can reveal the intricate dynamics and inherent interactions of biological systems by examining the biological time series data for these associations. To accomplish this goal, local similarity analysis algorithms and statistical methods that facilitate the local alignment of time series and assess the significance of the resulting alignments have been developed. Although these algorithms were initially devised for gene expression analysis from microarrays, they have been adapted and accelerated for multi-omics next generation sequencing datasets, achieving high scientific impact. In this review, we present an overview of the historical developments and recent advances for local similarity analysis algorithms, their statistical properties, and real applications in analyzing biological time series data. The benchmark data and analysis scripts used in this review are freely available at http://github.com/labxscut/lsareview.

INTRODUCTION

Identification of ecological and biochemical relationships in time-series data of a microbiome or a transcriptome is commonly achieved by correlation-based co-occurrence and co-expression analyses. However, global correlation measures such as Pearson's or Spearman's correlation are only effective when two time series are globally synchronized and associated. Yet, more complex relationships with dynamic changes are prevalent in real biological data, including localized and time-delayed associations, which have been observed in various fields such as microbial ecology [1–4], molecular biology [5, 6], and functional neuroscience [7, 8].

For instance, the regulatory mechanism in gene expression analysis between ARG2 (acetyl glutamate synthase) and CAR2 (ornithine aminotransferase) depends on the expression level of CPA2 (arginine specific carbamoyl-phosphate synthase). Nonetheless, the global Pearson's correlation coefficient (PCC) expression level between ARG2 and CAR2 is nearly zero [9], which emphasizes the insensitivity of global correlation measures to detect local associations. Although global correlations continue to be valuable tool, it is necessary to formulate more flexible alternative hypotheses to detect the intricate relationships that exist in biological time series.

To address the limitations of global correlations, researchers have introduced numerous methodologies, the prominent one of such is local similarity analysis (LSA) [10–13]. LSA is a local alignment-based method that identifies alignment configurations optimized for local similarity scores between two time series, which enables the detection of local and potentially time-delayed pairwise associations. A related method, known as local trend analysis (LTA) [14–16], a variant of LSA, identifies alignment configurations and correlations between pairs of trend-of-change series, representing up, down, and no-change statuses. Both LSA and LTA also evaluate the statistical significance of the identified associations. Notably, the term local similarity has not only been used in biology [17], but also in the fields of chemistry [18], public health [19], and power engineering [20], differing in meanings. We refer the readers to their respective papers and reviews.

This paper exclusively delves into the context of local similarity within biological time series. Within the realm of biological time series, local similarity analysis plays a pivotal role in uncovering potential relationships between microbial species richness and environmental factors. Local similarity becomes particularly relevant when actual correlations are limited to specific time subintervals or when correlations are subject to time delays [10]. The utilization of local similarity analysis has garnered extensive application in scrutinizing the spatiotemporal evolution of microbial communities across diverse environments. Importantly, the latent relationships between these microbial species cannot be effectively discerned through conventional correlation-based analysis methods.

Other methods exist for identifying diverse types of complex associations [9, 21–23]. For example, the liquid association (LA) technique [9, 22, 24] detects triplet rather than pairwise associations, while the maximal information coefficient (MIC) method identifies unspecified generic pairwise associations [25]. Granger causality analysis allows for identification of potential causal relationships in molecular ecological time-series data [26]. Furthermore, LA combined with LSA can be utilized to detect mediated covariant dynamics to discover the third-party mediator variable [27]. In this review, we focus on the advancement of LSA and LTA, including algorithms, statistical significance evaluation methods, software tools availability and their applications. Papers [28, 29] provide further information regarding related correlation measures and their applications.

LSA is founded on the concept of local sequence alignment, utilizing a dynamic programming algorithm, similar to the Smith-Waterman algorithm [17] to detect the optimal time series alignment configuration that maximizes the local similarity score. Qian et al. first introduced the method for gene expression profile analysis [10]. Based on this method, LSA has been extended to analyze the microbiome data including operational taxonomic unit (OTU) analysis by Ruan et al. [13] and metagenomic data analysis by Xia et al. [12]. Methodological advances have also been made to LSA, including support for replicate series data by Xia et al. [12], theoretical approximations of statistical significance (i.e., p-values) reported by Xia et al. [11] and Durno et al. [30], and more advanced methods for calculating p-values, such as moving block bootstrap LSA (MBBLSA) [31] and data-driven LSA (DDLSA) by Zhang et al. [32].

Accompanied by many other recently developed correlation measurement techniques for extracting meaningful association patterns from time series data, LSA by Ruan et al. [13] and the extended local similarity analysis tool, eLSA, a significantly accelerated implementation by Xia et al. [11, 12] remain advantageous. For instance, Wang et al. [33] developed a count statistic and evaluated its performance in comparison with LSA, the outcome revealed that as noise increased for pairs of correlated subsequences, LSA remained consistently robust and become more effective. Similarly, Tackmann et al. [34] introduced a Flashweave algorithm designed to infer high-resolution interaction networks while it was outperformed by eLSA when benchmarked with realistic data simulated by Weiss et al. [28].

Since 2011, the eLSA tool has been widely used in microbial data analyses, leading to numerous discoveries. For instance, Liu et al. [35] analyzed with eLSA the correlations between bacterial OTUs and phytoplankton in a reservoir and obtained a dynamical interaction network with significant local ecological associations. Similarly, Thiriet-Rupert et al. [36] constructed a gene co-expression network for a strain of microalgae to determine the relationships between differentially expressed transcription factors and mutant phenotypes by eLSA. Lee et al. [37] utilized eLSA to evaluate the time-dependent relationships among rhizome repair, environmental conditions, and bacterial genera in diesel-contaminated soil and obtained the correlation network.

Moreover, eLSA has been employed in studies of diverse subjects, including to investigate the correlation between marine paleontology and microbial communities by Parada et al. [38], to infer the diverse ecological relationships of microorganisms by Jones et al [39], to analyze the association network of microbial taxa of digested sludge by Liang et al. [40], and to identify the complex relationships among marine phytoplankton, bacteria and viruses [41, 42], between viral bacteriophages and their hosts [43] and between the abundance of antibiotic resistance genes and bacterial communities by many other researchers [44].

Additionally, eLSA's capability of identifying correlations with time shifts has also been widely used. For instance, Posch et al. [45] constructed an interaction network between ciliates and phytoplankton with eLSA to study the correlations with a time lag. They detected such associations among ciliates as well as between ciliates and algae. Džunková et al. [46] investigated the correlation among bacterial OTUs in the oral cavity of different individuals, considering the temporal delay. In another example, eLSA was utilized to evaluate the time-delayed correlation between microbial members of an activated sludge and environmental variables [47].

LSA related studies encouraged the development of methods sharing its rationale. One such method is local trend analysis (LTA), as introduced by He et al. [14] and Ji et al. [15], which analyzes the trend-of-change in gene expression profile series. In LTA, a local alignment algorithm is applied to the trend-transformed series of original data. Trend transformation involves taking the value differences of adjacent time points and discretizing them into {1, 0, −1} which represent {up, no-change, down} trending statuses, respectively. Xia et al. [16] developed a theory to approximate the statistical significance of LTA, which was further improved by Shan et al. known as the stationary theoretical local trend analysis (STLTA) [48].

Numerous studies have demonstrated the efficacy of LSA and LTA in inferring time-dependent associations in contrast to traditional global correlation measures. The added time dimension facilitates mining causality behind the delay. Additionally, the highly efficient theoretical approximation of statistical significance, i.e., p-value, has played a critical role in the widespread use of LSA and LTA analyses. These methods have been successfully applied to large-scale next-generation sequencing (NGS)-based datasets that would otherwise be impractical to analyze using much slower permutation-based p-value. In this paper, we provide an in-depth of the relevant algorithms and theories for LSA and LTA, as well as their recent improvements, including eLSA, fastLSA, MBBLSA, DDLSA, and STLTA.

LOCAL SIMILARITY ANALYSIS OF BIOLOGICAL SERIES DATA

Biological time series data

Time-series data are critical resources for studying the dynamics of biological systems. These data may be the gene expression levels in gene regulation studies or environmental and organismal factor levels in ecological studies. The nature of time-series data can be highly variable with values representing different measurements depending on the dataset.

Before conducting LSA on the original dataset, one needs to normalize the data. This step transforms raw values into standardized values with a mean of zero and a variance of one. However, one obstacle in standardizing count-based time series data is the high frequency of zeros. Excluding zero values from normalization procedures may be preferential since such values may be due to limited detection sensitivity or random dropouts other than true absence. Various normalization schemes are available, such as Z-score normalization and percentile normalization. Xia et al. used percentile normalization with sparsity adjustment, which generally provides the best and most robust result across datasets, even for datasets with non-linear associations (see paper [12] for the detailed procedure).

Local alignment and similarity score of biological series data

To explain the LSA process, we consider two factors with normalized time-series data at levels $X=(X_1,X_2\dots X_n)$ and $Y=(Y_1,Y_2\dots Y_n)$⁠. Without loss of generality, the local similarity (LS) score is determined by calculating the maximized absolute value of summation ${\sum }_{k=0}^{l-1}{X}_{i+k}{Y}_{j+k}$⁠, defined as $S=\underset{0\le i,j,l\le n;|i-j|\le \mathrm{D}}{max}\left|{\sum }_{k=0}^{l-1}{X}_{i+k}{Y}_{j+k}\right|$⁠. Here, $i+k\in I=[i,i+l-1]$ and $j+k\in J=[j,j+l-1]$ which represent the aligned intervals of series $X$ and $Y$⁠. The dynamic programming algorithm based on the Smith-Waterman algorithm [17] can identify these intervals, as shown in Algorithm 1.

This local alignment algorithm has unique features compared to Smith-Waterman's algorithm. First, the maximum time delay, $D$⁠, is introduced, which restricts the dynamic programming results to alignments with starting points in$X$ and $Y$ that are maximally $D$ time units apart. $D$ is typically small since delayed associations usually represent immediate responses to recent changes in other factors, with a lag time that is generally short. A small $D$ also effectively reduces the algorithm's search space, reducing the algorithm’s computational and space complexities from quadratic to linear of the input size. Second, other than optimizing one objective, two score matrices, $P$ and $N$⁠, are simultaneously updated for maximization: one for positive and the other for negative associations. The local similarity score, $S$⁠, is calculated as the maximum of all entries in $P$ and $N$ in absolute terms. $S$ is then standardized by taking the square root of the series $\sqrt{n}$⁠.

Local similarity score was initially scaled by $n$⁠, as in refs. [11–13]. However, in retrospect, $\sqrt{n}$ is more desirable when comparing the scores across varied series lengths [26]. Henceforth, we will refer to this $\sqrt{n}$ scaled alignment score as the local similarity score or LS score. It should be noted that unlike Pearson’s or Spearman's correlation, this LS score is not bound within the range of −1 to 1.

Local similarity score of replicated series data.

Assessing variability through replication is important for statistical inference [49–51]. The original LSA method only considered series data without replications. Xia et al. first proposed the extended local similarity analysis (eLSA) method to incorporate series data with replications [12]. The eLSA assumes that each sample has $m$ replications, and a function $F(\cdot )$ is employed to summarize repeated measurements. The modified eLSA dynamic programming algorithm for repeated measures is presented in Algorithm 2.

In this modified algorithm, the scalars $X_i$ and $Y_i$ from the original LSA algorithm are replaced with summarized and transformed values $F\left(X_i\right)$ and $F\left(Y_j\right)$⁠, respectively, which are applied to the $m$-replicated data vectors $X_i$ and $Y_j$ of size 1 by $m$⁠. Therefore, the eLSA framework considers the original LSA as its special case when $m=1$⁠.

Transformation functions to summarize replicates

To model the repeated measure and account for variability in eLSA, Xia et al. introduced the summarizing function $F(\cdot )$⁠. Researchers previously proposed several methods [51–53] to accomplish this task, including the simple mean method (‘Mean’), the Standard Deviation (SD)-weighted mean method (‘SD’) [52], and the multivariate correlation coefficient method [54], which is identical to the 'Mean' method. Moreover, compared to traditional methods, robust statistics do not need the normality assumption. They demonstrate effective resilience against the impact of outliers in the data, offering a more accurate reflection of the actual situation [54]. For instance, when employing the sample mean for overall mean estimation, the presence of an outlier significantly influences the estimation outcome. On the contrary, the sample median exhibits lower sensitivity to data outliers, rendering it considerably more robust than the sample mean [55]. Therefore, robust statistics, such as Median (Med) and Median Absolute Deviation (MAD) are robust replacement to the 'Mean' and 'SD' [56], and the resulting summary functions were included in eLSA and termed the median method ('Med') and the MAD-weighted median method ('MAD'), respectively. Table 1 shows the formula for common summarizing functions.

Bootstrap confidence interval for local similarity score

The bootstrap method provides an empirical confidence interval (CI) for local similarity score $S$⁠. For a given type-I error $\alpha$⁠, the estimated $1-\alpha$confidence interval will cover the true value of $S$ with probability $1-\alpha$. The bootstrap samples ${\stackrel{\sim }{X_{\mathrm{i}}}}^k$ and ${\stackrel{\sim }{Y_{\mathrm{j}}}}^k$ are generated by sampling with replacement from the original data vectors $X_i$ and $Y_j$⁠, respectively, for a pre-specified total bootstrap sample number $B$⁠. The eLSA calculation is performed on the bootstrap samples in the same way as it is performed on the original data, resulting in $B$ bootstrap LS scores. These scores are then sorted in ascending order, such that ${\stackrel{\sim }{S}}^{(k)}=S({\stackrel{\sim }{X}}^{(k)},{\stackrel{\sim }{Y}}^{(k)})$ represents the $k$-th ranked score. The 95% CI for the original score $S$ is inferred as $[{\stackrel{\sim }{S}}^{\left(\lfloor \frac{\alpha }{2}B\rfloor \right)},{\stackrel{\sim }{S}}^{\left(\lfloor 1-\frac{\alpha }{2}B\rfloor \right)}]$⁠.

A statistical theory for the significance of local similarity score

Permutation statistical significance

Non-replicated data

In the original non-replicated data LSA analysis, the statistical significance (p-value) of an LS score was computed by permutation. This process involved the following steps:

(i) For each pair of time series $X$ and $Y$ with LS score $S(X,Y)$⁠, permute the data values in each sequence to generate a permuted sample pair.

(ii) Compute a permuted LS score for the permuted pair with the LSA algorithm (Algorithm 1).

(iii) Repeat steps (i) and (ii) $L$ times to obtain $(X^1,X^2,\dots ,X^L)$ and $(Y^1,Y^2\dots Y^L)$⁠, which represent $L$ permuted samples of $X$⁠, $Y$ and compute their permuted LS scores.

(iv) The p-value is determined as the proportion of permuted LS scores that were at least as great as $S(X,Y)$⁠, i.e.,

where $I(\cdot )$ is the indicator function.

Replicated data

In the case of replicated data eLSA analysis, the permutation test is modified such that the series data $Y$ remains unchanged for each permuted sample while all columns of $X$ are reshuffled. For a fixed number of permutations $L$⁠, and $(X^1,X^2\dots X^L)$ are the permuted samples of $X$⁠, and the p-value $P_L$ is computed by:

Moving block bootstrap

Zhang et al. [31] developed a moving block bootstrap method (MBBLSA) to replace point-wise permutation and compute the statistical p-value of LSA. The MBBLSA calculation process is as follows:

(i) Calculate the LS scores $S(X,Y)$ for each pair of time series $X$ and $Y$⁠.

(ii) Divide the time series into $n-l+1$ continuous time block of equal length, assuming the time series length is $n$ and the block length is $l$ for each block.

(iii) Randomly select a block from all of the blocks and use its values as the first $l$ resamples. Then, continue selecting blocks randomly, using the next $l$ values as the next resample. If the $k$-th permuted sample $X^k=(X_1^k,X_2^k,\dots X_n^k)$ has already generated data points at $n$ time points, discard the data points after $X_n^k$ and proceed to step (iv). Otherwise, continue sampling.

(iv) Calculate the LS scores of $X^k$ and $Y$⁠, and record them as $S(X^k,Y)$⁠.

(v) Repeat steps (iii) and (iv) to obtain $(X^1,X^2,\dots ,X^L)$ and its LS score $(S(X^1,Y),S(X^2,Y),\dots ,S(X^L,Y))$⁠.

(vi) Calculate the bootstrap p-value as in Eq. (2).

Regarding the block length, please refer to the block length selector proposed by Sherman et al. [57]:

where $NI(\bullet )$ is the nearest integer function, and $\hat{a}$ is the estimated value of the autoregressive coefficient obtained by fitting the autoregressive (AR) (1) model with data.

Bottlenecks of permutation and bootstrap procedures

Although the p-value can be computed with any desired degree of accuracy using a large number of permutations, this approach can be computationally costly. The time complexity for calculating the LS score is $\Theta (n)$ at runtime for a single pair of series, where $n$ is the number of time points. The time complexity for computing bootstrap CI via $B$ bootstraps for the LS score is $\Theta (Bn)$⁠. The time complexity for estimating the statistical significance of a pair of series through $L$ permutations is $\Theta (Ln)$⁠. Thus, the time complexity for conducting a full eLSA analysis on a single pair of series is$\Theta (BLn)$⁠. With series data of $T$ factors, the total pairs are $\frac{T(T-1)}{2}$⁠. Then for an all-to-all analysis, the total time complexity will be $\Theta \left(T^{2}BLn\right)$⁠, which is impractical for large-scale datasets with thousands of factors, such as genome-wide multi-omics or metagenomics microbiome datasets. This motivates the demand for the development of theory-based approaches for efficient p-value approximation.

Theoretical statistical significance for local similarity scores

Feller’s theory for the range of random walk excursion

Due to computational complexity, permutation-based procedures are only practical for pairwise LSA analysis with a few hundred factors. Therefore, it is necessary to develop effective theoretical p-value approximations to alleviate this problem. One direction involves adapting Feller's theory for the range of random walk excursions to the settings of LSA. Feller [58] investigated the approximate distribution of the range of the sum of $\mathrm{n}$ random variables with an expectation of zero. Let $Z_i$ be independent and identically distributed (i.i.d.) random variables such that $E\left(Z_i\right)=0$ and $Var\left(Z_i\right)={\sigma }^2$⁠, where $E(\bullet )$ is the mathematical expectation, and $Var(\bullet )$ is the variance. Let $C_n=Z_1+Z_2+\cdots +Z_n$⁠, $M_n=max\{0,C_1,C_2,\cdots ,C_n\}$⁠, and $m_n=min\{0,C_1,C_2,\cdots ,C_n\}$⁠. The excursion range is then defined as $R_n=M_n-m_n$⁠. According to Feller [56], it can be concluded that:

A theory for local similarity score

Xia et al. [11] proposed the LS score, $S_n$⁠, which is computed through the dynamic alignment of a pair of standard normalized series with length n. It was established that the LS score follows exactly the distribution of $\frac{R_n}{\sigma \sqrt{n}}$ when $\sigma =1$⁠. Notably, Eq. (4) involves an infinite sum, and a stopping rule is required for numerical approximation. In Xia et al. [11], an upper bound was established to determine when to terminate the summation for practical p-value calculations. For a specified error limit, $\beta$⁠, one can choose a summation term upper bound, $K$⁠, such that:

Then, the tail probability $P(\frac{R_n}{\sigma \sqrt{n}}\ge \mathrm{x})$⁠, i.e., the p-value of $S_n$⁠, is approximated by:

A generalization of the theory of local similarity scoring

The approach proposed by Xia et al. [11] is effective only for independent and identically distributed (i.i.d.) sequences. However, in practice, many sequences are not i.i.d.. To solve this issue, Zhang et al. [32] introduced the long-term variance ${\omega }^2$ and replaced $\sigma$ with $\omega$ to approximate the p-value using Eq. (4).

Accurately estimating $\omega$ is crucial for evaluating the statistical significance of LS scores for dependent sequences. Zhang et al. [32] applied the data-dependent $AR(1)$ plug-in [59] method to estimate long-term variance as:

where ${\hat{\gamma }}_x(k)$ and ${\hat{\gamma }}_y(k)$ are the sample auto-covariance functions of $X$ and $Y$⁠, respectively. The bandwidth parameter, $b_{\omega }$⁠, should be solved from $b_{\omega }=1.447{\left(\hat{\tau }n\right)}^{\frac{1}{3}}$ [57], where:

$Z_t=X_t{Y}_t$⁠, t$=1,\dots ,n$⁠, and $\bar{Z}=\frac{1}{n}\sum _{i=1}^n{Z}_i$ is the mean of $Z_t$⁠.

An upper bound approach for statistical significance

Durno et al. [30] proposed an alternative upper bound for the LS score p-value with the fastLSA software. They slightly modified the presentation of the LSA algorithm (Algorithm 1) to enable derivation of their bound for the tail distribution of the LS score statistics. They defined $Z_k=X_{i+k}\ast Y_{j+k}$⁠, $k=0,\dots ,min\{n-i,n-j\}-1$⁠. They established the following equivalence to derive the tail probability bounds: $\{\left|S\right|>x\}=\{\left(U_{i,j}\{P_{i,j}>nx\}\right)\cup \left(U_{i,j}\{N_{i,j}>nx\}\right)\}$ (Lemma 1 in ref [30], where $U$ represents the union operation). This equivalence can be used to express the probability of the tail event {$\left|S\right|>x$} (i.e., the targeted p-value) in terms of the probabilities of the tail events {$P_{i,j}>x$} and {$N_{i,j}>x$}.

For aligned sequences starting from $X_i$ and $Y_j$⁠, Durno et al. found that $\left\{Z_k\right\}$ are i.i.d. random variables. Consequently, condition established by Lindeberg [15] was fulfilled for $\left\{Z_k\right\}$ (Lemma 2 in ref. [30]), which means that the variance of a distribution stabilizes as more variables are incorporated, and its tails are pinned down. Essentially, as the time series get longer, the upper bound of the distribution becomes more well-defined and calculable. In other words, if the random variables $\left\{Z_k\right\}$ have zero expectation and finite variance and the Lindeberg's condition holds, its distribution uniformly converges to that of a one-sided standard normal distribution as the length of the sequence becomes infinite (Theorem 3 in ref. [30]). By applying this theorem to the $\left\{Z_k\right\}$ sequence, they obtained the following:

where $G(x)=\sqrt{\frac{2}{\pi }}{\int }_0^x{e}^{-t^2/2}dt$ if $x\ge 0$and $0$ if $x<0$⁠, which represents the cumulative distribution function of one-sided standard normal distribution (Theorem 4 in ref. [30]).

Durno et al. acknowledged the limitations of their approach. Note that Eq. (8) is asymptotic, and $n$ must be substantially large for this bound to be applicable. They suggested that at least 30 time points are needed for reliable LS score p-value bound implementation. Additionally, it should be noted that convergence can vary across datasets, and applying fastLSA to shorter time series can produce non-negligible false positives, despite the bound being conservative.

A comparison of statistical significance approximation approaches

We compared five available open-source LSA software tools that implement either the permutation or theoretical approaches for estimating the statistical significance of LS scores, as presented in Table 2. The tools include eLSA [11, 12], fastLSA [30], LocSim [13], MBBLSA [31] and DDLSA [32]. Both eLSA and fastLSA are implemented in efficient C++ code, while LocSim, MBBLSA, and DDLSA are written in R. eLSA offer both theoretical (Xia et al. [11]) and permutation-based p-value calculations. FastLSA and DDLSA provide the theoretical p-value calculation. Meanwhile, LocSim and MBBLSA only implement the permutation approach. eLSA also includes additional features, such as trend series analysis, replicated data analysis, and bootstrap confidence intervals. Furthermore, we provided Figure 1 as a viable decision tree for selecting appropriate software tools for various tasks or datasets.

Furthermore, we compared the performance of p-value approximation approaches: theory or permutation. The theoretical methods were implemented in eLSA (Xia et al. [11]), fastLSA (Durno et al. [30]) and DDLSA (Zhang et al. [32]). The permutation approaches were implemented in LocSim (Ruan et al. [13]) and MBBLSA (Zhang et al. [31). The LocSim R tool was excluded from this comparison because its permutation approach is the same as that in eLSA but slower.

We obtained the benchmark data prepared by Weiss et al. from the GitHub (https://github.com/wdwvt1/correlations) [28]. For the benchmark purpose of this paper, we utilized a portion of their data which is shown in Table 3. Tables S1-4 (originally Tables 2.6–2.9 in Weiss et al.) simulated the Lotka-Volterra relationship of two species with different constraints. The Lotka-Volterra relationships describe $n$ species' abundances by a system of $n$ differential equations that mimic ecological relationships. The simulated positive Lotka-Volterra relationships in these tables were confounded with simulated random covariates (i.e., OTUs) from the lognormal and gamma distributions as negative controls.

After pre-processing the raw tables, we filtered out sequences with more than 80% time points being zero values. In the final analysis, Tables S1 and S3 contained 4770 pairs of sequences, including 10 pairs of time series with true correlation. Table S2 contained 4366 pairs of sequences, including 6 pairs of time series with true correlation, and Table S4 contained 4285 pairs of total sequences, including 5 pairs of time series with true correlation. We then input these series data into each LSA analytical software tool to calculate the pairwise LS scores between covariates.

We conducted benchmark tests on the two-species Lotka-Volterra (TLV) datasets with the eLSA, fastLSA, MBBLSA, and DDLSA tools. The performance of each tool was illustrated in Figure 2 by the Receiver Operating Characteristic (ROC) curve (the pROC R package). The curves showed that the eLSA theory (theo), eLSA permutation (perm), and DDLSA methods generally outperformed others in the tests. Specifically, the eLSA (theo) method demonstrated the highest overall accuracy compared to all other methods. Furthermore, we evaluated the tools with data simulating many-species Lotka-Volterra relationship and compared the resultant ROC curves, as shown in Supplementary Table S1 and Figure S1. We observed that the DDLSA methods yielded the best performance in those tests.

Furthermore, Table 4 shows a ranking of correctly identified sequence pairs by all the tested software tools. In the table, n represents the total number of OTU pairs calculated by the tool, and the accurate column recorded the number of correctly called pairs by the tool. Among these tools, we observed that the eLSA (theo) method exhibited the best performance and ranked the first in overall accuracy.

We also conducted a comparison of computational efficiency for the four LSA software tools in terms of the total running time, as measured in minutes, from the benchmark tests (see Table 5). The results indicated that all theoretical approaches - eLSA (theo), fastLSA and DDLSA, demonstrated good performance in terms of computational efficiency, showing more than hundreds of times acceleration than permutation. They finished the analysis within seconds when handling hundreds of covariates (OTUs). Our experience showed, the theory-based approaches can scale up to thousands of covariates using current single processor personal computer, while the permutation approaches become too slow to finish. Overall, the benchmark tests demonstrated that eLSA (theo) is the most recommended choice for LSA analysis, considering overall accuracy and efficiency.

Local trend analysis and statistical significance

Local trend series and local trend analysis

A related LSA analysis technique, local trend analysis (LTA) was also implemented in eLSA for time-dependent trend association mining. Recent studies have suggested that the similarity of increasing, stabilizing or decreasing trends along the timeline can be a strong indicator of associations. To address this, trend series analysis was thus developed to analyse the transformed series. Ji et al. [15] explored this approach by coding the trends of $n$ consecutive time points of gene expression profiles into an $n-1$ time point series representing the local trends of increase, no change, and decrease in gene expression levels. They analysed the transformed series using an exhaustive search algorithm to identify possible local associations of specific length or time span between genes. Later, He et al. [14] updated the analysis with a dynamic programming algorithm and a permutation p-value approach, identical to that of LSA, to compute similarity scores and evaluate the statistical significance. Therefore, we refer to the LSA technique when applied to the transformed trend series as local trend analysis or LTA, with its corresponding similarity measure as the local trend (LT) score.

In LTA, the initial step is to discretize the series profile using the changing trend alphabet $\Sigma$⁠, which reflects a set of symbols of interest representing distinct direction changes in states. Typically, a two-letter alphabet such as $\Sigma =\{-1,1\}$ for trend-down and trend-up states, or a three-letter alphabet such as $\Sigma =\{-1,0,1\}$ for trend-down, no-change and trend-up states, is applied. It is also achievable to discretize into a bigger-size alphabet.

For the given original time series $X=(X_1,X_2\dots X_n)$⁠, we transform the 1-by-n vector$X$ to a 1-by-$(n-1)$ trend vector $d_i^X,i=1,2,\cdots ,n-1$ by the following rules: When$X_i\ne 0$⁠,

where $t\ge 0$ is a threshold value for declaring changing trends. When$X_i=0,$  $d_i^X$is simplified as

Statistical significance of local trend score

Similar to LSA, a theory for the statistical significance of the LT score in LTA is necessary. After transformation, the resulting trend series is no longer independent and identically distributed (i.i.d.). As illustrated in Figure 3, it has a dependence structure wherein the dependence relationship is diminishing over time. Xia et al. [16] demonstrated that a first-order Markov chain provides a good approximation of the resulting dependent trend series in practice:

Therefore, the LT score $S$⁠, as determined by LTA, is identical to the LS score, except that it represents the maximum range of random walk excursion of the sum of Markov-like random variables.

The theory by Xia et al. for LS score can be extended to Markov random variables to approximate the p-value of LT score. Daudin et al. [87] explored the distribution of the maximum cumulative sum of first-order Markov chain with values taken from a finite subset of R. Let $\phi$ be the stationary distribution of such Markov random variables $Z_i,i=1,2,\cdots ,$ where $E\phi \left(Z_1\right)=0$ and

Therefore, following Eq. (6), the p-value approximation formula for LT score is:

Eq. (13) is equivalent to Eq. (6), where the only difference being that the standard deviation ${\sigma }_{\phi }$ is no longer 1. ${\sigma }_{\phi }$ must be calculated from the stationary distribution $\phi$ of trend series $d^X{d}^Y$⁠. In the case of a two-letter alphabet case (i.e., $t=0$ and $\Sigma =\{-1,1\}$⁠), an exact solution was obtained for ${\sigma }_{\phi }$ by Xia et al. [16], which is $\sqrt{1.25}$⁠. However, in the case of three-letter alphabet (i.e., $t>0$ and $\Sigma =\{-1,0,1\}$⁠), Monte Carlo simulations are necessary to approximate the stationary distribution$\phi$ numerically, and the spectral expansion technique is required to solve for ${\sigma }_{\phi }$⁠. Xia et al. [16] provides the details for such calculations, as well as the necessary numerical equations to handle multi-letter alphabets.

The analysis conducted by Xia et al. [16] assumed that the original sequence is i.i.d.. Shan et al. [48] mitigated this assumption by proposing the spectral decomposition theory of matrices, which enables the solution of ${\sigma }_{\varphi }^2$ to be expanded to a broader range of scenarios. When $t=0$⁠, it is assumed that $d_i^X$ and $d_i^Y$ are two independent two-state Markov chains, but the transition probability matrices of $d_i^X$ and $d_i^Y$ are assumed to be distinct 2-by-2 matrices $T_X$ and $T_Y$⁠:

where $a_X$ and $a_Y$ are the transition probabilities.

Under the condition that $X_i$ and $Y_i$ are independent, then:

was proposed by Shan et al. [48] as plug-in estimates to Eq. (13) for computing the LTA p-value. When $t\ne 0,d_i^X$ and $d_i^Y$ are three-state Markov chain, the 3-by-3 transition probability matrices of $d_i^X$ and $d_i^Y$ are:

where $b_i$⁠, $d_i$ and $c_i$ (⁠$i=X,Y)$ are the transition probabilities.

If $X_i$ and $Y_i$ are independent,

was proposed by Shan et al. [48] as plug-in estimates to Eq. (13) for computing the LTA p-value.

Shan et al. [48] also studied the mixed-state Markov chain model. When $d_i^X$ is a two-state Markov chain and $d_i^Y$ is a three-state Markov chain (referred to as a mixed-state Markov model), then:

the proposed plug-in estimates to Eq. (13) for computing the LTA p-value.

Shan et al. [48] referred to this method as STLTA. Through simulation, they demonstrated that the Type I error rate of the STLTA method closely aligns with the specified significance level as the number of time points increases, which indicates that the STLTA method is effective asymptotically.

Finally, Table 6 presented a summary of LSA and LTA software tools and their real applications in various data domains and specialties, aiming to provide practitioners with a convenient overview of the wide adoption of these tools.

DISCUSSION AND CONCLUSIONS

Local association analysis methods, such as local similarity and trend analysis (LSA and LTA), are powerful tools for identifying time-dependent associations in biological time series data. These methods are widely applied in hypothesis generation to determine the most relevant interactions in biological systems. Standard dynamic programming algorithms are often employed to find the optimal time series alignment and compute local similarity and trend scores. Advances in theoretical statistical significance approximation have addressed the common computational bottleneck faced by researchers, notably when analysing large-scale NGS datasets, LSA and LTA methods were brought under a common conceptual framework of random walk excursions for theoretical study. Nevertheless, computational and theoretical challenges remain in cases where local association has gaps, all-to-all comparing tens of thousands or even more variables, or analysing local association involves three or more co-factors. Further development of novel algorithms and statistical theories will be necessary enable these analyses in the future.

Currently, multiple software tools are accessible for executing LSA. Table 2 provides a clear comparison of the features of five such tools: eLSA, fastLSA, LocSim, MBBLSA, and DDLSA. Although the older R tool LocSim is limited to the permutation procedure for p-values, the recent Python and C++ implementation, i.e., eLSA, supports both local similarity and trend analysis and provides both permutation and theoretical approximation for p-values. eLSA is particularly suitable for large NGS-based datasets. eLSA also supports replicated series with multiple summarizing functions and bootstrap confidence intervals. Recently developed MBBLSA and DDLSA methods further improved the permutation and theoretical approximation of p-values. Additionally, another C++ implementation, fastLSA, trades off p-value accuracy for faster computing speed.

This benchmark study has shown that eLSA is the preferred software tool for the analysis of time series data in both sparse and dense data scenarios. Weiss et al. [28] also recommended eLSA, along with MIC [25], as the best analysis methods for dense cross-sectional data when the delay $D$ is set to zero. They highlighted the desirable robust properties of LSA (as implemented in eLSA), irrespective of the distribution and sparsity of the dataset. We recommend that readers refer to ref. [28] for a more in-depth understanding of the comprehensive evaluation therein. Moreover, ref. [29] provides a summary of various analytical methods for studying microbial community networks.

ABBREVIATIONS

DDLSA: Data-driven local similarity analysis

eLSA: Extended local similarity analysis.

i.i.d.: Independent and identically distributed.

LA: Liquid association.

LSA: Local similarity analysis.

LTA: Local trend analysis.

MBBLSA: Moving block bootstrap local similarity analysis.

MIC: Maximal information coefficient.

NGS: Next-generation sequencing.

OTU: Operational taxonomic unit.

PCC: Pearson's correlation coefficient.

STLTA: Stationary theoretical local trend analysis.

ACKNOWLEDGEMENTS

We thank Drs. Jed Fuhrman, Jacob Cram and Joshua Steele for helpful discussions. We thank Drs. R. Knight, S. Weiss, and Mr. Van Treuren for developing and sharing the benchmark data tables.

COMPETING INTERESTS

The authors declare no competing interests.

FUNDINGS

This work was supported by grants from the National Natural Science Foundation of China (61873027 to D.A. and L.C.X.), Open Project of the National Engineering Laboratory for Agri-product Quality Traceability (AQT-2020-YB6 to D.A.), and Guangdong Basic and Applied Basic Research Foundation (2022A1515011426 to L.C.X.).

DATA AVAILABILITY

The datasets analyzed by the study were made publicly available by Weiss et al. on GitHub (https://github.com/wdwvt1/correlations).

Author Biographies

Dongmei Ai is a professor in the School of Mathematics and Physics at University of Science and Technology Beijing. Her research is focused on bioinformatics and computational biology.

Lulu Chen is a master’s student in the School of Mathematics and Physics at University of Science and Technology Beijing. Her research interests include bioinformatics and computational biology.

Jiemin Xie is a master’s student in the School of Mathematics, South China University of Technology. His research interest is bioinformatics.

Longwei Cheng is a master’s student in the School of Mathematics and Physics at University of Science and Technology Beijing. His research interests include computational biology and machine learning.

Fang Zhang works as a machine learning specialist at the Shenwan Hongyuan Securities Co. Ltd. He received his PhD in Mathematics at the School of Mathematics, Shandong University. His research interests are statistics and machine learning.

Yihui Luan is a professor at the School of Mathematics at Shandong University. He works in the fields of biostatistics and time series analysis.

Yang Li is a master’s student in the School of Mathematics, South China University of Technology. His research interest is deep learning and its application to biomedicine.

Shengwei Hou is an assistant professor in the Department of Ocean Science and Engineering, Southern University of Science and Technology. His research interests are bioinformatics, microbial ecology and evolution.

Fengzhu Sun is a professor of Quantitative and Computational Biology and Mathematics at the University of Southern California. He works in computational biology, bioinformatics, statistical genetics, and mathematical modeling.

Li Charlie Xia is a professor in the Department of Statistics and Financial Mathematics, South China University of Technology. He specialized in statistical modeling and algorithm development for biomedical data.

References