# Brain-Inspired Spiking Neural Network for Online Unsupervised Time Series
Prediction

Biswadeep Chakraborty  Georgia Institute of Technology
Atlanta, USA
<EMAIL_ADDRESS>Saibal Mukhopadhyay  Georgia Institute of Technology
Atlanta, USA
<EMAIL_ADDRESS>
###### Abstract

Energy and data-efficient online time series prediction for predicting
evolving dynamical systems are critical in several fields, especially edge AI
applications that need to update continuously based on streaming data.
However, current Deep Neural Network (DNN)-based supervised online learning
models require a large amount of training data and cannot quickly adapt when
the underlying system changes. Moreover, these models require continuous
retraining with incoming data making them highly inefficient. We present a
novel Continuous Learning-based Unsupervised Recurrent Spiking Neural Network
Model (CLURSNN), trained with spike timing dependent plasticity (STDP) to
solve these issues. CLURSNN makes online predictions by reconstructing the
underlying dynamical system using Random Delay Embedding by measuring the
membrane potential of neurons in the recurrent layer of the recurrent spiking
neural network (RSNN) with the highest betweenness centrality. We also use
topological data analysis to propose a novel methodology using the Wasserstein
Distance between the persistent homologies of the predicted and observed time
series as a loss function. We show that the proposed online time series
prediction methodology outperforms state-of-the-art DNN models when predicting
an evolving Lorenz63 dynamical system.

###### Index Terms:

spiking neural network, recurrent, STDP, Wasserstein distance, persistent
homologies, online time series prediction

##  I Introduction

Figure 1: Block Diagram showing the Different Steps for Unsupervised
Prediction of Time Series using RSNN

Real-world systems are not static and keep evolving with time, thus creating
the need for an evolving learning system [1]. Therefore, developing online
learning methods that can make predictions by extracting the underlying
dynamics from the observed time series for real-time learning and prediction
of time-varying environments for machine learning (ML) models running at the
edge is critical [2], [3]. However, current supervised learning models predict
future timesteps based on past observations without any knowledge of the
dynamics of the underlying dynamical system [4], [5], [6]. Hence, these models
fail dramatically for systems continuously evolving, especially systems with
sudden data distribution changes. This is because these models require a large
amount of data for training and cannot adapt quickly to emergent internal
dynamics.

In this paper, we propose a completely unsupervised prediction of the time
series, which is quick and robust in adapting to new unseen dynamics. We use
recurrent spiking neural networks for extremely-energy efficient brain-
inspired models that can continually learn from streaming incoming data using
unsupervised plasticity rules. We focus on a recurrent spiking neural network
(RSNN) model trained using unsupervised spike timing dependent plasticity
rules (STDP) [7], [8] to capture the complex temporal correlations. Recent
works [7], [9], [10] have shown such brain-inspired recurrent spiking neural
networks trained with STDP to exhibit promising performance with very few
computations. Hence the model continually learns representations of the
underlying dynamical systems from which the data is generated. We predict by
reconstructing the underlying dynamical systems using random delay embedding
on the neuron activation time series data. This method of unsupervised
forecasting of time series data contrasts with current regression-type methods
used to predict time series data. It is not only data efficient but also helps
in the continual modeling of evolving intelligent agents. The overall model
can be described using the following four-step process:

  1. 1.

We observe the projected time series from an evolving dynamical system. The
time series is then fed into the encoder layer of an RSNN, where the data is
converted to spike streams that are, in turn, fed into the RSNN, a recurrently
connected network of Leaky Integrate and Fire (LIF) neurons, with the synapses
being continually updated using STDP.

  2. 2.

We sample neurons from the RSNN with the highest betweenness centrality
$\mathcal{C}_{b}$, and each of the observed spike streams is decoded to get a
high-dimensional time series data

  3. 3.

We then use the randomly distributed embedding [11] method to reconstruct the
underlying dynamical system and predict the system’s future steps from this
reconstruction.

The detailed process is shown in Fig. 1 and explained in greater detail in
Section III-A. The key contributions of this work are as follows:

  * •

We propose a completely unsupervised online learning model using a spiking
neural network using continual reconstruction of the underlying dynamical
system

  * •

We sample the nodes with the highest $\mathcal{C}_{b}$ for the reconstruction
and show that such sampling is better than using all the neurons. We relate
this sampling to the bottleneck layer of the autoencoder.

  * •

We proposed two different methodologies for the reconstruction - one where the
model is reconstructed using the RMSE error as the metric and the other where
the model is reconstructed using the 1-Wasserstein distance $d_{W,1}$ between
the persistent homologies between the predicted ($\hat{p}_{H}$) and observed
time series ($p_{H}$) as the metric. We show that the Wasserstein-based
reconstruction method always outperforms the root mean square
error(RMSE)-based reconstruction method, highlighting that the persistent
homologies can better capture the underlying dynamics of the system.

  * •

We showed that for a limited amount of data, the RSNN-based models outperform
the DNN-based models and converge faster. However, with multiple repetitions
of the same data, the DNNs outperform the RSNN models.

We evaluated the models on a synthetic evolving dynamical system dataset and
real-world datasets like YReal, and Dow Jones Industrial Average. The proposed
CLURSNN methods outperform the supervised DNN models for online time series
prediction tasks. The rest of the paper is organized as follows: Sections II
and III discusses the Background and the Methods used in the paper, while in
Section IV, we discuss the experiments performed and the results observed.
Finally, we present the conclusions from these observations in Section V.

##  II Background

Related Works:  Recent works on predictive models [4], [5], [6], [12]
demonstrate the ability to discover hierarchical latent representations and
complex dependencies. However, such studies focus on batch learning settings
where the entire training data set should be available a priori, meaning that
the relationship between inputs and outputs is always static. This assumption
is restrictive for real-world applications where data arrives in streams, and
the underlying dynamical system generating the data changes [13]. Cui et al.
[3] described the potential of ML methods for various applications like
traffic profiling, device identification, system security, Internet of Things
applications, edge computing, etc. However, training an artificial
intelligence model at the edge from scratch can be time-consuming and energy-
inefficient. Therefore, it is desirable to train deep predictors online using
only new samples to capture the changing dynamics of the environment [14],
[15].

Despite the ubiquity of online learning in many real-world applications,
training detailed predictors online remains a challenge. First, standard deep
neural networks converge slowly on the data stream due to the unavailability
of offline training advantages such as mini-batch and multi-epoch training
[16], [17]. Overall, deep neural networks have powerful representation
learning capabilities but lack mechanisms for successfully training data
streams. Therefore, online time series forecasting using deep models is a
promising but challenging problem. Several methods have been proposed to solve
the online learning problem. However, current deep learning-based models
require a huge amount of data for training and are also extremely energy
inefficient, making it difficult to deploy in edge devices. Brain-inspired
spiking neural networks have been proposed as the next generation of energy-
efficient neural network models that process information in the discrete spike
domain, unlike standard DNNs[18]. Most recent works on spiking neural networks
are targeted for time series or spatiotemporal classification [19], [9], [20].
Some recent works also use spiking neural networks for time series forecasting
[21], [10]. However, these models are trained using supervised learning
methods [21] or need a supervised readout layer [22], hence cannot adapt to
evolving systems. However, no work exists that uses a completely unsupervised
time series prediction method that reconstructs the underlying dynamical
system, which is essential for edge ML applications.

Recurrent Spiking Neural Network (RSNN)  consists of spiking neurons connected
with synapses. Our simulations use the Leaky Integrate and Fire (LIF) neuron
model. In the LIF neuron model, the membrane potential of the $i$-th neuron
$u_{i}(t)$ varies over time as:

$\tau_{m}\frac{dv_{i}(t)}{dt}=-\left(v_{i}(t)-v_{rest}\right)+I_{i}(t)$ (1)

where $\tau_{\mathrm{m}}$ is the membrane time constant, $v_{rest}$ is the
resting potential and $I_{\mathrm{i}}$ is the input current. When the membrane
potential reaches the threshold value $v_{\text{th }}$ a spike is emitted,
$v_{i}(t)$ resets to the reset potential $v_{\mathrm{r}}$ and then enters a
refractory period where the neuron cannot spike. Spikes emitted by the $j$th
neuron at a finite set of times $\left\\{t_{j}\right\\}$ can be formalized as
a spike train $\displaystyle S_{i}(t)=\sum\delta\left(t-t_{i}\right)$. We use
a STDP rule [7] for updating a synaptic weight ($\Delta w$) and is given by :

