In hallucigenia-sparsa/seqtime: Time Series Analysis of Sequencing Data

# Global options
library(knitr)
opts_chunk$set(fig.path="figure_network_inference/")

We start by loading the seqtime library.

library(seqtime)

Next, we generate a "known" interaction matrix randomly.

N=20
A=generateA(N, c=0.1)
rownames(A)=c(1:N)
colnames(A)=rownames(A)

Below, we plot our known interaction matrix. Green entries stand for a positive interaction of species i with species j, whereas red entries stand for a negative interaction of species i with species j (where i and j are row and column indices, respectively). For instance, a competition consists of two red entries for the negative impact of species i on j and species j on i. The diagonal is entirely red, since the interaction matrix is generated by default with negative values on the diagonal, which represent intra-species competition.

plotA(A, header="Known interaction matrix")

We can also visualize the interaction matrix as a network. In contrast to correlation networks, the network is directed. Thus, it can represent asymmetric ecological relationships. For example, consider a parasitic relationship, where species i benefits at the cost of species j. It would be represented in the matrix by a positive entry in row i and column j and a negative entry in row j and column i. In the directed network, this parasitic relationship would correspond to two arcs, one from species i to species j and a reverse one from species j to species i. The negative diagonal entries are represented as self-loops.

network=plotA(A,method="network")
# use igraph's function tkplot for manual layout of the network

Given our interaction matrix, we can now simulate a test time series with the Ricker community model.

out.ricker=ricker(N,A=A)
tsplot(out.ricker,main="Ricker")

Now, we run LIMITS to test how well it can infer the known interaction matrix from the time series. LIMITS was developed by Fisher and Mehta PLoS ONE 2014. LIMITS needs a minute to run and then returns the estimated interaction matrix.

Aest=limits(out.ricker)$Aest

To compare the known and inferred interaction matrix, we can plot them next to each other:

par(mfrow=c(1,2))
plotA(A,header="known")
plotA(Aest,header="inferred")
par(mfrow=c(1,1))

We can also compute the mean cross-correlation of the known and inferred interaction matrix:

crossCor=cor(A,Aest)
mean(diag(crossCor), na.rm=TRUE)

Finally, we can plot a few quality estimators for the interaction matrix inference. The quality plot displays the correlation between current and future time points one to five steps ahead (autocor) and the correlation between the original time series and a time series generated from the inferred interaction matrix step by step (cor). We can see that the predicted time series barely outperforms lag-one auto-correlation.

limitsqual=limitsQuality(out.ricker,A=Aest,plot=TRUE)

For comparison, we now generate a neutral time series with the Hubbell model. This model does not take species interactions into account explicitly (though implicitly there is a competition between all species) and hence does not take the interaction matrix A as an input.

out.hubbell=simHubbell(N=N, M=N,I=1500,d=N, m=0.1, tskip=500, tend=1000)
tsplot(out.hubbell,main="Hubbell")

We infer an interaction matrix from the neutral time series:

Aesth=limits(out.hubbell)$Aest

The quality plot shows that the inferred interaction matrix leads to a high cross-correlation between predicted and original time series, even though the dynamics in the neutral time series is not determined by species-specific interactions (all interactions are negative and of equal strength, since no species has an advantage over another in the competition for a space in the local community). Thus, a high correlation between predicted and original time series can be misleading as an indicator for the quality of the inferred interaction matrix.

limitsqualh=limitsQuality(out.hubbell,A=Aesth, plot=TRUE)

hallucigenia-sparsa/seqtime documentation built on Jan. 9, 2023, 11:53 p.m.