$\displaystyle\Delta w(\Delta
t)=\\{\begin{array}[]{l}A_{+}(w)e^{-\frac{|\Delta
t|}{\tau_{+}}}\quad\quad\quad\text{ if }\Delta t\geq 0\\\
-A_{-}(w)e^{-\frac{|\Delta t|}{\tau_{-}}}\quad\quad\quad\text{ if }\Delta
t<0\end{array}$ (4)
$\displaystyle\text{s.t.}A_{+}(w)=\eta_{+}(w_{\max}-w),A_{-}(w)=\eta_{-}(w-w_{\min})$

where $\Delta t=t_{\text{post }}-t_{\text{pre }}$ is the time difference
between the post-synaptic spike and the presynaptic one, with synaptic time-
constant $\tau_{\pm}$.The RSNN is made of three layers: (1) an input encoding
layer ($\mathcal{I}$), (2) a recurrent spiking layer ($\mathcal{R}$), and (3)
an output decoding layer ($\mathcal{O}$). The recurrent layer consists of
excitatory and inhibitory neurons, distributed in a ratio of
$N_{E}:N_{I}=4:1$. The PSPs of post-synaptic neurons produced by the
excitatory neurons are positive, while those produced by the inhibitory
neurons are negative. We used a LIF neuron model and trained the model using
STDP rules. Fig. 2 shows an RSNN with LIF neurons and STDP synapses.

Figure 2: Figure showing the Recurrent Spiking Neural Network with LIF neurons
and STDP synapses. The inset figures show the dynamics of the LIF neurons and
the STDP weight dynamics.

Randomly Distributed Embedding (RDE)  algorithm was proposed by Ma et al. [11]
to accurately predict future states based on short-term high-dimensional data.
RDE accurately predicts future dynamics based on short-term high-dimensional
data. To compensate for the limited number of time points, the RDE framework
creates the distribution of information from the interactions among high-
dimensional variables. Thus, instead of roughly predicting a single trial of
future values, RDE achieves accurate prediction using the distribution
information. The RDE framework randomly generates sufficient low-dimensional
non-delay embeddings from the observed data of high-dimensional variables.
Each of these non-delay embeddings is then mapped to a delay embedding
constructed from the data of the target variable. A target variable is a
variable that we aim to predict using the RDE framework. Any of these mappings
can perform as a low-dimensional weak predictor for future state prediction,
and all such mappings generate a distribution of predicted future states. This
distribution patches all pieces of association information from various
embeddings into the complete dynamics of the target variable. After being
operated by appropriate estimation strategies, it creates a stronger predictor
for achieving prediction in a more reliable and robust form. Our proposed
method performs the RDE-based state-space reconstruction method shown in Fig.
1 every time the batch error crosses some threshold. We use two metrics to
calculate batch error: the RMSE error between the predicted and the observed
time series and the Wasserstein distance between the persistent homologies of
the two $d_{W,1}(\hat{p}_{H},p_{H})$, hereafter denoted as $d_{W}$. The
details about the two loss metrics are further discussed in Section III-B.

###  II-A Topological Data Analysis and Persistent Homology:

Topological Data Analysis (TDA) [23] uses information from topological
structures in complex data for statistical analysis and learning to identify
shape-like structures in data. This paper uses a mathematical tool used in TDA
called persistent homology, a part of the computational topology that converts
data into simplicial complexes. A simplicial complex is a space with a
triangulation. Formally, a simplicial complex $K$ in $\mathbb{R}^{n}$ is a
collection of simplices in $\mathbb{R}^{n}$ such that (i) every face of a
simplex of $K$ is in $K$, and (ii) the intersection of any two simplices of
$K$ is a face of each of them. Thus, persistent homology computes the birth
and death of such topologies via a persistence diagram. The procedure to
compute persistent homology associated with the input point cloud data set
involves the construction of filtration of simplicial complexes, ordered with
respect to some resolution (scaling) parameter. A topological feature that
persists for a more extensive range of scales is a significant one. An
important reason why we use persistent homology is that it does not require an
artificial cutoff between signal and noise since all topological features from
the data are preserved, and weights are assigned according to their
persistence. The output of this filtration procedure generates the persistence
diagram. In the persistence diagram, the two coordinates of each point
represent the birth value and the death value of a k-dimensional hole. A key
property of persistent homology is that both persistence diagrams are robust
under perturbations of the underlying data. The shift in the persistence
diagram is directly proportional to any change in the dataset. These
properties are leveraged for developing statistical methods for data analysis
using persistent homologies. TDA provides a new type of analysis that
complements standard statistical measures like RMSE. We use TDA to study the
shape of the time series data by calculating $d_{W}$ to detect and quantify
topological patterns that appear in multidimensional time series [24].
Following the works of Gidea et al. [24], we apply TDA methods to explore the
temporal behavior of topological features in the dynamical system time series
data as the state of the dynamical system evolves. Using a sliding window, we
extract time-dependent point cloud data sets to which we associate a
topological space. We detect transient loops in this space and measure their
persistence, encoded in real-valued functions called a persistence landscape.

##  III Methods

###  III-A General Framework

There are many components in the unsupervised online learning methodology.
Figure 1 visualizes the complete data flow of the model to be implemented.
First, the evolving time series is derived by projecting the Lorenz63 (L63)
dynamical system [25] onto a single axis. The L63 system works in four modes:
fixed point, chaos, limit cycle, and normal. Further details about the L63
dynamical system are given in Section IVa. We consider the underlying
dynamical system shifting from one mode of operation to another to result in
an evolving time series as depicted in Fig. 4(a). This time series is then fed
as input into a layer of encoder neurons for converting the signal into spikes
using a temporal encoding scheme. These spike trains are given as input to the
recurrent spiking neural network (RSNN), which is trained continually using
spike timing dependent plasticity rule (STDP) [9], [10]. Hence, we use the CS-
HiBet method [26] to efficiently detect top-k nodes in networks with the
highest $\mathcal{C}_{b}$. CS-HiBet can perform as a distributed algorithm by
using only the local information at each node, which makes it highly scalable.
We record the membrane potential of these K nodes and reconstruct the
underlying dynamical system from this data using random delay embedding (RDE)
[11] from which we get the final predictions.

Encoder  encodes the time series data into a spike train. For this paper, we
use a population Poisson encoder, i.e., each value of the continuous input
time series signal is converted into a series of spike trains. The signal must
be represented as spikes for the RSNN to process our time series to higher
dimensional output. The spikes are only generated whenever the signal changes
in value. This paper uses a temporal encoding mechanism to represent the
signals efficiently. We implement the encoding mechanism based on the Step-
Forward (SF) algorithm [27]. The property of temporal encodings to transfer
required information faster than rate encoding can be helpful for online
learning. Online adapting models that have to make quick decisions based on
changes in the environment benefit most from temporal encoding.

Neuron Sampling: The proposed algorithm is based on the hypothesis that the
RSNN captures the dynamics of the underlying system from the observed input
time series data. Thus, a reconstruction of the dynamical system from the
neuron states can predict the future states of the dynamical system. In this
work, we use a neuron sampling algorithm based on $\mathcal{C}_{b}$ and select
the neurons with the highest information flow. $\mathcal{C}_{b}$ is used to
detect the degree of influence a node has over the flow of information in a
network and to find nodes that serve as a bridge between two different parts
of the network. Hence, we can interpret the nodes with high $\mathcal{C}_{b}$
similar to the bottleneck layer neurons in an autoencoder that showcase a
lower-dimensional embedding of the higher-dimensional input space.
$\mathcal{C}_{b}$ measures the proportion of shortest paths in the network
passing through a specific node. This is the fraction of times a node acts as
a bridge in transferring valuable information between any pair of nodes along
their shortest paths within the network. Hence, $\mathcal{C}_{b}$ is
correlated to the nodes’ importance from an information-flow standpoint in the
network. Following previous works [28], we define $\mathcal{C}_{b}$ of node
$u\in V$ as: $\displaystyle\mathcal{C}_{b}=\sum_{v,w,v\neq
w}\frac{\sigma_{vw}(u)}{\sigma_{vw}}$, where $\sigma_{vw}$ is the total number
of shortest paths between every $v,w\in V,v\neq w$, and $\sigma_{vw}(u)$ is
the number of such paths that pass through $u$.

Decoder We take the membrane potential outputs of the sampled nodes and feed
them into the decoder to output a high-dimensional time series. To represent
the RSNN state $r(t)$, we take the sum of all the spikes $s$ over the last
$\tau$ time steps into account. $\displaystyle
r_{i}^{M}(t)=\sum_{n=0}^{\tau}\gamma^{n}s_{i}(t-n)\quad\forall i\in E$. By
definition, this is an exponentially decreasing rate of decoding using a
sliding window. A time series that is highly volatile in a short time frame
requires the decoder to react fast to the changes. The parameters $\tau$ and
$\gamma$ are balanced to optimize the memory size of the stored data and its
containment of information, including adjusting $\tau$ to the pace at which
the temporal data is presented processed. Since we model a continually
evolving system, we chose the values of $\tau\leq 50$ and $\gamma^{\tau-1}\geq
0.1$.

###  III-B Real-time Error Computation

We note that the reconstruction of the dynamical system is done every time
some loss metric between the predicted and observed time series exceeds some
threshold. A more frequent reconstruction would make a system inefficient and
prone to noise and outliers. For the experiments, we used an RSNN with 5000
LIF neurons and a synapse connectivity probability of 20%. However, though the
model was the same for all the simulations, we used two different loss
functions for reconstructions as follows:

RMSE-CLURSNN: We use the RMSE loss as a metric to evaluate the loss between
the predicted and the observed time series. The RMSE loss is given as
$\text{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2}}$
where $y$ is the observed time series and $\hat{y}$ is the predicted time
series. We calculate the RMSE loss on a rolling window of the observed and
predicted time series.

Wass-CLURSNN:  Though the standard RMSE error is an efficient loss metric, it
fails to capture whether the predicted time series successfully captures the
shape of the data. When using RMSE, we can get a high value even though the
predicted time series captures the structure of the data. Since we aim to
reconstruct the underlying dynamical system, we are more interested in whether
we can capture the dynamics successfully, which might lead to an amplified or
time-shifted prediction. Though these types of distortions do not change the
fact that the model has successfully learned the underlying model, it leads to
a high RMSE loss, which is undesirable. Hence, to bypass this problem, we
introduce a novel loss metric in this paper. We introduce $d_{W}$ as an
alternate loss metric instead of the standard RMSE error. As described before,
the persistent homologies capture the topological features of the data and not
just the point values, unlike RMSE. $d_{W}$ is calculated on a rolling time
window. The detailed procedure for calculating $d_{W}$ is shown in Fig. 3. A
persistence diagram is a multiset of points in the extended plane,
$\overline{\mathbb{R}}^{2}$. Let $X$ and $Y$ be two persistence diagrams. To
define the distance between them, we consider bijections $\eta:X\rightarrow Y$
and take the sum of $q$-th powers of the $L_{\infty}$-distances between
corresponding points, and minimize overall bijections. Measuring the distance
between points $x=\left(x_{1},x_{2}\right)$ and $y=\left(y_{1},y_{2}\right)$
as
$\|x-y\|_{\infty}=\max\left\\{\left|x_{1}-y_{1}\right|,\left|x_{2}-y_{2}\right|\right\\}$
and taking the infimum over all bijections, we get the q-Wasserstein distance
$d_{W.q}$ as: $\displaystyle d_{W,q}(X,Y)=\left[\inf_{\eta:X\rightarrow
Y}\sum_{x\in X}\|x-\eta(x)\|_{\infty}^{q}\right]^{1/q}.$

The RSNN model is kept the same for either case; hence, the learning procedure
does not change. We only use the error function as a metric for efficient
reconstruction of the dynamical system attractor space based on the sampled
nodal observations of the RSNN.

Figure 3: Block Diagram showing the computation of Wasserstein Distance
between the persistent homologies of the observed and the predicted time
series

##  IV Experiments

###  IV-A Dataset

Synthetic: Evolving Lorenz63 (L63) System: We study an evolving Lorenz system
using the RSNN models. The classical three-dimensional L63 attractor [25] was
the first example of a low-dimensional system with chaotic solutions. The
following differential equations describe the L63 system:

$\frac{\mathrm{d}x}{\mathrm{~{}d}t}=\sigma(y-x),\frac{\mathrm{d}y}{\mathrm{~{}d}t}=x(\rho-z)-y,\frac{\mathrm{d}z}{\mathrm{~{}d}t}=xy-\beta
z$ (5)

We consider four different types of behavior of the L63 system depending on
the parameter values:

  * •

Fixed Points  (F.P) $(\rho=60,\sigma=20,\beta=8)$

  * •

Chaos  $(\rho=36,\sigma=8.5,\beta=3.5)$

  * •

Limit Cycles  (L.C.) $(\rho=35,\sigma=21,\beta=1)$

  * •

Normal  $(\rho=28,\sigma=10,\beta=2.66)$

Real-world Datasets:  YReal Datasets:  [29] The YReal time series is the real-
world timeseries from the A1Benchmark data in Yahoo’s S5 dataset. We use the
time series in A4Benchmark containing synthetic outliers and change points.
The time series represents the metrics of various Yahoo services with
different change points and outliers. Since our model is learning from
streaming data, the model’s performance on outliers can help us determine the
robustness of the methods.

Dow Jones Industrial Average: We use the daily index of the Dow Jones
Industrial Average (DJIA) from 1885–1962.

###  IV-B Methodology

Figure 4: (a) Fig showing the four modes of operation of the L63 dynamical
system. Each mode has 300 timesteps. (b) The predicted and the observed time
series data for the three models compared in the paper - RMSE-CLURSNN, Wass-
CLURSNN, and online DNN. (c) The RMSE loss between the predicted and the
observed time series computed in a rolling time window of 30 timesteps for the
three models (d) $d_{W}$ of the three models (e) The predicted and observed
time series of the time series with a sinusoidal trend. We only show the first
300 timesteps (corresponding to the FP of the L63 System). A similar process
was repeated for the other three modes of operation.

Baselines DNN Model: We compare the performance of the novel SNN-based
unsupervised learning methodology with a supervised DNN-based online learning
model. We choose the conventional real-time current learning (RTRL) of LSTM
neural networks[30], henceforth called online DNN. It updates the network
using the newly available data without considering the change points. The
online DNN model uses a stochastic gradient descent (SGD) based method to
learn the model from the streaming time series. With new observations, the
model parameters are updated online according to the gradients of the loss of
the newly available data. Online DNN explores the local features of time
series to weigh the learning method’s gradients with the local data’s
distributional properties. In particular, if a newly available observation is
a potential outlier behaving differently from the regular patterns, the
corresponding gradient will be down-weighted to avoid abruptly leading the
online model to drift from the underlying patterns.

###  IV-C Results

Prediction Performance:  We compare the online prediction performances of the
RSNN models with the baseline DNN models. Fig. 4(b) shows the prediction
results of the Wass-CLURSNN, RMSE-CLURSNN, and the online DNN models compared
to the observed time series. The Wass-CLURNN can follow the evolving system
better than the other two models. The online DNN model shows the worst
convergence among the three. These results are further validated by the plot
of RMSE and $d_{W}$ shown in Fig. 4(c, d), respectively. From these plots, we
see that the Wass-CLURSNN performs the best and converges the fastest among
the three, while the online DNN, i.e., the LSTM model, performs the worst and
converges the slowest. It is to be noted here that the RMSE and $d_{W}$ are
calculated here on a rolling window of 30 timesteps.

Adaptation to Timeseries with Sinusoidal Trend:  We evaluate the rate at which
the different methods used in this paper can adapt to a changing dynamical
system. To make the system more complex and study this adaptation behavior
further, we use a sinusoidal signal as a trend line for the time series from
the L63 system sinusoidal oscillating within a sinusoidal function. We call
this process the sinusoidal trend of the time series. The frequency of the
trendline sinusoidal signal determines the time series’ evolution rate, where
a higher frequency (smaller Time period) implies the time series changes
rapidly and the model needs to adapt more quickly to the changing system.
Studying the prediction performance of the three models (Wass-CLURSNN, RMSE-
CLURSNN, and online DNN) on this sine-modulated L63 time series, we test how
each of these models can adapt to changing dynamical systems. Fig 4(e) shows
an instance of the first 300 timesteps of the sine-modulated time series. The
figure shows where the L63 system in the fixed point mode is modulated using a
sine wave of the time period of 300 timesteps and an amplitude of 3. Table I
shows the average RMSE and average $d_{W}$ for each of the three methods. The
mean and standard deviation are taken over each of the four modes of
operation. First, we report the baseline case for an unmodulated time series.
We also give the results of sine-waves with periods of 100, 300, and 500 with
amplitude of 3 and 5 for each case. The results are summarized in Table I. We
see that Wass-CLURSNN outperforms RMSE-CLURSNN and online DNN models.

TABLE I: Table comparing the performances of the three models - RMSE-CLURSNN,
Wass-CLURSNN, and online DNN when using different sine-modulated time series
of the evolving L63 dynamical system. The experiment is repeated for each of
the four modes of operation, and the mean and standard deviation are
calculated from these observations.

Trendline Sine Amplitude (A) |  Trendline Sine Time Period | RMSE-CLURSNN | Wass-CLURSNN | online DNN
---|---|---|---|---
| Average
---
RMSE
| Average
---
$d_{W}$
| Average
---
RMSE
| Average
---
$d_{W}$
| Average
---
RMSE
| Average
---
$d_{W}$
0 | - | $0.324\pm 0.021$ | $0.301\pm 0.015$ | $0.206\pm 0.02$ | $0.166\pm 0.013$ | $0.913\pm 0.064$ | $0.628\pm 0.044$
3 | 100 | $0.615\pm 0.037$ | $0.352\pm 0.019$ | $0.441\pm 0.031$ | $0.198\pm 0.013$ | $1.554\pm 0.11$ | $0.954\pm 0.067$
300 | $0.596\pm 0.031$ | $0.326\pm 0.014$ | $0.322\pm 0.024$ | $0.183\pm 0.011$ | $1.388\pm 0.094$ | $0.799\pm 0.051$
500 | $0.357\pm 0.025$ | $0.311\pm 0.011$ | $0.217\pm 0.019$ | $0.174\pm 0.01$ | $1.068\pm 0.075$ | $0.701\pm 0.042$
5 | 100 | $0.719\pm 0.053$ | $0.434\pm 0.021$ | $0.507\pm 0.048$ | $0.276\pm 0.019$ | $1.994\pm 0.139$ | $0.973\pm 0.061$
300 | $0.655\pm 0.049$ | $0.388\pm 0.019$ | $0.389\pm 0.044$ | $0.221\pm 0.016$ | $1.635\pm 0.111$ | $0.882\pm 0.598$
500 | $0.433\pm 0.032$ | $0.335\pm 0.016$ | $0.294\pm 0.029$ | $0.195\pm 0.014$ | $1.121\pm 0.079$ | $0.866\pm 0.554$

Figure 5: Figure compares the convergence ability of the three models over
repeated iterations of the same input. The convergence is calculated using the
mean number of timesteps required for the model to converge in each case. The
experiment was repeated for each of the four modes of operation of the L63
system, and the mean and standard deviation from these observations are
plotted. We use the log-scale value for the y-axis.

Memory Retention and Convergence:  We know that the time series constitutes
four different modes of operation of the L63 dynamical system. The time series
observed for these different modes of operation is shown in Fig. 4(a). Let us
denote these behaviors as $a^{x}$ where $x$ can be $FP=$Fixed Point , $C$=
Chaotic, $LC$=Limit cycles, $N$= Normal. Also, let us denote the prediction
time series as $a_{n}$. We say that $a_{n}$ has converged to a value L if
$\forall\varepsilon>0\exists N$ such that $n\geq
N\Rightarrow|a-L|<\varepsilon$. For our simulations, we fix the value of
$\varepsilon=0.1$. We use a single training epoch for the supervised models to
keep the results comparable to unsupervised SNN models. Hence, to see how the
convergence behavior changes when the same input is repeated, we perform the
following experiment: since each iteration constitutes four modes of
operation, we take 300 timesteps of each mode and calculate the number of
timesteps the models require to converge for each mode. Let the number of
timesteps required for iteration $i$ for mode $x$ by the model $\mathcal{M}$
be denoted by $T_{x}^{i}(\mathcal{M})$. Hence, we compare the average of the
number of timesteps required to converge over all the modes of operation for
each iteration, i.e., $\sum_{x}T_{x}^{i}(\mathcal{M})$. The results are
plotted in Fig. 5. The figure shows that the Wass-CLURSNN outperforms the
other methods in the first few iterations and converges the fastest. We also
see that the online DNN model performs the worst among the three, indicating
that when training, the DNN-based model requires more training samples to
converge compared to the unsupervised RSNN-based models. However, as the
number of iterations increases, the online DNN takes much fewer timesteps.
From this observation, we can infer that DNN-based models can converge much
faster once the models are fully trained. Thus, the DNN-based models fail when
the underlying dynamical system is fast evolving and there is not a lot of
data available, and that is when the RSNN-based models perform the best. Thus,
these two models have two separate regions of operations.

Metric Model | RMSE-CLURSNN | Wass-CLURSNN
---|---|---
| All
---
Neurons
| Sampled
---
Neurons
| All
---
Neurons
| Sampled
---
Neurons
| RMSE
---
(Mean $\pm$ Std. Deviation)
| $0.544$
---
$\pm 0.112$
| $0.324$
---
$\pm 0.021$
| $0.411$
---
$\pm 0.132$
| $0.206$
---
$\pm 0.02$
| Wasserstein Distance
---
(Mean $\pm$ Std. Deviation)
| $0.457$
---
$\pm 0.854$
| $0.301$
---
$\pm 0.015$
| $0.351$
---
$\pm 0.773$
| $0.166$
---
$\pm 0.013$
| Timesteps to
---
Convergence
(Mean $\pm$ Std. Deviation)
| $178.39$
---
$\pm 46.75$
| $135.45$
---
$\pm 29.68$
| $153.15$
---
$\pm 41.34$
| $117.38$
---
$\pm 27.32$

TABLE II: Table comparing the RMSE loss, $d_{W}$ and Timesteps required to
converge for (i)RMSE-CLURSNN and (ii)Wass-CLURSNN using all the neurons in the
recurrent layer vs. sampling neurons with highest $\mathcal{C}_{b}$. The
experiment was repeated for the 4 modes of operation of L63 system, and the
mean and standard deviation is calculated from these observations.

Node Sampling:  We evaluate the model’s performance using all the neurons in
the recurrent layer compared to the case when we are sampling neurons based on
$\mathcal{C}_{b}$. As discussed earlier, the nodes are sampled according to
the ones with the highest $\mathcal{C}_{b}$ [26]. We calculate these two types
of predictions’ average RMSE and average $d_{W}$. We calculate the mean and
standard deviations over the four modes of operation of the L63 system. The
results are summarized in Table II. We observe that both the models perform
better and converge faster, with lesser standard deviation when using the
sampled neurons than all the neurons in the recurrent layer. We may interpret
this result as follows: since $\mathcal{C}_{b}$ captures a node’s role in
allowing information to pass from one part of the network to the other, it
behaves like the bottleneck layer in an autoencoder. Hence sampling neurons
with the highest $\mathcal{C}_{b}$, we can better capture the dynamical
system’s latent dynamics, leading to better performance.

Performance on Real-world Datasets:  We compared the performances of the
different models on some real-world datasets like YReal, and DJIA (as
described above). The results are tabulated in Table III. The proposed methods
RMSE-CLURSNN and Wass-CLURSNN outperform the other supervised DNN models.

TABLE III: Table comparing the RMSE performance of different models on real-
world datasets

Models Dataset | YReal | DJIA
---|---|---
DNN (Supervised) | SR-LSTM | 13816.99 | 1.76
|  | RLSTM
---
(online DNN)
10104.24 | 1.85
SNN (Unsupervised) |  | RMSE-CLURSNN
---
(Ours)
8938.25 | 1.41
|  | Wass-CLURSN
---
(Ours)
8875.73 | 1.31

##  V Conclusions

This paper proposes a novel methodology for online prediction of the time
series from an evolving dynamical system using a completely unsupervised bio-
plausible recurrent spiking neural network model trained with spike-timing-
dependent plasticity. We leveraged recent works on topological data analysis
and embedding theory to develop a new methodology for continuous
reconstruction of the underlying dynamical system based on membrane potential
measurements of nodes in the RSNN with the highest betweenness centrality,
which behaves similarly to the bottleneck layer in an autoencoder. We proposed
the RMSE-CLURSNN and the Wass-CLRSNN model, where the prediction error is
calculated using the RMSE error and the Wasserstein distance between the
persistent homologies of the predicted and observed time series $d_{W}$,
respectively. We compared this new methodology’s performance and convergence
behavior with the online DNN model. We observed that the proposed method
outperforms the standard DNN model in terms of prediction performance and
convergence for an evolving dynamical system with few data points for
training. However, with repeated iterations, we see that DNN-based models
converge faster and show better performance. Thus, these DNN models need a
large amount of data to learn. RSNN-based models, on the other hand, can adapt
very quickly to changing environments.

## Acknowledgement

This work is supported by the Army Research Office and was accomplished under
Grant Number W911NF-19-1-0447. The views and conclusions contained in this
document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the Army
Research Office or the U.S. Government.

## References

  * [1] Daniel Leite, Igor Škrjanc, and Fernando Gomide.  An overview on evolving systems and learning from stream data.  Evolving systems, 11:181–198, 2020.
  * [2] Zhuoqing Chang, Shubo Liu, Xingxing Xiong, Zhaohui Cai, and Guoqing Tu.  A survey of recent advances in edge-computing-powered artificial intelligence of things.  IEEE Internet of Things Journal, 8(18):13849–13875, 2021.
  * [3] Laizhong Cui, Shu Yang, Fei Chen, Zhong Ming, Nan Lu, and Jing Qin.  A survey on application of machine learning for internet of things.  International Journal of Machine Learning and Cybernetics, 9:1399–1417, 2018.
  * [4] Fotios Petropoulos, Daniele Apiletti, Vassilios Assimakopoulos, Mohamed Zied Babai, Devon K Barrow, Souhaib Ben Taieb, Christoph Bergmeir, Ricardo J Bessa, Jakub Bijak, John E Boylan, et al.  Forecasting: theory and practice.  International Journal of Forecasting, 2022.
  * [5] Aadyot Bhatnagar, Paul Kassianik, Chenghao Liu, Tian Lan, Wenzhuo Yang, Rowan Cassius, Doyen Sahoo, Devansh Arpit, Sri Subramanian, Gerald Woo, et al.  Merlion: A machine learning library for time series.  arXiv preprint arXiv:2109.09265, 2021.
  * [6] Shiyang Li, Xiaoyong Jin, Yao Xuan, Xiyou Zhou, Wenhu Chen, Yu-Xiang Wang, and Xifeng Yan.  Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting.  Advances in neural information processing systems, 32, 2019.
  * [7] R Rossi Pool and Germán Mato.  Spike-timing-dependent plasticity and reliability optimization: the role of neuron dynamics.  Neural computation, 23(7):1768–1789, 2011.
  * [8] Biswadeep Chakraborty and Saibal Mukhopadhyay.  Characterization of generalizability of spike timing dependent plasticity trained spiking neural networks.  Frontiers in Neuroscience, 15:695357, 2021.
  * [9] Biswadeep Chakraborty and Saibal Mukhopadhyay.  Heterogeneous recurrent spiking neural network for spatio-temporal classification.  arXiv preprint arXiv:2211.04297, 2022.
  * [10] Biswadeep Chakraborty and Saibal Mukhopadhyay.  Heterogeneous neuronal and synaptic dynamics for spike-efficient unsupervised learning: Theory and design principles.  In Submitted to The Eleventh International Conference on Learning Representations, 2023.  Accepted.
  * [11] Huanfei Ma, Siyang Leng, Kazuyuki Aihara, Wei Lin, and Luonan Chen.  Randomly distributed embedding making short-term high-dimensional data predictable.  Proceedings of the National Academy of Sciences, 115(43):E9994–E10002, 2018.
  * [12] Zhihan Yue, Yujing Wang, Juanyong Duan, Tianmeng Yang, Congrui Huang, Yunhai Tong, and Bixiong Xu.  Ts2vec: Towards universal representation of time series.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 8980–8987, 2022.
  * [13] João Gama, Indrė Žliobaitė, Albert Bifet, Mykola Pechenizkiy, and Abdelhamid Bouchachia.  A survey on concept drift adaptation.  ACM computing surveys (CSUR), 46(4):1–37, 2014.
  * [14] Oren Anava, Elad Hazan, Shie Mannor, and Ohad Shamir.  Online learning for time series prediction.  In Conference on learning theory, pages 172–184. PMLR, 2013.
  * [15] Chenghao Liu, Steven CH Hoi, Peilin Zhao, and Jianling Sun.  Online arima algorithms for time series prediction.  In Thirtieth AAAI conference on artificial intelligence, 2016.
  * [16] Doyen Sahoo, Quang Pham, Jing Lu, and Steven CH Hoi.  Online deep learning: Learning deep neural networks on the fly.  arXiv preprint arXiv:1711.03705, 2017.
  * [17] Rahaf Aljundi, Klaas Kelchtermans, and Tinne Tuytelaars.  Task-free continual learning.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11254–11263, 2019.
  * [18] Daehyun Kim, Biswadeep Chakraborty, Xueyuan She, Edward Lee, Beomseok Kang, and Saibal Mukhopadhyay.  Moneta: A processing-in-memory-based hardware platform for the hybrid convolutional spiking neural network with online learning.  Frontiers in Neuroscience, 16, 2022.
  * [19] Xueyuan She, Saurabh Dash, and Saibal Mukhopadhyay.  Sequence approximation using feedforward spiking neural network for spatiotemporal learning: Theory and optimization methods.  In International Conference on Learning Representations, 2021.
  * [20] Biswadeep Chakraborty, Xueyuan She, and Saibal Mukhopadhyay.  A fully spiking hybrid neural network for energy-efficient object detection.  IEEE Transactions on Image Processing, 30:9014–9029, 2021.
  * [21] David Reid, Abir Jaafar Hussain, and Hissam Tawfik.  Financial time series prediction using spiking neural networks.  PloS one, 9(8):e103656, 2014.
  * [22] Dan-Andrei Margin, Iustin-Alexandru Ivanciu, and Virgil Dobrota.  Deep reservoir computing using echo state networks and liquid state machine.  In 2022 IEEE International Black Sea Conference on Communications and Networking (BlackSeaCom), pages 208–213. IEEE, 2022.
  * [23] Peter Bubenik.  The persistence landscape and some of its properties.  In Topological Data Analysis, pages 97–117. Springer, 2020.
  * [24] Marian Gidea and Yuri Katz.  Topological data analysis of financial time series: Landscapes of crashes.  Physica A: Statistical Mechanics and its Applications, 491:820–834, 2018.
  * [25] Edward N Lorenz.  Deterministic nonperiodic flow.  Journal of atmospheric sciences, 20(2):130–141, 1963.
  * [26] Hamidreza Mahyar, Rouzbeh Hasheminezhad, Elahe Ghalebi, Ali Nazemian, Radu Grosu, Ali Movaghar, and Hamid R Rabiee.  Compressive sensing of high betweenness centrality nodes in networks.  Physica A: Statistical Mechanics and its Applications, 497:166–184, 2018.
  * [27] Balint Petro, Nikola Kasabov, and Rita M Kiss.  Selection and optimization of temporal spike encoding methods for spiking neural networks.  IEEE transactions on neural networks and learning systems, 31(2):358–370, 2019.
  * [28] Linton C Freeman.  A set of measures of centrality based on betweenness.  Sociometry, pages 35–41, 1977.
  * [29] N Laptev, S Amizadeh, and Y Billawala.  S5-a labeled anomaly detection dataset, version 1.0 (16m), 2015.
  * [30] Ronald J Williams and David Zipser.  A learning algorithm for continually running fully recurrent neural networks.  Neural computation, 1(2):270–280, 1989.


Illumina Sequencing Artifacts Revealed by Connectivity Analysis of Metagenomic
Datasets

Adina Chuang Howe1,2, Jason Pell3, Rosangela Canino-Koning3, Rachel
Mackelprang4, Susannah Tringe4, Janet Jansson4,5 , James M. Tiedje1,2, and C.
Titus Brown1,3∗

1 Microbiology and Molecular Genetics, Michigan State University, East
Lansing, MI, USA
2 Plant, Soil, and Microbial Sciences, Michigan State University, East
Lansing, MI, USA
3 Computer Science and Engineering, Michigan State University, East Lansing,
MI, USA
4 Department of Energy (DOE) Joint Genome Institute, Walnut Creek, CA, USA
5 Lawrence Berkeley National Laboratory, Genomics Division, Berkeley, CA, USA
$\ast$ E-mail<EMAIL_ADDRESS>
## Abstract

Sequencing errors and biases in metagenomic datasets affect coverage-based
assemblies and are often ignored during analysis. Here, we analyze read
connectivity in metagenomes and identify the presence of problematic and
likely a-biological connectivity within metagenome assembly graphs.
Specifically, we identify highly connected sequences which join a large
proportion of reads within each real metagenome. These sequences show
position-specific bias in shotgun reads, suggestive of sequencing artifacts,
and are only minimally incorporated into contigs by assembly. The removal of
these sequences prior to assembly results in similar assembly content for most
metagenomes and enables the use of graph partitioning to decrease assembly
memory and time requirements.

## Introduction

With the rapid decrease in the costs of sequencing, we can now achieve the
sequencing depth necessary to study microbes from even the most complex
environments [1, 2]. Deep metagenomic sequencing efforts in permafrost soil,
human gut, cow rumen, and surface water have provided insights into the
genetic and biochemical diversity of environmental microbial populations [1,
3, 2] and their involvement in responding to environmental changes [4]. These
metagenomic studies have all leveraged _de novo_ metagenomic assembly of short
reads for functional and phylogenetic analyses. _De novo_ assembly is an
advantageous approach to sequence analysis as it reduces the dataset size by
collapsing the more numerous short reads into fewer contigs and enables
improved annotation-based approaches by providing longer sequences [5, 6].
Furthermore, it does not rely on the a priori availability of reference
genomes to enable identification of gene content or operon structure [1, 3].

Although _de novo_ metagenomic assembly is a promising approach for
metagenomic sequence analysis, it is complicated by the variable coverage of
sequencing reads from mixed populations in the environment and their
associated sequencing errors and biases [7, 8]. Several metagenome-specific
assemblers have been developed to deal with variable coverage communities,
including Meta-IDBA [9], MetaVelvet [10], and SOAPdenovo [11]. These
assemblers rely on analysis of local sequencing coverage to help build
assemblies and thus are sensitive to the effects of sequencing errors and
biases on coverage estimations of the underlying dataset. The effects of
sequencing errors on _de novo_ assembly has been demonstrated in simulated
metagenomes [12, 7, 8] and isolate genomes [13, 14], but these datasets do not
necessarily represent real metagenomic data. Specifically, these models
exclude the presence of known non-biological sequencing biases which hinder
assembly approaches [15, 16, 17].

In this study, we examine metagenomic datasets for the presence of artificial
sequencing biases that affect assembly graph structure, extending previous
work to large and complex datasets produced from the Illumina platform. We
characterize sequence connectivity in an assembly graph, identifying potential
sequencing biases in regions where numerous reads are connected together.
Within metagenomic datasets, we find that there exist highly connected
sequences which partially originate from sequencing artifacts. Moreover, these
sequences limit approaches to divide or partition large datasets for further
analysis, and may introduce artifacts into assemblies. Here, we identify and
characterize these highly connected sequences and examine the effects of
removing these sequences on downstream assemblies.

## Results

### Connectivity analysis of metagenome datasets

#### Presence of a single, highly connected lump in all datasets

We selected datasets from three medium to high diversity metagenomes from the
human gut [2], cow rumen [1], and agricultural soil (SRX099904 and SRX099905)
(Table 1). To evaluate the effects of sequencing coverage, we included two
subsets of the 520 million read soil metagenome containing 50 and 100 million
reads. We also included a previously published error-free simulated metagenome
based on a mixture of 112 reference genomes [8].

We evaluated read connectivity by partitioning reads into disconnected
components with a de Bruijn graph representation [18]. This approach
guarantees that reads in different partitions do not connect to each other and
permits the separate assembly and analysis of each partition. For each
metagenome, regardless of origin, we found a single dominant, highly connected
set of sequencing reads which we henceforth refer to as the “lump” of the
dataset (Table 1). This lump contained the largest subset of connected
sequencing reads and varied in size among the datasets, ranging from 5% of
total reads in the simulated metagenome to 75% of total reads in the human gut
metagenome. For the soil datasets, as sequencing coverage (e.g., the fraction
of reads mapped to an assembly) increased from 1.4 to 4.7 to 5.6%, the lump
size increased more dramatically from 7 to 15 to 35% of all reads, indicating
increasingly larger connectivity between sequences with more sequencing.

#### Characterizing connectivity in the dominant partition

We characterized the connectivity of sequences within each lump by estimating
the average local graph density from each k-mer (k=32 unless otherwise stated)
in the assembly graph (see Methods). Here, local graph density is a
measurement of total connected reads within a fixed radius. Sequences in the
identified metagenomic lumps were characterized by very high local graph
densities: between 22 to 50% of the total nodes in metagenomic lump assembly
graphs had average graph densities greater than 20 (Table 1). This indicates
that these nodes were in very nonlinear portions of the assembly graph and had
high connectivity. In comparison, 17% of the total nodes in the simulated lump
had an average local graph density greater than 20, and fewer than 2% of the
nodes in the entire simulated data (all partitions) set had an average graph
density higher than 20.

We next assessed the extent to which graph density varied by position along
the sequencing reads. The degree of position-specific variation of graph
densities was estimated by calculating the average local graph density within
ten steps of every k-mer by position in each read. In all environmental
metagenomic reads, we observed variation in graph density at the 3’-end region
of reads (Fig. 1). In soil metagenomes, we observed the most dramatic
variation with local graph density increasing in sequences located at the
3’-end of the reads. Notably, this trend was not present in the simulated
dataset.

Next, we performed an exhaustive traversal of the assembly graph and
identified the specific sequences within dense regions of the assembly graph
which consistently contributed to high connectivity. We observed that this
subset of sequences was also found to exhibit position-specific variation
within sequencing reads, with the exception of these sequences in the
simulated dataset (Fig. 2, solid lines). As with local density trends,
position-specific trends in the location of these sequences also varied
between metagenomes. As sequencing coverage increased among metagenomes, the
amount of 3’-end variation appeared to decrease (e.g., the soils) or increase
(e.g., rumen and human gut).

### Effects of removing highly connected sequences on assembly

#### Removal of highly connected sequences enables graph partitioning of
metagenome

Since these highly connected sequences exhibited position-specific variation
indicative of sequences of non-biological origin, we removed them and assessed
the effect of their removal on assembly (see Methods). We found that by
removing these k-mers, we could effectively break apart metagenomic lumps, and
the resulting largest partition of connected reads in each metagenome was
reduced to less than 7% of the total reads in the lump. Partitioning also had
the effect of significantly decreasing assembly time and memory usage [18].

#### Removing highly connected sequences resulted in minimal losses of
reference genes

We explored the extent to which the identified highly connected sequences
impacted assembly by first evaluating the effects of the removal of these
sequences from the simulated lump. The assembly of the reads in the original,
unfiltered simulated lump and that of the reads remaining after removing
highly connected sequences (the filtered assembly) were compared for three
assemblers: Velvet [19], Meta-IDBA [9], and SOAPdenovo [11]. Based on the
total assembly length of contigs greater than 300 bp, filtered assemblies of
the simulated metagenome resulted in a loss of between 4 - 16% of total
assembly length (Table 2). In general, the filtered assemblies contained fewer
total contigs than unfiltered assemblies, while the maximum contig size
increased in the Velvet assembly but decreased in the Meta-IDBA and SOAPdenovo
assemblies. Direct comparisons of the unfiltered and filtered simulated
metagenome assemblies found that the filtered assemblies comprised on average
89% of the unfiltered assemblies, and the unfiltered assemblies contained
nearly all (97%) of the filtered assembled sequences. Despite the removal of
over 3% of the total unique 32-mers in the simulated metagenome, the resulting
filtered assemblies lost only 3-15% of annotated original reference genes
(Table 3).

We next evaluated the effects of removing highly connected sequences in real
datasets. Similar to the simulated assemblies, the removal of highly connected
sequences for all metagenomes and assemblers resulted in a decrease of total
number of contigs and assembly length (Table 2). In general, filtered
assemblies were largely contained within unfiltered assemblies and comprised
51-87% of the unfiltered assembly. The observed changes in metagenomic
assemblies were difficult to evaluate as no reference genomes exist, and a
decrease in assembly length may actually be beneficial if it eliminates
contigs that incorporate sequencing artifacts. To aid in this evaluation, we
used the previously published set of rumen draft genomes from _de novo_
assembly efforts of high abundance sequences in the rumen metagenome [1].
Overall, we found that removal of highly connected sequences from the rumen
dataset resulted in 9-13% loss of sequences present in draft reference genomes
(Table 3).

#### Unfiltered assemblies contained only a small fraction of highly connected
sequences

To further study the effects of highly connected sequences, we examined their
incorporation into unfiltered assemblies. Except in the human gut sample,
fewer than 2% of highly connected sequences were incorporated by any assembler
(Table 4). Each assembled contig was divided into percentile bins and examined
for the presence of the previously identified highly connected sequences. We
found that contigs, especially in assemblies from Velvet and Meta-IDBA,
incorporated a larger fraction of these sequences at their ends relative to
other positions (Fig 3). The SOAPdenovo assembler incorporated fewer of the
highly connected sequences into its assembled contigs; in the simulated data
set, none of these sequences were assembled, and in the small soil data set
only 41 were assembled. For the human gut metagenome assemblies, millions of
the highly connected sequences were incorporated into assembled contigs,
comprising nearly 4% of all assembled sequences on Velvet contig ends (Fig 3).

#### Identifying origins of highly connected sequences in known reference
databases

For the simulated metagenome, we could identify the source of highly connected
k-mers using available reference genomes. Reference genes with multiple
perfect alignments to highly connected k-mers present in the dataset a minimum
of 50 times were identified (Table 5). Many of these sequences were from well-
conserved housekeeping genes involved in protein synthesis, cell transport,
and signaling. To determine possible biological sources of highly connected
sequences within real metagenomes, we compared the sequences shared between
the soil, rumen, and human gut metagenomes (a total of 241 million 32-mers).
Among these 7,586 shared sequences, we identified the closest reference
protein from the NCBI-nr database requiring complete sequence identity. Only
1,018 sequences (13%) matched existing reference proteins, and many of the
annotated sequences matched to genes conserved across multiple genomes. The
most abundant proteins conserved in greater than 3 genomes are shown in (Table
6), and largely encode for genes involved in protein biosynthesis, DNA
metabolism, and biochemical cofactors.

One potential cause of artificial high connectivity within metagenomes is the
presence of high abundance subsequences. Thus, we identified the subset of
highly connected k-mers which were also present with an abundance of greater
than 50 within each metagenome and their location in sequencing reads (Fig 2,
dotted lines). These high abundance k-mers comprised a very small proportion
of the identified highly connected sequences, less than 1% in the soils, 1.5%
in the rumen, and 6.4% in the human gut metagenomes, but the position-specific
variation of these sequences was very similar to the variation in the larger
set of highly connected k-mers.

We attempted to identify patterns within the sequences causing position
specific variation by examining the abundance distribution of 5-mers within
the highly abundant/highly connected 32-mers. There were significantly fewer
5-mers in the simulated sequences compared to those in metagenomes: 336 total
5-mers in the simulated data and from 425,572 to 221,085,228 total 5-mers in
the small soil and human gut datasets, respectively. In the simulated dataset,
the top ten most abundant unique k-mers made up 75% of the total 5-mers; in
contrast, in the metagenomes, k-mers were more evenly distributed: the top ten
most abundant 5-mers comprised less than 10% of the total 5-mers. The
cumulative abundance distribution of the ranked 5-mers shown in Fig. 4 shows
this even distribution in all of the real metagenomes. This suggests that
there is no single, easily-identifiable set of sequences at the root of the
highly connected component observed in real metagenomes.

## Discussion

### Sequencing artifacts are present in real metagenomes

Through assessing the connectivity of reads in several metagenomes, we
identified a disproportionately large subset of reads connected together
within an assembly graph, which we refer to as the “lump.” The total number of
reads in metagenomic lumps (7-75% of reads) was significantly larger than that
of simulated dataset (5% of reads) (Table 1). In the simulated data, this
component consists of reads connected by sequences conserved between multiple
genomes (Table 5). The larger size of this component within the soil, rumen,
and human gut metagenomes suggests that anomalous, non-biological connectivity
may be present within these lumps. Moreover, in the soil metagenomes, we
observed that for a 5% increase in sequencing coverage of the assembled
contigs, the amount of connectivity nearly doubled. While sequencing coverage
of assembled contigs increased slightly from 4.7 to 5.6% in the medium and
large soil metagenomes, the number of reads contained in the lump grew
significantly from 15 million to 182 million. Given the very high diversity
and very low coverage of these soil samples, the magnitude of the observed
increases in connectivity cannot be due simply to increased coverage. This
suggests the presence of non-biological features that falsely connect reads.

The superlinear increase in connectivity exhibited in these data sets
indicates that a form of preferential attachment is occurring in the graphs
[20]. This graph-theoretic phenomenon describes the results of a process where
highly connected nodes in a graph preferentially acquire new edges –
colloquially known as “the rich get richer.” In assembly, any systematic bias
towards producing specific subsequences from shotgun sequencing would lead to
a tendency to connect otherwise unrelated graph components; such a bias could
be biological (e.g. repeat present in multiple genomes or other highly
conserved DNA sequences), or non-biological (e.g., inclusion of sequencing
primers in reads or even a low-frequency trend towards producing specific
subsequences [21, 22, 23]).

We believe a significant component of the high connectivity that we see is of
non-biological origin. Shotgun sequencing is a random process and consequently
any position-specific variation within sequencing reads is unexpected and
probably originates from bias in sample preparation or the sequencing process
[15, 24, 16]. For the metagenomes studied here, we used two approaches to
examine characteristics of connectivity correlated to specific positions
within sequencing reads. First, we measured the connectivity of sequences at
specific positions within reads by calculating local graph density. Next, we
identified the specific k-mers which were consistently present in highly dense
regions of the assembly graph and evaluated their location within sequencing
reads. When these approaches were applied to the simulated dataset, we
observed no position-specific trends when assessing either local graph density
(Fig 1) or highly connected k-mers (Fig 2) as is consistent with the lack of
sequencing errors and variation in this dataset. In all real metagenomes,
however, we identified position-specific trends in reads for measurements of
both local graph density and the location of highly connected sequences,
clearly indicating the presence of sequencing artifacts. Although present in
all metagenomes, the direction of the variation varied between soil, rumen,
and human gut datasets, especially for the position-specific presence of
identified highly connected sequences. It is likely that there is a larger
presence of indirectly preferentially attached reads which are connected to
high coverage sequences of biological origins in higher coverage datasets,
such as the rumen and human gut. This preferential attachment of such reads
would result in increasing the number of total reads and consequently the
decrease the total fraction of highly connected k-mers (Fig 2, y-axis). This
trend is observed in the decreasing fractions of highly connected sequences at
the 3’ end of reads as sequencing coverage increased in the small, medium, to
large soil metagenomes and in the soil, rumen, to human gut metagenomes (Fig
2).

### Highly connected sequences are of unknown non-biological origin

We attempted to identify biological characteristics of highly connected
sequences. Among the highly connected sequences in the simulated dataset and
those shared by all metagenomes, we identified only a small fraction (13% in
simulated and less than 7% in metagenomes) which matched reference genes
associated with core biological functions (Table 5 and 6). This suggests that
the remaining sequences are either not present in known reference genes (i.e.,
repetitive or conserved non-coding regions) or originate from non-biological
sources. This supports the removal of these sequences for typical assembly and
annotation pipelines, where assembly is often followed by the identification
of protein coding regions.

Speculating that many of the highly connected sequences originated from high
abundance reads, we examined the most abundant subsequences. We found that
these subsequences (k-mers present more than 50x in the data set) displayed
similar trends for position-specific variation compared to their respective
sets of highly connected subsequences (Fig 2), indicating that they contribute
significantly to position-specific variation. We attempted to identify
signatures in these abundant, highly connected sequences from the simulated
and metagenomic datasets by looking at shorter k-mer profiles. In the
simulated dataset, we found that the total number of unique 5-mers was
significantly lower than in metagenomes and that the most abundant of these
5-mers comprised the large majority of the total. This result is consistent
with the presence of conserved biological motifs in the simulated dataset
which would result in a small number of highly abundant sequences; it would
also be consistent with the inclusion of sequencing primers in the data, were
this a real data set.

In contrast, within real metagenomic data, we found that the 5-mers are evenly
distributed and exhibit no specific sequence properties (Fig 4), making them
difficult to identify and evaluate. Most importantly, we were unable to
identify any characteristics that would explain their origin. In addition, a
G-C content analysis of the highly connective k-mers did not reveal any
systematic differences between the highly connected k-mers and the background
k-mer distribution.

When we reviewed the literature on random and systematic sequencing errors in
Illumina sequencing, we found many different types of sequencing errors: PCR
amplification errors prior to and during cluster generation; random sequencing
errors e.g. from miscalls of bases; sequencing errors triggered by specific
sequence motifs [25]; adaptor contamination; and post-adaptor read through. Of
these errors, only random sequencing errors and adaptor contamination and
readthrough would be biased towards the 3’ end of the read. However, random
sequencing error does not contribute to aberrant de Bruijn graph connectivity
[18], while adaptor contamination and readthrough would yield a sharply biased
5-mer distribution. The observed artifactual sequences thus do not match any
known set of random or systematic errors in Illumina sequencing.

Our current working hypothesis is that a low rate of false connections are
created by a low-frequency tendency towards producing certain k-mers in the
Illumina base calling software, as signal intensities decline. We cannot
verify this without access to the Illumina software or source code.

### Highly connected sequences are difficult to assemble

Not all of the observed connectivity within real metagenomes is artificial,
and our approach cannot differentiate between sequencing artifacts and real
biological connectivity. Therefore, removing highly connected sequences could
remove real biological signal in addition to sequencing artifacts. However, we
suspected that assemblers would be unable to generate contigs from highly
connected graph regions, and so even the real sequences would be
underrepresented in the assembly.

Indeed, very few highly connected sequences with abundances greater than 50
were incorporated into contigs (Table 4). Moreover, those which were assembled
were often disproportionately placed at the ends of contigs (Fig 3),
demonstrating that they terminated contig assembly. Although this trend was
observed for all three assemblers, it was more prevalent in the Velvet and
Meta-IDBA assemblers, highlighting differences in assembler heuristics.

One initial concern upon discovering this false connectivity was that these
artifacts might nucleate false assemblies, e.g. as seen with the cow rumen
[1]. While we only examined their effects on contig assembly and ignored
scaffolding issues, misassembly is an especially significant concern for
environmental metagenomics, where experimental validation of assemblies is
virtually impossible. However, our observation that k-mers from highly
connected graph regions are generally not present in contigs suggests that
these particular artifacts do not create false assemblies at a high rate.

### Filtered assemblies retained most reference genes

The advantages of removing highly connected sequences must be balanced against
consequences to resulting assemblies. We compared several metagenome
assemblies before and after the removal of these sequences. In comparing the
simulated dataset’s assemblies, the removal of highly connected sequences
resulted in very little loss of annotated reference genes (less than 1% total)
and some loss of assembled contigs ($\sim$ 15% of the final assembly). For the
rumen metagenome, we performed a partial evaluation of the assemblies using
available draft reference genomes. Similar to the simulated assemblies, we
observed only a small loss (less than 3% total) of rumen reference genomes
assembled (Table 3). In general, for all metagenomes, we observed $\sim$ 25%
loss in assembly after removing highly connected sequences, much more than
observed in assemblies of reference genes and genomes in the simulated and
rumen datasets. Some of this loss could be beneficial, resulting from removal
of sequencing artifacts. In addition, trimming or removing reads from already
low-coverage data sets could result in the loss of contigs due to the length
cutoff we use; as observed above, most of these highly connected k-mers fall
at the ends of contigs, and removing them could shorten the contigs enough to
lose them from our assembly. It is also possible, of course, that our approach
removes sequences which can accurately be assembled, but we cannot evaluate
this in the absence of reference genomes.

### Filtered reads can be assembled more efficiently

Our original motivation for analyzing connectivity was to assess the
practicality of partitioning metagenomes for later assembly, which can
dramatically reduce the memory requirements for assembly [18]. Not only did
removal of highly connected sequences eliminate problematic sequences, but it
resulted in the dissolution of the largest component and allowed us to
partition the metagenomes.

We compared the combined assembly of the partitioned sets of filtered reads to
the original lump dataset, for several assemblers. For the partitioned reads,
we were able to assemble subsets of reads in parallel, resulting in
significantly reduced time and memory requirements for assembly (Table 2). In
the case of the largest soil metagenome (containing over 500 million reads),
we could not complete the Meta-IDBA assembly of the unfiltered reads in even
100 GB of memory, but after removing highly connected sequences and
partitioning, the assembly could be completed in less than 2 GB of memory.
Using partitioned sets of reads for all metagenomes, we were also able to
efficiently complete multiple k-mer length assemblies (demonstrated with
Velvet) and subsequently merge resulting assembled contigs. For unfiltered
datasets, this was either impossible (due to memory limitations) or
impractical (due to excessive processing time).

## Conclusion

In this study, we characterize the connectivity of sequences in several
metagenomes to detect and characterize a set of likely sequencing artifacts.
These artifacts are surprisingly abundant, comprising 5-13% of the total
unique sequence in real data sets. Moreover, they have a significant impact on
the overall graph connectivity of the data sets, leading to components
containing as many as 75% of the reads in the human gut data set.

Despite the prevalence and impact of these sequences, removal of the highly
connected k-mers at the heart of these components leads to assemblies that are
significantly but not catastrophically smaller than the original assemblies.
The original unfiltered assemblies contain the majority of the filtered
assemblies, while the filtered assemblies generally contain 70-94% of the
unfiltered assemblies. The variability in these statistics between the
different assemblers (Table 2) demonstrates that the assemblers have at least
as large an effect on the content of the assemblies as our filtering
procedure.

We cannot reach strong conclusions about the impact of these highly connected
sequences on the correctness of the assembled contigs. In particular, in the
absence of complex metagenomes that have been characterized by approaches
other than short-read shotgun approaches, we have no very high quality
positive control. However, we present evidence that these highly connected
sequences primarily affect the ends of contigs and thus are unlikely to cause
incorrect contig assembly with high frequency.

Our original motivation in exploring metagenome connectivity was to enable
partitioning, an approach that leads to substantially greater scalability of
the assembly procedure. In this respect, we were successful. By applying
partitioning to filtered metagenome data, we were able to reduce the maximum
memory requirements of assembly (including the filtering stage) to well below
48 GB of RAM in all cases. This enables the use of commodity “cloud” computing
for all of our samples [26]. The decreased computational requirements for
assembly also enabled ready evaluation of different assemblers and assembly
parameters; as metagenome datasets grow increasingly larger, this ability to
efficiently analyze datasets and evaluate multiple assemblies will be
increasingly important.

More generally, our results demonstrate that metagenome assembly is still at
an early stage of technology development, particularly for low-coverage data
sets. Different assemblers and different filtering techniques yield
substantially different assembly statistics on these data sets. Moreover,
there are likely to be additional sources of sequencing artifacts lurking
within large sequencing data sets, suggesting that more and better
computational filtering and validation approaches need to be developed as
environmental metagenomics moves forward. Evaluating the assembly graph
connectivity created by reads will be a useful approach in the future.

## Methods

### Metagenomic datasets

All datasets, with the exception of the agricultural soil metagenome,
originate from previously published datasets. Rumen-associated sequences
(Illumina) were randomly selected from the rumen metagenome (read length 36 -
125 bp) available at ftp://ftp.jgi-psf.org/pub/rnd2/Cow_Rumen [1]. Human-gut
associated sequences (Illumina) of samples MH0001 through MH0010 were obtained
from
ftp://public.genomics.org.cn/BGI/gutmeta/ Raw_Reads [2] (read length 44 bp).
The simulated high complexity, high coverage dataset was previously published
[8]. Soil metagenomes (read lengths 76-113 bp) are in the SRA (SRX099904 and
SRX099905). All reads used in this study, with the exception of those from the
simulated metagenome, were quality-trimmed for Illumina’s read segment quality
control indicator, where a quality score of 2 indicates that all subsequent
regions of the sequence should not be used. After quality-trimming, only reads
with lengths greater than 30 bp were retained. All quality trimmed datasets,
including the previously unpublished agricultural soil metagenome, are
available on a public Amazon EC2 snapshot (snap-ab88dfdb). The sequencing
coverage of each metagenome was estimated as the fraction of reads which could
be aligned to assembled contigs with lengths greater than 500 bp. For the
coverage estimates, an assembly of each metagenome was performed using Velvet
(v1.1.02) with the following parameters: K=33, exp cov=auto, cov cutoff=0, no
scaffolding. Reads were aligned to assembled contigs with Bowtie (v0.12.7),
allowing for a maximum of two mismatches.

### de Bruijn graph analysis and partitioning software

We used the probabilistic de Bruijn graph representation previously described
by [18] to store and partition the metagenome assembly graphs. The khmer and
screed software packages are required for the analysis, and the versions used
for this publication are available at https://github.com/ged-
lab/khmer/tree/2012-assembly-artifacts and https://github.com/ged-
lab/screed/tree/2012-assembly-artifacts.

For metagenomes in this study, we used 4 x 48e9 bit bloom filters (requiring
24 GB RAM) to store the assembly graphs. The data processing pipeline used for
this analysis is available for public use on the Amazon Web Services public
EBS snapshot snap-ab88dfdb: data-in-paper/lumps and method-
examples/0.partitioning-into-lump.

The local graph density was calculated as the number of k-mers within a
distance of N nodes divided by N. In this study, N was equal to 10. For the
largest metagenomes, the human gut and large soil datasets, local graph
density was calculated on a randomly chosen subset of reads because of
computational limitations.

To identify specific highly connected sequences within the lump assembly
graphs, graph traversal to a distance of 40 nodes was attempted from marked
waypoints. If more than 200 k-mers were found within this traversal were
identified (i.e. a graph density $>5$), all k-mers within this traversal were
marked. If the same k-mers were consistently identified in other graph
traversals, up to five times, the k-mer was flagged as a highly connected
sequence. Aligning these k-mers to original sequencing reads, we identified
the position-specific location of these k-mers. Data and examples of scripts
used for this analysis are available on the Amazon EC2 public snapshot: data-
in-paper/density-bias, data-in-paper/hc-kmer-bias, method-examples/1.density-
analysis, method-examples/2.identifying-hc-kmers, and method-examples/3.hc-
kmer-analysis.

We identified the sources of highly connected k-mers from the simulated
metagenome by aligning them against reference genes originating from the 112
source genomes using Bowtie (v0.12.7) requiring exact matches. Highly
connected k-mers shared between all the metagenomes were also aligned against
the NCBI non-redundant genome database (ftp://ftp.ncbi.nih.gov/blast/db,
March, 1, 2011) using blastn [27], requiring an exact match over the entire
k-mer.

We also identified the subset of highly connected 32-mers which were present
greater than 50 times within lumps. Data used for this analysis are available
on the Amazon EC2 public snapshot: data-in-paper/lumps/HC-kmers/HA-HC-kmers
and method-examples/4.abundant-hc-kmers. These high abundance, highly
connected sequences were aligned to sequencing reads to demonstrate position
specific variation as described above. We evaluated the existence of short
k-mer (k=5) motifs within high abundance, highly connected k-mers which did
not have an exact match to the NCBI non-redundant database. Each identified
32-mer was broken up into shorter 5-mers, and the frequency of each unique
5-mer was calculated. Next, each unique 5-mer was ranked based on its
abundance, from high to low, and the cumulative percentage of total 5-mers is
shown in the resulting rank-abundance plot (Fig 4).

###  _De novo_ metagenomic assembly

The lump within each dataset was assembled and referred to as the “unfiltered
assembly”. Additionally, highly connected sequences identified as described
above were were trimmed from sequencing reads and the remaining reads
partitioned and assembled, resulting in the “filtered assembly”. _De novo_
metagenomic assembly of reads was completed with Velvet (v1.1.02) with the
following parameters: velveth -short -shortPaired (if applicable to the
dataset) and velvetg -exp_cov auto -cov_cutoff 0 -scaffolding no [19]. For the
small and medium soil, rumen, and simulated datasets, Velvet assemblies were
performed at K=25-49, resulting contigs were dereplicated to remove contigs
with 99% similarity using CD-HIT (v 4.5.6, [28]), and final contigs were
merged with Minimus (Amos v3.1.0, [29]). For the largest soil and human gut
metagenomes, assemblies were performed at only K=33 due to the size of the
datasets and memory limitations. Additional assemblies were performed with
Meta-IDBA (v0.18) [9] : –mink 25 –maxk 50 –minCount 0 and with SOAPdenovo: -K
31 -p 8 max_rd_len=200 asm_flags=1 reverse_seq=0. After removal of highly
connected k-mers in metagenomic lumps, each filtered lump was partitioned into
separate disconnected subgraphs. Multiple subgraphs were grouped together such
that assembly could be performed in parallel on groups of sequences. Identical
assembly parameters and methods as described above were used for these
assemblies. Unfiltered and filtered assemblies were compared using the total
number of contigs, total assembly length, and maximum contig size. Additional,
the coverage of each assembly was calculated through estimating the average
base pair coverage of the BLAST alignment of each assembly to one another
(E-value less than 10-5) or, in the case of the simulated and rumen
assemblies, to reference genomes. The simulated and rumen reference genomes
were previously published in [1] and [8], respectively. Resulting assemblies
are available on the Amazon EC2 public snapshot: /data-in-paper/assembly*.

We examined incorporation and the location of the identified high abundant,
highly connecting k-mers within assembled contigs. Incorporation of these
sequences was evaluated by dividing assembled contigs into words of 32 bp
length and identifying exact matches between sequences and contig fragments.
The location of these k-mers within assembled unfiltered contigs was examined
by dividing each contig into 100 equally-sized regions. The fraction of highly
connecting k-mers which aligned exactly to each region was calculated for each
metagenome. Data and examples of scripts used for this analysis are available
on the Amazon EC2 public snapshot: method-examples/5.hc-kmer-contigs/.

## Acknowledgements

This project was supported by Agriculture and Food Research Initiative
Competitive Grant no. 2010-65205-20361 from the United States Department of
Agriculture, National Institute of Food and Agriculture and National Science
Foundation IOS-0923812, both to C.T.B. A.H. was supported by NSF Postdoctoral
Fellowship Award #0905961 and the Great Lakes Bioenergy Research Center
(Department of Energy BER DE-FC02-07ER64494). The work conducted by the U.S.
Department of Energy Joint Genome Institute is supported by the Office of
Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.
We thank Nick Loman and Lex Nederbragt for their helpful comments on the
paper.

## References

  * 1. Hess M, Sczyrba A, Egan R, Kim TW, Chokhawala H, et al. (2011) Metagenomic discovery of biomass-degrading genes and genomes from cow rumen.  Science 331: 463–7.
  * 2. Qin J, Li R, Raes J, Arumugam M, Burgdorf KS, et al. (2010) A human gut microbial gene catalogue established by metagenomic sequencing.  Nature 464: 59–65.
  * 3. Iverson V, Morris RM, Frazar CD, Berthiaume CT, Morales RL, et al. (2012) Untangling genomes from metagenomes: revealing an uncultured class of marine euryarchaeota.  Science 335: 587–90.
  * 4. Mackelprang R, Waldrop MP, DeAngelis KM, David MM, Chavarria KL, et al. (2011) Metagenomic analysis of a permafrost microbial community reveals a rapid response to thaw.  Nature 480: 368–71.
  * 5. Miller JR, Koren S, Sutton G (2010) Assembly algorithms for next-generation sequencing data.  Genomics 95: 315–27.
  * 6. Pop M (2009) Genome assembly reborn: recent computational challenges.  Briefings in bioinformatics 10: 354–66.
  * 7. Mende DR, Waller AS, Sunagawa S, Järvelin AI, Chan MM, et al. (2012) Assessment of metagenomic assembly using simulated next generation sequencing data.  PLoS ONE 7: e31386.
  * 8. Pignatelli M, Moya A (2011) Evaluating the fidelity of de novo short read metagenomic assembly using simulated data.  PLoS ONE 6: e19984.
  * 9. Peng Y, Leung HCM, Yiu SM, Chin FYL (2011) Meta-idba: a de novo assembler for metagenomic data.  Bioinformatics 27: i94–101.
  * 10. Namiki T, Hachiya T, Tanaka H, Sakakibara Y (2012) MetaVelvet: an extension of Velvet assembler to de novo metagenome assembly from short sequence reads.  Nucleic Acids Research .
  * 11. Li R, Zhu H, Ruan J, Qian W, Fang X, et al. (2010) De novo assembly of human genomes with massively parallel short read sequencing.  Genome Res 20: 265–72.
  * 12. Mavromatis K, Ivanova N, Barry K, Shapiro H, Goltsman E, et al. (2006) Use of simulated data sets to evaluate the fidelity of metagenomic processing methods.  Nature Methods 4: 495–500.
  * 13. Morgan J, Darling A, Eisen J (2010) Metagenomic sequencing of an in vitro-simulated microbial community.  PLoS ONE 5: e10209.
  * 14. Chitsaz H, Yee-Greenbaum JL, Tesler G, Lombardo MJ, Dupont CL, et al. (2011) Efficient de novo assembly of single-cell bacterial genomes from short-read data sets.  Nature Biotechnology 29: 915–921.
  * 15. Gomez-Alvarez V, Teal TK, Schmidt TM (2009) Systematic artifacts in metagenomes from complex microbial communities.  The ISME Journal 3: 1314–7.
  * 16. Keegan KP, Trimble WL, Wilkening J, Wilke A, Harrison T, et al. (2012) A platform-independent method for detecting errors in metagenomic sequencing data: Drisee.  PLoS Comput Biol 8: e1002541.
  * 17. Niu B, Fu L, Sun S, Li W (2010) Artificial and natural duplicates in pyrosequencing reads of metagenomic data.  BMC Bioinformatics 11: 187.
  * 18. Pell J, Hintze A, Canino-Koning R, Howe A, Tiedje JM, et al. (2012) Scaling metagenome sequence assembly with probabilistic de Bruijn graphs.  Proceedings of the National Academy of Sciences of the United States of America 109: 13272–13277.
  * 19. Zerbino DR, Birney E (2008) Velvet: algorithms for de novo short read assembly using de bruijn graphs.  Genome Res 18: 821–9.
  * 20. Barabási A, Albert R (1999) Emergence of scaling in random networks.  Science 286: 509.
  * 21. Hansen KD, Brenner SE, Dudoit S (2010) Biases in Illumina transcriptome sequencing caused by random hexamer priming.  Nucleic Acids Research 38: e131–e131.
  * 22. Minoche AE, Dohm JC, Himmelbauer H (2011) Evaluation of genomic high-throughput sequencing data generated on Illumina HiSeq and Genome Analyzer systems.  Genome Biology 12: R112.
  * 23. Dohm JC, Lottaz C, Borodina T, Himmelbauer H (2008) Substantial biases in ultra-short read data sets from high-throughput DNA sequencing.  Nucleic Acids Research 36: e105–e105.
  * 24. Haas BJ, Gevers D, Earl AM, Feldgarden M, Ward DV, et al. (2011) Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons.  Genome Research 21: 494–504.
  * 25. Meacham F, Boffelli D, Dhahbi J, Martin DI, Singer M, et al. (2011) Assessment of metagenomic assembly using simulated next generation sequencing data.  BMC Bioinformatics 12.
  * 26. Angiuoli SV, White JR, Matalka M, White O, Fricke WF (2011) Resources and Costs for Microbial Sequence Analysis Evaluated Using Virtual Machines and Cloud Computing.  PLoS ONE 6: e26624.
  * 27. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ (1990) Basic local alignment search tool.  J Mol Biol 215: 403–10.
  * 28. Li W, Jaroszewski L, Godzik A (2001) Clustering of highly homologous sequences to reduce the size of large protein databases.  Bioinformatics 17: 282–3.
  * 29. Sommer DD, Delcher AL, Salzberg SL, Pop M (2007) Minimus: a fast, lightweight genome assembler.  Bmc Bioinformatics 8: 64.

Table 1: The original size and proportion of highly connective 32-mers in the largest subset of partitioned reads (“lump”) in several medium to high complexity metagenomes. Read coverage was estimated with the number of aligned sequencing reads to Velvet-assembled contigs (K=33). The dominant lump, or largest component of each metagenome assembly graph, was found to contain highly connecting (HC) k-mers responsible for high local graph density. High density nodes refer to nodes with graph density greater than 20. | Sm Soil | Med Soil | Large Soil | Rumen | Human Gut | Sim
---|---|---|---|---|---|---
Total Reads (millions) | 50.0 | 100.0 | 520.3 | 50.0 | 350.0 | 9.2
Mapped to assembly (percent) | 1.4 | 4.7 | 5.6 | 10.3 | 3.5 | 14.8
Reads in Lump (millions) | 3.0 | 15.0 | 182.2 | 10.3 | 263 | 0.5
Lump Fraction (%) | 7% | 15% | 35% | 21% | 75% | 5%
HC 32-mers (millions) | 6.4 | 33.3 | 230.4 | 25.4 | 136.6 | 0.4
Total 32-mers (million) | 84.9 | 326.5 | 2,198.1 | 201.5 | 860.6 | 11.6
Fraction of HC 32-mers (%) | 8% | 10% | 10% | 13% | 16% | 3%
High Density Nodes (%) | 50% | 37% | 40% | 22% | 28% | 17%
Table 2: Total number of contigs, assembly length, and maximum contig size was
estimated for metagenomic datasets with multiple assemblers, as well as memory
and time requirements of unfiltered read assembly (UF). Filtered reads (F)
were processed in 24 GB of memory, and after filtering required less than 2 GB
of memory to assemble. Velvet assemblies of the unfiltered human gut and large
soil datasets (marked as *) could only be completed with K=33 due to
computational limitations. The Meta-IDBA assembly of the large soil metagenome
could not be completed in less than 100 GB.

| UF Assembly | F Assembly | UF Requirements
---|---|---|---
| (contigs / length / max size) | (contigs / length / max size) | Memory (GB)/Time (h)
_Velvet_ |  |  |
Small Soil | 25,470 / 16,269,879 / 118,753 | 17,636 / 10,578,908 / 13,246 | 5 / 4
Medium Soil | 113,613 / 81,660,678 / 57,856 | 79,654 / 54,424,264 / 23,663 | 18 / 21
Large Soil | 554,825 / 306,899,884 / 41,217 | 290,018 / 159,960,062 / 41,423 | 33 / 12*
Rumen | 92,044 / 74,813,072 / 182,003 | 72,705 / 49,518,627 / 34,683 | 11 / 14
Human Gut | 543,331 / 234,686,983 / 85,596 | 203,299 / 181,934,800 / 145,740 | 76 / 8*
Simulated | 11,204 / 6,506,248 / 5,151 | 9,859 / 5,463,067 / 6,605 | <1 / <1

_MetaIDBA_ |  |  |
---|---|---|---
Small Soil | 15,739 / 9,133,564 / 37,738 | 12,513 / 7,012,036 / 17,048 | <1 / <1
Medium Soil | 76,269 / 45,844,975 / 37,738 | 52,978 / 30,040,031 / 18,882 | 2 / 2
Large Soil | 395,122 / 228,857,098 / 37,738 | N/A | >116 / incomplete
Rumen | 60,330 / 47,984,619 / 54,407 | 48,940 / 33,276,502 / 22,083 | 12 / 3
Human Gut | 173,432 / 211,067,996 / 106,503 | 132,614 / 142,139,101 / 85,539 | 58 / 15
Simulated | 8,707 / 4,698,575 / 5,113 | 7,726 / 4,078,947 / 3,845 | <1 / <1

_SOAPdenovo_ |  |  |
---|---|---|---
Small Soil | 14,275 / 7,100,052 / 37,720 | 12,801 / 6,343,110 / 13,246 | 3 / <1
Medium Soil | 66,640 / 33,321,411 / 28,695 | 56,023 / 27,880,293 / 15,721 | 10 / <1
Large Soil | 412,059 / 215,614,765 / 32,514 | 334,319 / 171,718,154 / 41,423 | 48 / 11
Rumen | 62,896 / 40,792,029 / 22,875 | 55,975 / 34,540,861 / 19,044 | 5 / <1
Human Gut | 190,963 / 171,502,574 / 57,803 | 161,795 / 139,686,630 / 56,034 | 35 / 5
Simulated | 6,322 / 2,940,509 / 3,786 | 6,029 / 2,821,631 / 3,764 | <1 / <1

Table 3: Comparison of unfiltered (UF) and filtered (F) assemblies of various
metagenome lumps using Velvet, SOAPdenovo, and Meta-IDBA assemblers.
Assemblies were aligned to each other, and coverage was estimated (columns
1-2). Simulated and rumen assemblies were aligned to available reference
genomes (RG) (columns 3-4).

Velvet Assembler
Cov. of UF by F Cov. of F by UF Cov. of RG by UF Cov. of RG by F Simulated
85.4% 99.4% 5.4% 4.6% Small Soil 74.7% 98.8% - - Medium Soil 75.6% 98.4% - -
Large Soil 50.9% 86.6% - - Rumen 75.9% 98.8% 17.5% 14.8% Human Gut 80.0% 89.1%
- - Meta-IDBA Assembler
Cov. of UF by F Cov. of F by UF Cov. of RG by UF Cov. of RF by F Simulated
87.4% 94.4% 4.7% 4.1% Small Soil 75.7% 94.2% - - Medium Soil 67.7% 94.8% - -
Large Soil N/A N/A - - Rumen 70.8% 95.0% 17.5% 14.8% Human Gut 74.4% 99.4% - -
SOAPdenovo Assembler
Cov. of UF by F Cov. of F by UF Cov. of RG by UF Cov. of RF by F Simulated
94.0% 97.0% 3.0% 2.9% Small Soil 86.8% 96.1% - - Medium Soil 82.4% 96.0% - -
Large Soil 78.9% 94.5% - - Rumen 85.2% 97.8% 14.9% 13.6% Human Gut 85.4% 99.3%
- -

Table 4: Total number of abundant (greater than 50x) highly connective sequences incorporated into unfiltered assemblies | Velvet | SOAPdenovo | MetaIDBA
---|---|---|---
Small Soil | 0 (0.0%) | 41 (0.0%) | 8,717 (0.1%)
Medium Soil | 32,328 (0.1%) | 852 (0.0%) | 23,881 (0.1%)
Large Soil | 643,071 (0.3%) | 279,519 (0.1%) | N/A
Rumen | 45,721 (0.2%) | 14,858 (0.1%) | 33,046 (0.1%)
Human Gut | 4,661,447 (3.4%) | 1,749,347 (1.3%) | 5,528,054 (4.0%)
Simulated | 5,118 (1.4%) | 0 (0.0%) | 5,480 (1.5%)
Table 5: Annotations (against 112 reference genomes) of highly-connecting (HC) sequences identified in the simulated metagenome. | Number of HC sequences with annotation
---|---
ABC transporter-like protein | 306
Methyl-accepting chemotaxis sensory transducer | 210
ABC transporter | 173
Elongation factor Tu | 94
Chemotaxis sensory transducer | 51
ABC transporter ATP-binding protein | 44
Diguanylate cyclase/phosphodiesterase | 36
ATPase | 36
S-adenosyl-L-homocysteine hydrolase | 36
Adenosylhomocysteine And downstream NAD binding | 36
Ketol-acid reductoisomerase | 34
S-adenosylmethionine synthetase | 34
Elongation factor G | 34
ABC transporter ATPase | 33
Table 6: Annotations (against NCBI-nr database) of highly-connecting (HC) sequences identified in thee three soil, rumen, and human gut metagenomes. | Number of HC sequences with annotation
---|---
Translation elongation factor/GTP-binding protein LepA | 11
S-adenosylmethionine synthetase | 8
Aspartyl-tRNA synthetase | 8
Malate dehydrogenase | 7
V-type H(+)-translocating pyrophosphatase | 6
Acyl-CoA synthetase | 6
NAD synthetase | 5
Ribonucleotide reductase of class II | 4
Ribityllumazine synthase | 4
Heavy metal translocating P-type ATPase, copA | 3
GyrB | 3
Glutamine amidotransferase chain of NAD synthetase | 3
ChaC family protein | 3
Figure 1: The extent to which average local graph density varies by read
position is shown for the lump of various datasets.

Figure 2: The extent to which highly connecting k-mers (solid lines) and the
subset of highly abundant (greater than 50) k-mers (dashed lines) are present
at specific positions within sequencing reads for various metagenomes. Figure
3: When incorporated into an assembly, abundant (greater than 50 times),
highly connecting sequences (k-mers) were disproportionately present at the
ends of contigs. The total fraction of highly connecting k-mers which are
incorporated into each contig binned region. Figure 4: Rank abundance plot of
5-mers present in abundant, highly connected sequences in various datasets.