# Wikidata-lite for
Knowledge Extraction and Exploration

Phuc Nguyen  National Institute of Informatics, Japan
<EMAIL_ADDRESS>Hideaki Takeda  National Institute of Informatics, Japan
<EMAIL_ADDRESS>
###### Abstract

Wikidata is the largest collaborative general knowledge graph supported by a
worldwide community. It includes many helpful topics for knowledge exploration
and data science applications. However, due to the enormous size of Wikidata,
it is challenging to retrieve a large amount of data with millions of results,
make complex queries requiring large aggregation operations, or access too
many statement references. This paper introduces our preliminary works on
Wikidata-lite, a toolkit to build a database offline for knowledge extraction
and exploration, e.g., retrieving item information, statements, provenances,
or searching entities by their keywords, attributes. Wikidata-lite has high
performance and memory efficiency, much faster than the official Wikidata
SPARQL endpoint for big queries. The Wikidata-lite repository is available at
https://github.com/phucty/wikidb.

###### Index Terms:

Wikidata-lite, Wikidata, Knowledge Graphs

##  I Introduction

Wikidata is the largest collaborative general knowledge graph supported by a
worldwide community. As of November 2022, Wikidata contains 100 million items,
1.4 billion statements and covers many general topics for knowledge
exploration and data science applications111Wikidata Statistics:
https://www.wikidata.org/wiki/Wikidata:Statistics. Wikidata also contains
references that provide additional information, such as sources and URL links
to such statements (1 billion statement references as in the dump November
2022 of Wikidata).

Unfortunately, due to the enormous size of Wikidata, extracting and exploring
knowledge from Wikidata are challenging even using the standard public
Wikidata Query Service (WDQS)222Wikidata Query: https://query.wikidata.org/ or
parsing official released dumps333Wikidata Dumps:
https://dumps.wikimedia.org/wikidatawiki/entities/ of Wikidata. The amount of
information in Wikidata has been growing rapidly 39 times over seven years,
from 3 GB (gzip compressed dump) in 2015 to 118 GB in 2022, and it keeps
increasing as more and more data. Because of the execution timeout limit
problem of Wikidata online query service [1], users cannot make large queries
with millions of results, execute complex queries requiring large aggregation
operations, or access huge statement references. Moreover, hosting a local
copy of Wikidata requires powerful servers (minimum hardware requirement: 16
CPU cores, 128GB RAM, and 1.5T SSD storage [2].)

This paper introduces our preliminary works on Wikidata-lite, a toolkit to
help build a Wikidata database offline from the Wikidata dump to extract,
explore, and analyze Wikidata knowledge [3]. Wikidata-lite aims to be fast and
easy for Wikidata knowledge extraction and exploration for data scientists and
application developers. With these objectives, we implement Wikidata-lite with
key-value stores (lmdb444lmdb: http://www.lmdb.tech/doc/, and python binding
[4]) as the underneath database. It allows fast reading and saves memory since
the entire database is available in a memory map of a hard disk; there is no
need to do extra copy data to memory. The tool also supports fully concurrent
read access from multiple processes and threads.

Wikidata-lite is an ongoing development project, the current release features
are as follows:

  * •

Retrieving information about 100 million entities and 1.4 billion statements
for knowledge exploration. The information includes item labels, descriptions,
aliases in multilingual, and statement information.

  * •

Accessing more than 1 billion references for fact-checking applications.

  * •

Searching entities by their attributes for subset item extraction. For
example, it takes 0.02 seconds to get 5,868,897 males (wd:Q6581072) using
Wikidata-lite while getting a timeout limit error using the online query
service of Wikidata.

The rest of the paper is structured as follows. We describe the related work
in Section II. We introduce the Wikidata-lite toolkit and the current
representative use cases in Section III. Finally, in Section IV, we present
conclusions and discuss the upcoming features of Wikidata-lite.

Figure 1: Wikidata-lite workflow and features

##  II Related Work

Many tools are available to parse, query, and explore knowledge from Wikidata
dumps. KGTK is a knowledge graph toolkit for creating, transforming, merging,
and exploring knowledge graphs [5]. KGTK model knowledge graphs in tab-
separated (TSV) files with four columns: edge-identifier, head, edge-label,
and tail. Wikidata dumps could be converted to KGTK formats with
transformation operations. KGTK provides many graph manipulation operations
such as intersection, subtraction, and joining and other analytics operations
such as calculating node popularities, degrees, shortest paths, or graph
embeddings. The TSV files could be imported to the Kypher processor to query
Wikidata knowledge [6]. The Kyper underneath technique translates the TSV file
into SQL and executes it on SQLite so that it also consumes an enormous memory
if executing a large query. However, the current KGTK release does not support
accessing statement references. Moreover, Wikidata-lite is implemented using
lmdb on key-value storage with memory maps with lmdb, making it faster and
simpler than using SQLite (one more abstraction of SQL data model on top of
storage).

WDumper is a toolkit for creating a topical subset (set of entities related to
particular topics and the relationships between them) from the JSON dump of
Wikidata [7]. Since WDumper allows users to extract and filter a small
sampling of Wikidata, exploring subsets on small devices such as laptops or
personal computers is feasible. The WDumper output is an RDF dump (N-Triple).
Unlike the subset reusability objective of WDumper, Wikidata-lite focuses on
extracting and exploring whole Wikidata knowledge.

There are many other tools to query, transform and analyze RDF knowledge
graphs, such as Neo4J555Neo4j: https://neo4j.com/, Graphy666Graphy:
https://graphy.link/ and RDFlib777RDFlib:
https://rdflib.readthedocs.io/en/stable/. In this work, we only focus on the
JSON dump of Wikidata knowledge graph.

##  III Wikidata-lite

In this section, we first present Wikidata-lite in Section III-A, then the
use-cases are introduced in Section III-B.

###  III-A Architecture

Wikidata-lite helps user access Wikidata offline with high-speed performance
and memory efficiency usage. Figure 1 depicts an overview of Wikidata-lite.

Given a JSON dump file of Wikidata, Wikidata-lite parses the dump as a JSON
object for each item. Then, the item information is transformed into a key-
value object, where keys are item IDs. The values are JSON objects about item
labels, descriptions, aliases, site links, and all statements and references
related to this item. In the data indexing module, the key-value objects are
directly written to the memory map of a hard disk using lmdb library. The
entire Wikidata is stored in memory maps allowing direct access to item
information without copying this information to memory. We also build an
inverted index of Wikidata items for item search with attribute features and
calculate item popularity using the PageRank algorithms. The processing and
indexing steps take two days, and the final lmdb file size is around 200GB
(the compressed size for sharing is 35GB).

The Wikidata retrieval module helps to fast access item information about the
specific query as explained in Section III-B1. Users also can search items by
labels (keywords) as mentioned in Section III-B2. The item search by
attributes allows users to find a subset of entities based on their attributes
as introduced in Section III-B3. We released the source code, a transformed
database, and example codes in this preliminary work. The API and interfaces
will be released soon in future work.

###  III-B Use-cases

This section explains Wikidata-lite use cases, i.e., item information
retrieval, item search by keywords, and item search by attributes.

####  III-B1 Item Information Retrieval

As all item information is stored in memory maps, it is really fast to
retrieve items. The item information includes:

  * •

Labels, descriptions, aliases: Users could get item labels, descriptions, and
aliases in one specific language or all languages available.

  * •

Site links: Users can access all Wikipedia titles and other projects linked to
the item.

  * •

Statements and references: This function will return all statements and
references of the query item.

  * •

Others: Users also can retrieve item types, inverse relations such as the
instance of relationships, the subclass of relationships, or getting
relationships between two Wikidata items.

The item information retrieval of Wikidata-lite is also easy to integrate into
other downstream applications such as entity linking on text or tabular data
[8], [9].

####  III-B2 Item Search by Keywords

We use the MTabES approach [10] to build item search by keywords. MTabES is a
search tool that allows re-ranking between rankings of the algorithms BM25
(using ElasticSearch888ElasticSearch: https://www.elastic.co/), fuzzy search
with edit distances (using SymSpell: Symmetric Delete algorithm [11]), and
PageRank scores. We index all item labels and aliases in multilingual. Users
could access item search APIs and demos at https://mtab.kgraph.jp/mtabes.

####  III-B3 Item Search by Attributes

We implement the item search by attributes inspired from the haswbstatement
feature of Wikibase Cirrus Search999haswbstatement fetures:
https://www.mediawiki.org/wiki/Help:Extension:WikibaseCirrusSearch. The module
returns items with specific values in the statement with specific properties.
However, Wikidata-lite could directly search items with/or without specific
properties as the inputs. For example, getting all scholarly articles that
have relations with DNA, or X-ray diffraction, or molecular geometry, and
Francis Crick and Nature. The corresponding query input is a combination of
query statements of (AND, ”P31”, ”Q13442814”), (OR, None, ”Q7430”), (OR, None,
”Q12101244”), (OR, None, ”Q911331”), (AND, None, ”Q123280”), (AND, None,
”Q180445”). Wikidata-lite returns 46 items in 0.01754 sections for the query.
A query statements contain three parts: Boolean operations (AND/OR/NOT),
property IDs (If there is no property given, the value of property ID is
None), and item IDs.

Another example of the queries that return millions of results is getting all
male (AND, None, Q6581097). It takes 0.02 seconds to get 5,868,897 items,
while the online query service of Wikidata will get a timeout limitation for
this query. Other examples of item search by attributes can be found at
https://github.com/phucty/wikidb/blob/master/example.ipynb.

##  IV Conclusion and Future Work

This paper introduces our preliminary work on Wikidata-lite, a toolkit for
offline extracting and exploring Wikidata knowledge. Wikidata-lite has high
performance and memory efficiency, allowing users to retrieve item
information, statements, and references from Wikidata dump. Users also can
search items by keywords or attributes.

One extension of Wikidata-lite is the ability to handle attribute item search
with natural language text by integrating keyword item search to find correct
item IDs beforehand. We also plan to use Wikidata-lite for data summarization
tasks and natural language generation from Wikidata statements.

On the consumer side, we plan to make Wikidata-lite available as a service
released with REST APIs and an interface that allows users to use all
Wikidata-lite operations directly.

## Acknowledgements

The research was supported by the Cross-ministerial Strategic Innovation
Promotion Program (SIP) Second Phase, “Big-data and AI-enabled Cyberspace
Technologies” by the New Energy and Industrial Technology Development
Organization (NEDO).

## References

  * [1] D. Vrandecic and M. Krötzsch, “Wikidata: a free collaborative knowledgebase,” _Commun. ACM_ , vol. 57, no. 10, pp. 78–85, 2014. [Online]. Available: https://doi.org/10.1145/2629489
  * [2] A. Sanchez, “[wikidata] minimal hardware requirements for loading wikidata dump in blazegraph,” 2019. [Online]. Available<EMAIL_ADDRESS>  * [3] P. Nguyen, “WikiDB: Build a db (key-value store - lmdb style) from wikidata dump,” 2021. [Online]. Available: https://github.com/phucty/wikidb
  * [4] N. Watson, “Universal python binding for the lmdb ’lightning’ database,” 2013\. [Online]. Available: https://github.com/jnwatson/py-lmdb/
  * [5] F. Ilievski, D. Garijo, H. Chalupsky, N. T. Divvala, Y. Yao, C. M. Rogers, R. Li, J. Liu, A. Singh, D. Schwabe, and P. A. Szekely, “KGTK: A toolkit for large knowledge graph manipulation and analysis,” in _ISWC 2020_ , ser. Lecture Notes in Computer Science, vol. 12507. Springer, 2020, pp. 278–293. [Online]. Available: https://doi.org/10.1007/978-3-030-62466-8_18
  * [6] H. Chalupsky, P. A. Szekely, F. Ilievski, D. Garijo, and K. Shenoy, “Creating and querying personalized versions of wikidata on a laptop,” in _Wikidata Workshop@(ISWC 2021)_ , ser. CEUR Workshop Proceedings, vol. 2982, 2021. [Online]. Available: http://ceur-ws.org/Vol-2982/paper-4.pdf
  * [7] S. A. H. Beghaeiraveri, A. J. G. Gray, and F. J. McNeill, “Experiences of using wdumper to create topical subsets from wikidata,” in _Workshop on Knowledge Graph Construction@(ESWC 2021)_ , ser. CEUR Workshop Proceedings, vol. 2873, 2021. [Online]. Available: http://ceur-ws.org/Vol-2873/paper13.pdf
  * [8] V. Cutrona, J. Chen, V. Efthymiou, O. Hassanzadeh, E. Jiménez-Ruiz, J. Sequeda, K. Srinivas, N. Abdelmageed, M. Hulsebos, D. Oliveira, and C. Pesquita, “Results of semtab 2021,” in _SemTab@(ISWC 2021)_ , ser. CEUR Workshop Proceedings, vol. 3103. CEUR-WS.org, 2021, pp. 1–12. [Online]. Available: http://ceur-ws.org/Vol-3103/paper0.pdf
  * [9] P. Nguyen, I. Yamada, N. Kertkeidkachorn, R. Ichise, and H. Takeda, “Semtab 2021: Tabular data annotation with mtab tool,” in _SemTab@(ISWC 2021)_ , ser. CEUR Workshop Proceedings, vol. 3103. CEUR-WS.org, 2021, pp. 92–101. [Online]. Available: http://ceur-ws.org/Vol-3103/paper8.pdf
  * [10] P. Nguyen, I. Yamada, and H. Takeda, “MTabES: Entity search with keyword search, fuzzy search, and entity popularities,” in _SWO 35 (2021)_. Annual Conference of the Japanese Society for Artificial Intelligence,, 2021, pp. 1N4IS1a02–1N4IS1a02. [Online]. Available: https://www.jstage.jst.go.jp/article/pjsai/JSAI2021/0/JSAI2021_1N4IS1a02/_pdf
  * [11] W. Garbe, “SymSpell: Symmetric delete algorithm,” 6 2012. [Online]. Available: hhttps://github.com/wolfgarbe/SymSpell


# Estimating sample-specific regulatory networks

Marieke Lydia Kuijjer 1,2,†, Matthew Tung 1,2,† GuoCheng Yuan1,2, John
Quackenbush1,2,3, Kimberly Glass4111contact<EMAIL_ADDRESS>
† these authors contributed equally to this work 1Department of Biostatistics
and Computational Biology, Dana-Farber Cancer Institute, Boston, MA, USA
2Department of Biostatistics, Harvard School of Public Health, Boston, MA, USA
3Department of Cancer Biology, Dana-Farber Cancer Institute, Boston, MA, USA
4Channing Division of Network Medicine, Department of Medicine, Brigham and
Women’s Hospital, Harvard Medical School, Boston, MA

###### Abstract

Biological systems are driven by intricate interactions among the complex
array of molecules that comprise the cell. Many methods have been developed to
reconstruct network models that attempt to capture those interactions. These
methods often draw on large numbers of measured expression samples to tease
out subtle signals and infer connections between genes (or gene products). The
result is an aggregate network model representing a single estimate for edge
likelihoods. While informative, aggregate models fail to capture the
heterogeneity that is often represented in a population. Here we propose a
method to reverse engineer sample-specific networks from aggregate network
models. We demonstrate the accuracy and applicability of our approach in
several datasets, including simulated data, microarray expression data from
synchronized yeast cells, and RNA-seq data collected from human subjects. We
show that these sample-specific networks can be used to study the evolution of
network topology across time and to characterize shifts in gene regulation
that may not be apparent in the expression data. We believe the ability to
generate sample-specific networks will revolutionize the field of network
biology and has the potential to usher in an era of precision network
medicine.

##  I Introduction

In recent years we have come to appreciate that in many instances, and in
particular in analyzing complex traits and diseases, a single gene or pathway
cannot fully explain a particular cellular state. Instead, biological
processes are better characterized by complex networks whose structures are
altered as the phenotype changes. Studying the pattern of connections between
biological components, and how these structures change between cell states,
can yield new insights into the mechanisms driving disease. However,
accurately reconstructing these networks in a way that captures both the
properties and complexities of each phenotype remains a significant challenge.

Estimating an informative network often relies upon combining information from
large quantities of data (most commonly gene expression data) representing a
potentially very diverse underlying sample population. This means that even in
datasets where samples may represent a spectrum of phenotypes, most widely-
used network inference methods use gene expression data to estimate an
“aggregate” network De Smet and Marchal (2010); Marbach et al. (2012). These
methods generally begin by calculating a statistic for each gene pair based on
shared information across a set of input samples. These scores are frequently
augmented in some way and are then used to infer the presence or absence of
edges between the genes in the network. In this context, the heterogeneity in
the underlying samples is generally seen as essential to correctly estimating
a (single) network model given the variance in the data. However this
consensus averaging may ignore the multiple different underlying regulatory
networks that might be active in the individual samples. Recognizing this
diversity is important. It is crucial that we examine the collection of
networks to fully appreciate the diversity of processes at work in biological
systems.

Here we propose a method to reverse engineer sample-specific networks from an
aggregate network model. The basic approach depends on comparing a network
deduced from a population to a related network for the population missing a
single sample and then estimating the contribution of that sample. What this
means is that we can use aggregate network models to “extract” specific
networks for each of the individual input samples. We call this approach
LIONESS (Linear Interpolation to Obtain Network Estimates for Single Samples).
As will be shown in the examples below, this approach allows us to
characterize the regulatory processes active in individual samples and to
model how regulatory networks change over time, both of which have many
potential applications in precision medicine and health and biomedical
research. LIONESS is based on relatively straight-forward linear algebra and
is independent of the network inference method used, suggesting the approach
can be generalized to other single sample inference problems.

###  I.1 LIONESS: Linear Interpolation to Obtain Network Estimates for Single
Samples

We begin by suggesting that an “aggregate” network predicted from a set of N
samples can be thought of as the average of individual component networks
reflecting the contributions from each member in the input sample set.
Mathematically, this means that the weight of an edge, ${e}_{ij}^{({\alpha})}$
between two nodes ($i$ and $j$) in an aggregate network derived using all
samples ($\alpha$) can be modelled as the linear combination of the weight of
that edge across a set of networks, where each network in the set corresponds
to a sample used to reconstruct the aggregate model, and
${w}_{s}^{({\alpha})}$ represents the relative contribution of sample ($s$) in
that model:

${e}_{ij}^{(\alpha)}=\displaystyle\sum_{s=1}^{N}w_{s}^{(\alpha)}e_{ij}^{(s)},\text{
where }\displaystyle\sum_{s=1}^{N}w_{s}^{(\alpha)}=1$ (1)

Similarly, the edge weight in a second network reconstructed using all but one
of the samples ($\alpha-q$), can be written as:

${e}_{ij}^{(\alpha-q)}=\displaystyle\sum_{s\neq
q}^{N}w_{s}^{(\alpha-q)}e_{ij}^{(s)},\text{ where }\displaystyle\sum_{s\neq
q}^{N}w_{s}^{(\alpha-q)}=1$ (2)

By comparing these two aggregate network models, we can solve for the network
for an individual sample $q$. Subtracting the above equations we find:

$\begin{split}{e}_{ij}^{(q)}&={e}_{ij}^{(\alpha)}-{e}_{ij}^{(\alpha-q)}\\\
&={w}_{q}^{(\alpha)}{e}_{ij}^{q}+\displaystyle\sum_{s\neq
q}^{N}(w_{s}^{(\alpha)}-w_{s}^{(\alpha-q)}){e}_{ij}^{(s)}\\\
&={w}_{q}^{(\alpha)}{e}_{ij}^{q}-{w}_{q}^{(\alpha)}\displaystyle\sum_{s\neq
q}^{N}{w}_{s}^{(\alpha-q)}{e}_{ij}^{(s)}\end{split}$ (3)

The network specific to sample $q$ in terms of the aggregate networks is then:

${e}_{ij}^{(q)}=\frac{1}{w_{q}^{(\alpha)}}(e_{ij}^{(\alpha)}-e_{ij}^{(\alpha-q)})+e_{ij}^{(\alpha-q)}$
(4)

In summary, we find that the edge scores for a given individual network ($q$)
are equal to the difference in edge scores for an aggregate network
constructed using all the samples, and an aggregate network reconstructed
using all but the sample of interest, multiplied by a scaling factor and added
to the edge scores of the network reconstructed using all but the sample of
interest (Figure 1). In the following analysis we give samples equal weight
(${w}_{q}^{(\alpha)}={1}/{N}$) although one could, in principle, weight
samples differently based on the quality of the data for individual samples or
some other measure. A more detailed version of the LIONESS derivation is
provided in the Supplemental Material (Section VI).

Figure 1: Visual illustration of how LIONESS estimates the network for a
single sample based on two aggregate network models, one reconstructed using
all biological samples in a given dataset, and one using all except the sample
of interest ($q$).

##  II Results

###  II.1 LIONESS accurately and reproducibly predicts networks using in
silico data

To test LIONESS we created a dataset where the underlying networks
corresponding to each input expression sample are known (Supplemental Figure
1). To begin, we created 10000 “gold standard” networks with 40 nodes and
random edges but the same degree distribution. For each network, we generated
1000 random initial states and applied a Boolean model (see Methods) to
determine the corresponding attractors Wuensche (1998). We then averaged over
the all states defined within these 1000 attractors to generate “expression”
values for the 40 nodes (which represent genes) in the network.

We used this in silico data to evaluate the ability of LIONESS to accurately
predict an individual sample network and explore how this prediction might be
influenced by the number of “background” samples included when reconstructing
the aggregate models. To begin, we used Pearson correlation to calculate
aggregate networks based on the in silico expression data. Each of these
models is a complete graph with edge-weight estimates for every possible
network interaction. We then applied LIONESS (Equation 29) to reconstruct all
individual sample networks, and evaluated each of the predicted networks by
comparing with the original “gold standard” networks and calculating the Area
Under the Receiver Operator Characteristic curve (AUCROC, or more simply AUC).
We repeated this multiple times and using a varying number of “background”
samples. We observe that as we include larger numbers of samples in the
aggregate network model, the accuracy of the single-sample networks predicted
by LIONESS increases, and that these are consistently more similar to the
“true” network than the aggregate models from which they were derived (Figure
2A).

Figure 2: Evaluation of LIONESS’ ability to recover known single-samples
networks in an in silico dataset. (A) A plot of the mean and standard
deviation of AUC values calculated when estimating single-sample networks from
aggregate models that include information from varying numbers of background
samples. (B) The median and inter-quartile range of the similarity (calculated
using the Spearman correlation) between pairs of networks, as a function of
the number of background samples used to calculate the aggregate model. For
the plots in (A) and (B) each point represents the results across $1000$
random selections of samples. (C) The distribution of AUC values that result
when we evaluate each of the $10000$ LIONESS-predicted single-sample networks
using their corresponding original “gold-standard” networks. The values are
almost always greater than $0.5$, indicating overall accurate prediction of
the true sample-specific networks. (D) The difference in AUC values when
evaluating the $10000$ LIONESS-predicted networks versus evaluating the
aggregate network reconstructed using all samples. This difference is almost
always greater than zero, indicating that the LIONESS-predicted single-sample
networks are more reflective of the underlying true networks than the
aggregate network from which they are derived.

We also tested whether the LIONESS networks are reproducible. To do this we
compared two single-sample networks that represent the same expression sample,
but which were derived from two aggregate network models that incorporate
fully independent sets of background samples. We observe increasing
reproducibility as we increase the number of background samples used to
reconstruct the aggregate models, with near identity for large background
sample sizes (Figure 2B). As a control, we also compared two single-sample
networks representing different gene expression samples that were derived
using the same set of background samples in the aggregate model. We see almost
no similarity between these networks, especially with increasing numbers of
background samples.

Finally, we investigated the generalizability of LIONESS by estimating single-
sample networks from aggregate models derived using several common network
reconstruction approaches, including Pearson correlation co-expression, PANDA
Glass et al. (2013), mutual information, and CLR Faith et al. (2007) (for more
information, see Methods). We show the distribution in AUC values for the
single-sample network predictions (Figure 2C) and the difference in this AUC
value compared to the AUC of the aggregate model evaluated using the same
“gold standard” (Figure 2D). We also report the median of these values and the
percent of single-sample networks whose AUC is greater than the AUC of the
overall aggregate model (Figure 2C–D). We find that LIONESS is able to
accurately predict single-sample networks for all four of the approaches, with
median AUC values between $0.6$ and $0.7$ (Figure 2C). In all cases the
majority of the single-sample networks reconstructed using LIONESS are more
predictive of the true underlying networks than the aggregate from which they
were derived. Of the four network inference methods, PANDA networks were the
most specific, with over $99\%$ of the single-sample network models more
representative of the true underlying network than the aggregate network
(Figure 2D).

###  II.2 Estimating single-sample networks using experimental data from yeast

We next tested LIONESS using experimental data from cell-cycle synchronized
yeast cells. We downloaded gene expression data (GEO accession, GSE4987;
Pramila et al. (2006)) consisting of dye-swap technical replicates measured
every five minutes for 120 minutes and ma-normalized each replicate Yang et
al. (2007); we then removed probe sets with missing data, batch-corrected
using ComBat Johnson et al. (2007), averaged probe sets mapping to the same
ORF annotation and quantile-normalized the resulting gene-by-sample matrix of
expression values. We note that the 105 minute time point was excluded in both
replicates due to poor hybridization performance Pramila et al. (2006).

We used four different approaches (Pearson Correlation, PANDA Glass et al.
(2013), mutual information, and CLR Faith et al. (2007)) to reconstruct
aggregate networks for this dataset and applied LIONESS to estimate the
networks for each of the individual samples. The correlation between each pair
of the estimated sample-specific networks is shown in the first column of
Figure 3A (R1&R2-from-R1&R2). We see that network estimates for the same
technical replicate are highly similar, as evidenced by the strong diagonal in
the upper-right and lower-left square of each comparison; additional structure
is also evident in off-diagonal similarities that reflect the fact that the
time course data represents more than one cell cycle.

Figure 3: Analysis of LIONESS networks predicted for $48$ expression samples
collected across a yeast cell-cycle time course experiment. (A) LIONESS was
used to predict networks for each sample in the expression dataset by applying
four different aggregate network reconstruction approaches. For each approach
we built the aggregate models either using all samples (R1&R2 from R1&R2), or
only the samples from the same technical replicate (R1-from-R1 & R2-from-R2).
The Spearman correlation was used to evaluate how similar these networks are
to each other. (B) A bar graph showing the change in AUC values for the
aggregate models predicted for each approach, compared to a baseline. For
R1&R2-from-R1&R2 (red bars) we have a single aggregate network. For R1-from-R2
& R2-from-R2 we have two aggregate models (one for each replicate), which we
average before computing the AUC (blue bars). In this analysis, the baseline
AUC values for the Pearson correlation, MI and CLR networks was set equal to
$0.5$. Since PANDA begins with a prior network derived from motif structure,
we use that AUC value ($0.68$) as the baseline in order to have a more fair
comparison with the other approaches.

To test if strong reproducibility was due to including replicates in the
expression data, we also ran LIONESS separately on each individual replicate.
This analysis produced 24 single-sample networks estimated using only the data
in replicate one, and $24$ single-sample networks estimated using only the
data in replicate two (R1-from-R1 & R2-from-R2). The correlation between these
networks is shown in the second column of Figure 3A. As before, we observe
strong reproducibility between technical replicates. Interestingly, even
though we have corrected for batch effects in our expression data, several of
the methods, especially CLR, appear to be sensitive to the “background” data
used.

Although we have shown the networks estimated by applying LIONESS to each of
four reconstruction approaches are reproducible, it is also important that
these networks are also accurate. To test the accuracy of each of the four
network reconstruction methods, we compared the aggregate networks to ChIP-
chip information, and tested whether the aggregate network models could
predict transcription factor binding. For three of the four network
reconstruction approaches we observe AUC values close to $0.5$ (Figure 3B).
The exception is PANDA, which had an AUC above $0.7$.

PANDA differs from the other approaches in that it begins with a “prior” gene
regulatory structure (based on motif data; see Methods) and uses the Pearson
correlation as a proxy to gene co-regulation by the same upstream regulator,
rather than a direct regulatory interaction between a regulator and a
downstream target. We also note that, in addition to having overall higher AUC
values, the improvement in AUC values compared to a null state (the AUC of the
original “seed” motif prior for PANDA and AUC$=0.5$ for the other approaches)
is greatest for PANDA (Figure 3B). Based on these results, in the following
analysis we focus on the single-sample networks derived using PANDA as the
underlying aggregate approach. Results for the other reconstruction approaches
are presented in Supplemental Figure 2.

###  II.3 Single-sample networks show periodic structure across the cell cycle

Next we averaged sample networks representing the same time point in each of
the two replicate network reconstructions and compared the results to the
ChIP-chip data to evaluate the accuracy of each of the resulting network
models. We recognize that this is not a true evaluation of the sample-specific
networks, since we expect the network structure to change throughout the cell
cycle. In Figure 4A we show the difference in the AUC values for the each of
the PANDA+LIONESS-predicted single-sample networks when compared to the AUC of
the motif prior network that was given as an input to PANDA. Although the
increases are modest, each of the single-sample networks has an overall AUC
greater than the prior network. There appears to be an oscillatory pattern in
these AUC values across the time course, with higher values at greater times
post the initial synchronization. This makes sense as cells will gradually
become unsynchronized, better representing the cell population used for the
ChIP-chip data.

Figure 4: Characterizing networks across the yeast cell-cycle. (A) The
difference in the AUC for each sample-specific network predicted by
PANDA+LIONESS versus the AUC of the motif prior used to seed the PANDA
reconstruction approach. The same ChIP-chip data was used as the basis for
each evaluation. (B) A heat map of the edge weights for the $1000$ most
variable edges across the sample-specific network models. Rows are Z-score
normalized for visualization purposes. (C) The average expression of genes
targeted by the four transcription factors that were identified as regulatory
nodes of the $1000$ top most variable edges as well as the average weight of
high-confidence edges that extend between those transcription factors and
their target genes.

We next tested whether these single-sample networks could provide insight into
gene regulation and dynamic cellular network properties. We identified the
$1000$ edges with the highest variability across the individual networks and
plotted those edge weights in Figure 4B. We again observe strong oscillatory
patterns that correspond closely with the oscillations in AUC. It is worth
noting that all these highly variable edges originate from one of four
transcription factors, each of which is important in regulating the yeast cell
cycle (MBP1, SWI4, SWI6, and STB1) Ho et al. (1999).

We examined the genes for which there is strong evidence of targeting by these
transcription factors (average edge weight across all LIONESS networks greater
than zero). In Figure 4C we plot the average weight of these high-evidence
interactions for each regulating transcription factor and the average
expression of their target genes. It is immediately apparent that oscillation
in edge weights occurs at exactly twice the frequency of the oscillation in
gene expression, and that the expression oscillation has a period
approximately equal to the yeast cell cycle.

To understand this result we point out that PANDA interprets correlation in
target gene behavior as indicative of upstream regulation, emphasizing edges
between a regulator and its targets either when those targets are all
coordinately expressed (activated), or coordinately not-expressed (de-
activated or repressed). Thus high edge weight can be interpreted as evidence
for information flow, which may be the result of the physical presence of a
transcription factor (TF) regulating its downstream targets, but could also
reflect a strong lack of that TF. Indeed, by using correlation-based metrics
to infer network relationships many reconstruction approaches operate under a
similar mathematical framework. Although this observation may change how we
interpret the edges in the estimated single-sample networks, the “turn on/turn
off” behavior is exactly what one would predict given how PANDA estimates
networks and is further evidence that LIONESS is extracting realistic single-
sample networks.

###  II.4 Reconstructing single-sample networks for human lymphoblastoid cell-
lines

Lastly, we wanted to test LIONESS in the context of modeling human gene
regulatory networks. Based on the results we observed when applying four
different aggregate reconstruction approaches in conjunction with LIONESS (see
Figures 2–3), we chose to focus on PANDA as our aggregate network
reconstruction approach in this analysis. PANDA incorporates prior network
information, which we constructed by using Haystack Pinello et al. (2014) to
scan the promoter regions of human genes for JASPAR transcription factor
binding motifs Mathelier et al. (2013).

For our expression data we used a set of RNA-seq experiments performed across
multiple immortalized lymphoblastoid cell lines representing different
individuals Degner et al. (2012). We downloaded raw fastq files from the
Pritchard lab website (http://eqtl.uchicago.edu/) and aligned samples to hg19
using Bowtie Langmead et al. (2009); subsequent quality control analysis using
RNA-SeQC DeLuca et al. (2012) excluded two individuals from further analysis
due to low expression profile efficiency scores. This left us with a final set
of $153$ RNA-seq experiments that includes replicates and represents $65$
distinct individuals. We normalized this data using DEseq2 Love et al. (2014).
For additional data processing and normalization information, see Online
Methods.

We used PANDA to reconstruct aggregate gene regulatory network models that
integrate transcription factor binding information and expression data. We
then applied LIONESS to these aggregate models, resulting in $153$ single-
sample networks, one for each of the RNA-seq expression samples. A
hierarchical clustering (complete linkage, Spearman Correlation) of the
network edge weights demonstrates that networks for the same individual nearly
always cluster more strongly with each other than with networks representing
different individuals (Supplemental Figure 3). We also evaluated these
networks using a ChIP-seq “gold standard” and found $91.5\%$ had a higher AUC
compared to the “seed” motif network used by PANDA (see Online Methods). This
analysis demonstrates that even when constructing networks using biological
data from higher-order organisms such as human, the sample-specific networks
predicted by LIONESS are highly reproducible and predictive of known
interactions.

###  II.5 Complex relationships between network targeting and gene expression

In our previous analysis in yeast (Figure 4C) we observed that edges in
single-sample networks reflect strong information flow, which does not
necessarily correspond to physical binding events. To explore this further, we
investigated the relationship between gene targeting and expression in human
networks. First, we averaged single-sample networks that represent the same
individual, resulting in 65 “person-specific” regulatory networks. We then
selected high-evidence regulatory interactions for each transcription factor
(average edge-weight across all single-sample networks greater than zero), and
directly compared the mean edge-weight for these interactions in each of the
single-sample networks to the average expression of the targeted genes in the
original expression samples.

We consistently observe a non-linear relationship between targeting and
expression, with the highest average edge weights occurring when target genes
have either overall high or low expression levels (Figure 5A); this is
consistent with what we observed in our yeast analysis (Figure 4C). Coloring
by the transcription factor expression level in each sample reveals further
patterns; some transcription factors appear to be acting primarily as
activators (increased target gene expression upon increased TF expression and
targeting) or as repressors (decreased target gene expression upon increased
TF expression and targeting). However, the relationship between transcription
factor and target genes is not always simple, indicating that other regulatory
mechanisms, such as co-activators, post-translational modifiers, or epigenetic
mechanisms, are likely playing an important role in mediating these regulatory
events.

Figure 5: Comparison of gene regulation, gene expression and DNase
hypersensitivity data. (A) For six representative transcription factors, the
mean expression of target genes and the mean weight of the edges targeting
those genes across the $65$ samples are plotted. For each sample, the
expression of the TF is shown as a color, scaled to the normal distribution
for visualization purposes. (B) A cartoon illustrating how high edge weights
and thus regulatory activity is not necessarily equivalent to the presence of
a physical interaction. (C) The distribution of the Spearman correlation
values when comparing gene-targeting (calculated by combining LIONESS
predictions with TF expression; $k^{(L+e)}$) and the significance level of
DNase hypersensitivity in a gene’s promoter across all the samples. We also
show the percentage of genes whose targeting positively correlates with DNase
hypersensitivity when targeting is calculated using only the LIONESS-predicted
edge weight ($k^{(L)}$: no expression considered) or a combination of
expression and motif information ($k^{(m+e)}$). These percentages were
calculated either using all $12424$ genes included in our network model, or
for the set of $3488$ genes with a DNase-peak called in all $65$ samples (open
genes).

###  II.6 Increased network targeting corresponds to open chromatin

There is available DNase hypersensitivity profiling of cell lines representing
these same $65$ individuals Degner et al. (2012). As we did not use this data
when modeling the LIONESS networks, we next sought to leverage this
information to investigate how network structures might reflect epigenetic
state. To begin, we downloaded the data from the Pritchard lab website
(http://eqtl.uchicago.edu/) and called DNase “peaks” for each sample. When a
peak fell within the promoter region of a gene, we associated that gene-sample
pair with a score reflecting the significance level of the associated peak
call. We found $12424$ genes with a DNase promoter-peak in at least one sample
and $3488$ with a promoter-peak in all samples. For more information on the
DNase data processing, see Online Methods.

A DNase hypersensitivity peak represents a region of open chromatin that is
often presumed to be occupied by a set of regulatory proteins, including
transcription factors. Consequently, we wanted to determine if differences in
chromatin state between the $65$ individuals is reflected by alterations in
transcription factor targeting across our single-sample networks. Based on the
analysis presented in Figure 4C and 5A, we hypothesized that physical binding
of a transcription factor to the promoter region of a target gene would be
best represented by combining: (1) the edge weight of that interaction in our
single-sample network models, since this value indicates whether information
is flowing between that transcription factor and target gene; and (2) the
expression level of the transcription factor, since this value indicates
whether the TF is physically present in the cell (Figure 5B). Therefore, we
combined this information (see Online methods), resulting in a set of
expression-modified edge-weights for each sample.

We used the sum of the edge-weights associated with each gene to estimate the
number of transcription factors regulating that gene in each of the $65$
single-sample networks ($k^{(L+e)}$). As a control we also calculated gene-
targeting two other ways: (1) using LIONESS edge-weight estimates in the
absence of gene expression information ($k^{(L)}$) and (2) using gene
expression information in the absence of LIONESS-predicted edge-weights
($k^{(m+e)}$); for the second measure we combined transcription factor
expression in each sample with the motif information used for PANDA’s prior.
Finally, to evaluate the association of network targeting with chromatin
state, for each gene we calculated the Spearman correlation between these
targeting-values across the networks and the significance scores of that
gene’s promoter-DNase across the corresponding cell-lines. We find that gene-
targeting in the expression-modified LIONESS model ($k^{(L+e)}$) is very
strongly correlated with promoter-DNase events (Figure 5C). This association
is greater than when using only expression and motif information
($k^{(m+e)}$), demonstrating that the LIONESS approach provides additional
information on chromatin state not apparent in the data used to seed the
algorithm.

###  II.7 Differential-targeting of genes highlights important biological
processes

Finally, we wanted to determine if there are common structures across these
single-sample regulatory networks that might be reflective of important
biological processes. We performed a hierarchical clustering (complete
linkage; Spearman Correlation) on the edge-weights in the $65$ single-sample
networks and observe clear sets of samples defined by this network information
(Figure 6A). To evaluate if this clustering would also be identified using the
expression data alone, we also performed a hierarchical clustering using gene
expression values (Figure 6B). We observe a strikingly different grouping of
individuals compared to the one obtained when clustering the networks.

Figure 6: (A) A hierarchical clustering on the edge-weights for $65$
regulatory networks, one for each subject included in the RNA-seq dataset. (B)
An equivalent hierarchical clustering on gene expression values for these $65$
individuals. This clustering is distinct from the one based on the network
edge weights. Subject labels are colored based on the network clustering. (C)
Reactome pathways enriched based on a GSEA analysis using gene-targeting
instead of gene expression and comparing the right and left clusters of
networks presented in (A). No Reactome pathways were identified when comparing
the expression values of genes in the different groups defined by the
hierarchical clustering presented in (B).

To determine if the expression-based clustering highlights important
biological information we compared the expression levels of genes between
groups of individuals and performed Gene Set Enrichment Analysis Subramanian
et al. (2005). Interestingly, although we find many genes that are
differentially expressed ($2620$ with FDR $<0.01$), running GSEA finds no
enrichment for known biological functions. Separately, we ran GSEA to evaluate
if biological functions are associated with differential-targeting patterns
between the two groups of networks shown in Figure 6A Glass et al. (2014).
This analysis identified enrichment for many cellular processes related to
cell proliferation (in the smaller “orange” cluster; $n=18$) and immune
function (in the larger “blue” cluster; Figure 6C).

Unfortunately, there is little publicly available phenotypic information for
the $65$ individuals in our dataset, and the characteristics that we do know
Choy et al. (2008) are not highly associated with either the expression or the
network clusters. However, given our functional enrichment results, we believe
that the regulatory differences we observe between the network clusters is
likely related to differences in cellular growth rate induced by variable
Epstein-Barr Virus (EBV) levels in the cell-lines. EBV is used to transform
human B-cells into immortalized lymphoblastoids and is known to activate
NF$\kappa$B transcriptional response Cahir et al. (1999). Consistent with this
hypothesis, we find the signature “Activation of NF-kappaB in B-cells” highly
targeted in the small, “cell proliferation” cluster (ES$=0.5$,
FDR$=2*10^{-3}$).

These results indicate that evaluating single-sample networks can lend insight
into alterations in the biological processes active in different individuals
even when a similar analysis of the gene expression data does not. Further,
the relevance of the processes identified through differential targeting in
the model suggests that LIONESS+PANDA are capturing essential biological
processes relevant to each sample.

##  III Discussion

Network models that use aggregate statistics to describe regulation across a
population are important and have contributed to our understanding of a wide
range of phenotypes. However they are limited to modeling the biological
processes common across samples, and thus cannot fully characterize any
potential heterogeneity in regulatory processes within set of underlying
samples. This heterogeneity is important. If we think about the simplest
system—a collection of cells from a single organ, we expect each individual
cell to have its own unique gene expression levels and its own unique active
gene regulatory networks. Understanding the diversity and complexity of these
networks is essential if we are to fully understand the population dynamics of
the system. Similarly, in analyzing disease phenotypes, we recognize that each
individual has his or her own unique manifestation of the disease and that
this should be reflected in individual-specific gene regulatory networks. Our
current aggregate network inference methods do not allow us to fully explore
the complexity of the systems we are hoping to understand.

We propose LIONESS as a theoretical framework to reconstruct networks at the
level of the individual sample by taking advantage of the aggregate models
predicted by widely-used network reconstruction methods. The LIONESS approach
not only addresses the problem of phenotype heterogeneity, but also provides a
means of estimating networks when samples of a particular phenotype or disease
subtype are rare. One can pool a diverse collection of samples and then
extract individual sample networks, grouping these afterward based on network
edge weights.

One implicit assumption in the LIONESS framework is that the samples in a
given collection represent the full span of possible regulatory interactions.
Namely, LIONESS assumes that an aggregate model reconstructed without the
sample of interest represents a common network, and that the difference in
this model with the one constructed using all samples represents a
perturbation to that common network. Importantly, in our analysis we observe
that increasing the number of samples used as a “background” when
reconstructing an aggregate model is beneficial to the single-sample networks.

We recognize that the quality and interpretation of the LIONESS networks is
highly dependent on the underlying aggregate network reconstruction algorithm
used. We emphasize that there are many existing network reconstruction methods
and there remains much debate in the field as to the “correct” or “best” one
to use. In this paper we show only four representative gene expression network
reconstruction approaches: Pearson correlation, Mutual Information, Context
Likelihood of Relatedness, and PANDA. These were chosen not because they are
necessarily the best, but because they illustrate network reconstruction
methods that use either a linear (Pearson) or non-linear (mutual information)
correlation measure, and the extensions of those measures to better capture
true regulatory interactions instead of simple correlative effects. However,
within this limited sample, our analysis suggests that applying LIONESS to
aggregate networks reconstructed using PANDA has the greatest potential for
reconstructing accurate network models able to identify relevant network
biology. Finally, although we tested our approach in the context of using gene
expression to reverse-engineer regulatory networks, the linear framework at
the heart of LIONESS is generalizable and could be applied to infer sample-
specific networks of other types or using other data sources and network
inference methods.

Looking ahead, we believe that accurately estimating single-sample networks
provides a unique opportunity to investigate and identify patient-specific
disease processes and has the potential to open the way for more personalized
network medicine. In network medicine applications regulatory models are often
constructed for a disease subtype. Applying LIONESS to these models could
identify the sub-networks that are active in an individual patient, providing
an opportunity to select specific therapeutic interventions based on the
unique characteristics of the gene regulatory networks active in each
individual. Finally and critically for medical applications, we point out that
by reconstructing patient-specific networks we will be able to directly
correlate network responses and features with clinical information, opening
the door for investigating the complexity of disease features within the
context of network structures.

## References

  * De Smet and Marchal (2010) R. De Smet and K. Marchal, “Advantages and limitations of current network inference methods,” Nature Reviews Microbiology 8, 717 (2010).
  * Marbach et al. (2012) D. Marbach, J. C. Costello, R. Küffner, N. M. Vega, R. J. Prill, D. M. Camacho, K. R. Allison, M. Kellis, J. J. Collins, G. Stolovitzky, et al., “Wisdom of crowds for robust gene network inference,” Nature methods 9, 796 (2012).
  * Wuensche (1998) A. Wuensche, “Discrete dynamical networks and their attractor basins,” (1998).
  * Glass et al. (2013) K. Glass, C. Huttenhower, J. Quackenbush, and G.-C. Yuan, “Passing messages between biological networks to refine predicted interactions,” PloS one 8, e64832 (2013).
  * Faith et al. (2007) J. J. Faith, B. Hayete, J. T. Thaden, I. Mogno, J. Wierzbowski, G. Cottarel, S. Kasif, J. J. Collins, and T. S. Gardner, “Large-scale mapping and validation of escherichia coli transcriptional regulation from a compendium of expression profiles,” PLoS biology 5, e8 (2007).
  * Pramila et al. (2006) T. Pramila, W. Wu, S. Miles, W. S. Noble, and L. L. Breeden, “The forkhead transcription factor hcm1 regulates chromosome segregation genes and fills the s-phase gap in the transcriptional circuitry of the cell cycle,” Genes & development 20, 2266 (2006).
  * Yang et al. (2007) Y. Yang, A. Paquet, and S. Dudoit, “marray: Exploratory analysis for two-color spotted microarray data,” Version 1.16 (2007).
  * Johnson et al. (2007) W. E. Johnson, C. Li, and A. Rabinovic, “Adjusting batch effects in microarray expression data using empirical bayes methods,” Biostatistics 8, 118 (2007).
  * Ho et al. (1999) Y. Ho, M. Costanzo, L. Moore, R. Kobayashi, and B. J. Andrews, “Regulation of transcription at thesaccharomyces cerevisiae start transition by stb1, a swi6-binding protein,” Molecular and cellular biology 19, 5267 (1999).
  * Pinello et al. (2014) L. Pinello, J. Xu, S. H. Orkin, and G.-C. Yuan, “Analysis of chromatin-state plasticity identifies cell-type–specific regulators of h3k27me3 patterns,” Proceedings of the National Academy of Sciences 111, E344 (2014).
  * Mathelier et al. (2013) A. Mathelier, X. Zhao, A. W. Zhang, F. Parcy, R. Worsley-Hunt, D. J. Arenillas, S. Buchman, C.-y. Chen, A. Chou, H. Ienasescu, et al., “Jaspar 2014: an extensively expanded and updated open-access database of transcription factor binding profiles,” Nucleic acids research p. gkt997 (2013).
  * Degner et al. (2012) J. F. Degner, A. A. Pai, R. Pique-Regi, J.-B. Veyrieras, D. J. Gaffney, J. K. Pickrell, S. De Leon, K. Michelini, N. Lewellen, G. E. Crawford, et al., “Dnase [thinsp] i sensitivity qtls are a major determinant of human expression variation,” Nature 482, 390 (2012).
  * Langmead et al. (2009) B. Langmead, C. Trapnell, M. Pop, S. L. Salzberg, et al., “Ultrafast and memory-efficient alignment of short dna sequences to the human genome,” Genome Biol 10, R25 (2009).
  * DeLuca et al. (2012) D. S. DeLuca, J. Z. Levin, A. Sivachenko, T. Fennell, M.-D. Nazaire, C. Williams, M. Reich, W. Winckler, and G. Getz, “Rna-seqc: Rna-seq metrics for quality control and process optimization,” Bioinformatics 28, 1530 (2012).
  * Love et al. (2014) M. I. Love, W. Huber, and S. Anders, “Moderated estimation of fold change and dispersion for rna-seq data with deseq2,” Genome biology 15, 550 (2014).
  * Subramanian et al. (2005) A. Subramanian, P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, S. L. Pomeroy, T. R. Golub, E. S. Lander, et al., “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles,” Proceedings of the National Academy of Sciences of the United States of America 102, 15545 (2005).
  * Glass et al. (2014) K. Glass, J. Quackenbush, E. K. Silverman, B. Celli, S. I. Rennard, G.-C. Yuan, and D. L. DeMeo, “Sexually-dimorphic targeting of functionally-related genes in copd,” BMC systems biology 8, 118 (2014).
  * Choy et al. (2008) E. Choy, R. Yelensky, S. Bonakdar, R. M. Plenge, R. Saxena, P. L. De Jager, S. Y. Shaw, C. S. Wolfish, J. M. Slavik, C. Cotsapas, et al., “Genetic analysis of human traits in vitro: drug response and gene expression in lymphoblastoid cell lines,” PLoS genetics 4, e1000287 (2008).
  * Cahir et al. (1999) M. E. Cahir, K. M. Izumi, and G. Mosialos, “Epstein-barr virus transformation: involvement of latent membrane protein 1-mediated activation of nf-kappab.,” Oncogene 18, 6959 (1999).
  * Albert (2005) R. Albert, “Scale-free networks in cell biology,” Journal of cell science 118, 4947 (2005).
  * Clauset et al. (2009) A. Clauset, C. R. Shalizi, and M. E. Newman, “Power-law distributions in empirical data,” SIAM review 51, 661 (2009).
  * Leek et al. (2014) J. T. Leek, W. E. Johnson, H. S. Parker, A. E. Jaffe, and J. D. Storey, “Package ’sva’,” (2014).
  * Harbison et al. (2004) C. T. Harbison, D. B. Gordon, T. I. Lee, N. J. Rinaldi, K. D. Macisaac, T. W. Danford, N. M. Hannett, J.-B. Tagne, D. B. Reynolds, J. Yoo, et al., “Transcriptional regulatory code of a eukaryotic genome,” Nature 431, 99 (2004).
  * Lab (2004) F. Lab, “Regulatory Map formatted for spreadsheet import,”, http://fraenkel.mit.edu/Harbison/release_v24/txtfiles/ (2004), [Online; accessed Jul 2011].
  * Zhang et al. (2008) Y. Zhang, T. Liu, C. A. Meyer, J. Eeckhoute, D. S. Johnson, B. E. Bernstein, C. Nusbaum, R. M. Myers, M. Brown, W. Li, et al., “Model-based analysis of chip-seq (macs),” Genome Biol 9, R137 (2008).
  * Heinz et al. (2010) S. Heinz, C. Benner, N. Spann, E. Bertolino, Y. C. Lin, P. Laslo, J. X. Cheng, C. Murre, H. Singh, and C. K. Glass, “Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and b cell identities,” Molecular cell 38, 576 (2010).
  * Liao et al. (2013) Y. Liao, G. K. Smyth, and W. Shi, “The subread aligner: fast, accurate and scalable read mapping by seed-and-vote,” Nucleic acids research p. gkt214 (2013).
  * Durinck et al. (2005) S. Durinck, Y. Moreau, A. Kasprzyk, S. Davis, B. De Moor, A. Brazma, and W. Huber, “Biomart and bioconductor: a powerful link between biological databases and microarray data analysis,” Bioinformatics 21, 3439 (2005).
  * Durinck et al. (2009) S. Durinck, P. T. Spellman, E. Birney, and W. Huber, “Mapping identifiers for the integration of genomic datasets with the r/bioconductor package biomart,” Nature protocols 4, 1184 (2009).
  * Meyer et al. (2008) P. E. Meyer, F. Lafitte, and G. Bontempi, “minet: Ar/bioconductor package for inferring large transcriptional networks using mutual information,” BMC bioinformatics 9, 461 (2008).
  * Smyth (2004) G. K. Smyth, “Linear models and empirical bayes methods for assessing differential expression in microarray experiments,” Statistical applications in genetics and molecular biology 3, 1 (2004).
  * Liberzon et al. (2011) A. Liberzon, A. Subramanian, R. Pinchback, H. Thorvaldsdóttir, P. Tamayo, and J. P. Mesirov, “Molecular signatures database (msigdb) 3.0,” Bioinformatics 27, 1739 (2011).

##  IV Online methods

###  IV.1 Network reconstruction approaches

Pearson Correlation: Pearson correlation evaluates the degree of a linear
relationship between two variables. Regulatory networks can be reconstructed
by calculating the Pearson correlation coefficient between expression levels
of each TF and each target gene. These coefficients are measures of whether
TFs and target genes are being co-expressed, which may indicate a regulatory
event.

Passing Attributes between Networks for Data Assimilation (PANDA): PANDA
builds regulatory networks by starting with a prior of possible interactions
between TFs and target genes, for example TF motif binding information. PANDA
integrates this regulatory prior with gene expression information and protein-
protein interaction data, using a message passing approach to determine
information flow between the different data types. The message passing
algorithm that is used in PANDA is based on the fact that if expression levels
of two genes correlate, those genes are more likely to be regulated by similar
sets of TFs than two genes that do not show correlation in expression.
Similarly, TFs that can form complexes are more likely to regulate similar
sets of target genes than TFs that cannot bind together.

Mutual Information (MI): MI compares the joint probability distribution of two
variables to the products of their corresponding marginal distributions.
Similar to Pearson correlation, MI is a measure of association between TFs and
target genes, which may indicate a regulatory event. This method does not
assume linearity or continuity of the data used for building the network.

Context Likelihood of Relatedness (CLR): CLR is based on MI, but applies a
double z-score transformation to the MI scores, normalizing each TF-gene
interaction to the background distribution of MI values for each gene, and to
the background distribution of MI values for each TF.

###  IV.2 Generation of the in silico expression data and regulatory networks

To test LIONESS’s ability to reconstruct sample-specific network models we
generated a set of network models and a corresponding associated set of gene
expression profiles. To begin, we created a single “seed” network model with
forty nodes. To approximate the structure of biological networks Albert
(2005), the out-degree of nodes in this seed model were given a power-law
degree distribution (generated using the approach published in Clauset et al.
Clauset et al. (2009), with $\alpha=5$) with their targets selected randomly.
We ensured that the out-degree and in-degree of all nodes in this seed network
was greater than zero and less than forty. The resulting “seed” network had
$81$ total edges and a density around $5\%$.

Next we randomized this “seed” network model, holding the degree fixed. Then
we generated a set of initial Boolean states for each node in the network, and
determined the subsequent states of the nodes using Stouffer’s Z-score method:

$S_{j}^{(t+1)}=round\Bigg{[}CDF^{-1}\Bigg{(}\frac{\displaystyle\sum_{i}Z_{ij}S_{i}^{(t)}}{\displaystyle\sum_{i}S_{i}^{(t)}}\Bigg{)}\Bigg{]}$
(5)

where $CDF^{-1}$ is the inverse cumulative distribution function for the
normal distribution, $Z_{ij}$ is the z-scored weight of an edge from node $i$
to node $j$, and $S_{i}^{(t)}$ is the state of node $i$ at time $t$. This
Boolean model was run until an attractor solution was found. In total we
generated $1000$ random initial Boolean states for each randomized network,
resulting in $1000$ attractor solutions. The expression level of a node in the
randomized network was then estimated as the average across these steady-state
solutions. This entire process was then repeated for other randomized versions
of the “seed” network model. In total we performed $10000$ network
randomizations and generated $10000$ corresponding expression states.

###  IV.3 Processing the yeast cell cycle expression data

GPR files were downloaded from the Gene Expression Omnibus (GEO; accession
GSE4987). Each of two replicates were separately ma-normalized using the
maNorm() function in the “marray” library in R/Bioconductor Yang et al.
(2007). The data was batch-corrected using the ComBat() function in the “sva”
library Leek et al. (2014) and probe-sets mapping to the same gene were
averaged, resulting in expression values for $5088$ genes across fifty
conditions. Two samples (corresponding to the $105$ minute time point) were
excluded for data-quality reasons, as noted in the original publication, and
genes without motif information (see below) were then removed, giving a final
expression data-set containing $48$ samples and $3551$ genes. This data was
quantile-normalized and used in all subsequent analysis.

###  IV.4 Generating the yeast motif prior data for PANDA and yeast ChIP-chip
gold standard

PANDA requires a prior regulatory network structure in addition to gene
expression information. To construct a motif prior network for yeast we
downloaded predicted binding sites for $204$ yeast transcription factors
Harbison et al. (2004); Lab (2004). This data includes $4360$ genes with
tandem promoters. $3551$ of these genes are also covered on the gene
expression array (see above). $105$ total transcription factors in this
dataset target the promoter of one of these $3551$ genes. The motif map
between these $105$ transcription factors and $3551$ target genes was used as
a prior regulatory network input to the PANDA algorithm. A subset of $65$ of
these transcription factors that also had expression information was used to
reduce the size of the Pearson, MI and CLR predicted networks to edges that
extend between transcription factors and genes.

We used ChIP-chip data from Harbison et al. Harbison et al. (2004) to
calculate AUC values for the predicted single-sample yeast regulatory
networks. Targets of transcription factors in this dataset were defined using
the criterion $p<0.001$.

###  IV.5 Processing the human DNase hypersensitivity data

Raw DNAse hypersensitivity data was downloaded from the Pritchard lab website
(http://eqtl.uchicago.edu/dsQTL_data/RAW_DATA/, accessed June 2014). $204$
different samples corresponding to $70$ different cell lines were available
for download. Data were aligned to the hg19/GRCh37 reference genome using
Bowtie Langmead et al. (2009), with options -v 1 and -m 10, allowing for not
more than one mismatch, and suppressing any alignments for reads having more
than ten reportable alignments. Quality control using Bowtie output revealed
that $16$ samples had a high percentage (greater than $80\%$) of failed reads.
We called DNase hypersensitivity peaks using MACS Zhang et al. (2008). Peaks
with significance score of less than $10^{-5}$ were mapped to the nearest gene
using HOMER Heinz et al. (2010). When the peak fell within the promoter region
of the gene, we assigned a score to that sample-gene pair equal to
$-10*log_{10}(p)$; otherwise the sample-gene pair was given a default score of
zero. We then removed the $16$ poor quality samples, as well as samples for
which we did not have good quality RNA-seq data (see below), leaving us with
$177$ samples corresponding to $65$ cell-lines.

To obtain a promoter-DNase score specific for each cell-line, we averaged
technical replicates. We then filtered this data to include only genes for
which we also had RNA-seq data available (see below) and genes that were
present in our motif prior. The result was a matrix of DNase-promoter scores
that included $12424$ genes that had a DNase promoter-peak in at least one
sample (at least one row-entry greater than zero) and $3488$ had a promoter-
peak in all samples (all row-entries greater than zero).

###  IV.6 Processing the human RNA-seq data

RNA sequencing (RNA-Seq) data were downloaded from the Pritchard lab website
(http://eqtl.uchicago.edu/, accessed April 2014). $173$ different samples
corresponding to $74$ different cell lines were available for download. We
aligned all samples to the hg19/GRCh37 reference genome using Bowtie Langmead
et al. (2009), with options -n 3 and -m 1, allowing for not more than three
mismatches in the seed ($28$ bases on the high-quality end of the read, the
default in Bowtie), and suppressing all non-unique alignments. We used RNA-
SeQC DeLuca et al. (2012) to determine the quality of the reads, using an
expression profiling efficiency cut-off of $0.75$. Cell lines NA19119 and
NA18853 fell below this cut-off. Next we used subread Liao et al. (2013) to
count reads, and the subread algorithm featureCounts to assign and summarize
counts to genes. Finally, we removed samples with poor quality reads, and
samples for which we did not have good quality DNase hypersensitivity data
available (see above), leaving us with $153$ samples corresponding to $65$
different cell lines.

We used the DESeq2 Love et al. (2014) package to analyze read counts for these
$153$ samples and adjusted for different library sizes using the
estimateSizeFactors function. Only genes that had raw counts in at least
$50\%$ of all $153$ samples ($21516/57820$, or $37\%$ of all genes) were
retained for further analysis. Correction for gene length was also performed;
each intensity value was divided by the length of the corresponding gene
(defined as the total length of the genomic region covered by the
features/exons) in the “gene_id” meta-feature in featureCounts. Finally,
Ensembl gene ids were converted to HGNC gene symbols ($16901/21516$ genes)
using R package biomaRt Durinck et al. (2005, 2009). This gene list was
subsetted to only include genes for which we found at least one transcription
factor binding site in the adjacent promoter (see below), and for which we had
at least one DNase hypersensitivity peak-call in the adjacent promoter (see
above). This resulted in a matrix of $12424$ HGNC symbols by $153$ expression
samples.

###  IV.7 Generating the human motif prior data for PANDA and human ChIP-seq
gold standard

We downloaded JASPAR motifs (http://jaspar.genereg.net/; Mathelier et al.
(2013)) and then used Haystack Pinello et al. (2014) to scan the entire hg19
genome for these motifs. Of the $205$ motifs in the JASPAR database, only
$158$ had genomic hits that met our significance threshold ($p<10^{-6}$). We
used HOMER (http://homer.salk.edu/homer/ngs/index.html; Heinz et al. (2010))
to get the distances of these motif hits to the nearest transcriptional start
site (TSS). We used these reported distances to parse the motif hits based on
their TSS proximity, keeping only those hits within the “promoter” region of a
gene, which we define as [$-750,+250$] around the TSS. We then filtered for
genes for which we had DNase hypersensitivity data and RNA-seq data available.
This information was used to make a prior transcription factor to gene map
that we used when running PANDA.

To generate a gold-standard to evaluate the predicted network models, we
downloaded ChIP-seq bed files assessing the localization of $39$ transcription
factors which were included in our motif network and have been assayed in
GM12878 by ENCODE (http://genome.ucsc.edu/ENCODE/). We used this location
information to create a “gold standard” transcription factor to gene network,
using HOMER (http://homer.salk.edu/homer/ngs/index.html; Heinz et al. (2010))
to map TF peaks to the nearest gene in a process identical to our mapping of
TF-motif location information, described above.

###  IV.8 Running LIONESS on networks reconstructed using mutual information
and CLR

In order to calculate the mutual information and reconstruct networks based on
the Context Likelihood of Relatedness, we used the build.mim() and clr()
functions (default parameters) within the “minet” package in R/Bioconductor
Meyer et al. (2008). It is worth noting that these algorithms create symmetric
gene by gene matrices of predicted edge scores. Therefore, in the context of
reconstructing single-sample yeast regulatory networks, we reduced these
aggregate networks by selecting only the portion of the predicted gene-by-gene
matrix that corresponds to edges from a transcription factor to a target gene.

###  IV.9 Running LIONESS on networks reconstructed using PANDA

PANDA requires a prior regulatory network structure in addition to gene
expression information. For the in silico data we used an identity matrix
(corresponding to each transcription factor targeting only itself) as the
prior regulatory network. For the yeast and human data a prior regulatory
network was constructed based on transcription sequence-motif information (see
above). It is worth noting that although PANDA can optionally take protein-
protein interaction (PPI) information as an input, to generate a more fair
comparison with other network reconstruction approaches PPI data was never
used in any of the networks generated as a part of this study.

###  IV.10 Calculating gene degree and comparing with DNase hypersensitivity
data

When comparing with the DNase information, we calculated gene degree three
different ways:

$\begin{split}{k}_{j}^{(L)}&=\displaystyle\sum_{i}p_{ij}^{(L)};\\\
{k}_{j}^{(L+e)}&=\displaystyle\sum_{i}{p}_{ij}^{(L)}{p}_{i}^{(e)};\\\
{k}_{i}^{(m+e)}&=\displaystyle\sum_{i}{p}_{ij}^{(m)}{p}_{i}^{(e)}\end{split}$
(6)

where $p_{i}^{(e)}$ is the “probability” that TF $i$ is expressed, calculated
by taking the inverse CDF of the z-score of the TF’s expression compared to
the background of its expression in all other samples; $p_{ij}^{(L)}$ is the
“probability” of the edge from the LIONESS networks, found by taking the
inverse CDF of the predicted edge-weight score (which is in z-score units);
$p_{ij}^{(m)}$ is the “probability” of the edge from the motif data, either
$0$ or $1$ based on whether the motif of TF $i$ was found in the promoter of
gene $j$.

###  IV.11 Clustering networks/expression and running GSEA

Clusters of networks and expression samples were defined based on a
hierarchical clustering in which we row-normalized edge-weights (or gene
expression) across samples, calculated distance based on the Spearman
correlation, and performed a complete-linkage clustering. When performing Gene
Set Enrichment Analysis (GSEA) on the network targeting we calculated a score
for each gene by summing up all predicted edge-weight scores for each gene. We
then used LIMMA analysis Smyth (2004) to determine the differences between the
clusters based on these scores and ranked genes based on the log-fold changes.
We ran a pre-ranked GSEA with $1000$ iterations and observed $68$ significant
Reactome pathway signatures Liberzon et al. (2011) with an FDR $<0.1$ and a
gene set size less than forty. An equivalent clustering/LIMMA/GSEA analysis on
the gene expression data resulted in no significant Reactome pathways.

##  V Supplemental Figures

Supplemental Figure 1: A schematic overview of how we generated in silico
expression data for a set of known underlying gene regulatory network models.
Supplemental Figure 2: Results when applying LIONESS to yeast cell cycle data
using various other aggregate reconstruction approaches. Although LIONESS is a
generalizable approach, these results highlight the importance of selecting a
robust underlying aggregate reconstruction algorithm when applying LIONESS.
Supplemental Figure 3: A hierarchical clustering of the $153$ single-sample
networks predicted for LIONESS in lymphoblastoid cell-lines. Close inspection
reveals that single-sample networks derived from separate expression samples
(technical replicates), but measuring the same cell-line, almost always
cluster together. The AUC of these networks, when using ChIP-seq data for the
GM12878 cell-line as a gold-standard, is shown directly below the clustering.

##  VI Supplemental Material

##  Appendix A Derivation to Find Regulatory Networks for Individual Samples
in a Collection

To begin, we assume that the value of a given edge (${e}_{ij}^{({\alpha})}$)
from a transcription factor ($i$) to a gene ($j$) predicted by a network
reconstruction algorithm using a collection of samples ($\alpha$) is the
linear combination of the value of that edge across networks specific to each
of the input samples ($e_{ij}^{(s)}$), where ${w}_{s}^{({\alpha})}$ represents
the relative contribution of sample ($s$):

${e}_{ij}^{({\alpha})}=\displaystyle\sum_{s=1}^{N}w_{s}^{(\alpha)}e_{ij}^{(s)},\text{
where }\displaystyle\sum_{s=1}^{N}w_{s}^{(\alpha)}=1$ (7)

Given this assumption, we can calculate two “aggregate” networks, one using
all samples (${e}_{ij}^{(\alpha)}$), as described above, and the other using
all but one of the samples (${e}_{ij}^{(\alpha-q)}$):

${e}_{ij}^{(\alpha-q)}=\displaystyle\sum_{s\neq
q}^{N}w^{(\alpha-q)}_{s}e_{ij}^{(s)},\text{ where }\displaystyle\sum_{s\neq
q}^{N}w^{(\alpha-q)}_{s}=1$ (8)

Now, subtracting these two “aggregate” network estimates we get:

$\displaystyle{e}_{ij}^{(\alpha)}-{e}_{ij}^{(\alpha-q)}$
$\displaystyle=\displaystyle\sum_{s=1}^{N}w_{s}^{(\alpha)}e_{ij}^{(s)}-\displaystyle\sum_{s\neq
q}^{N}w^{(\alpha-q)}_{s}e_{ij}^{(s)}$ (9)
$\displaystyle=w^{(\alpha)}_{q}e_{ij}^{(q)}+\displaystyle\sum_{s\neq
q}^{N}w_{s}^{(\alpha)}e_{ij}^{(s)}-\displaystyle\sum_{s\neq
q}^{N}w^{(\alpha-q)}_{s}e_{ij}^{(s)}$ (10)
$\displaystyle=w^{(\alpha)}_{q}e_{ij}^{(q)}+\displaystyle\sum_{s\neq
q}^{N}(w^{(\alpha)}_{s}-w^{(\alpha-q)}_{s})e_{ij}^{(s)}$ (11)

We can then solve for the network specific to a the single sample,
$e_{ij}^{(q)}$:

$\displaystyle e_{ij}^{(q)}$
$\displaystyle=\frac{1}{w^{(\alpha)}_{q}}\Bigg{[}{e}_{ij}^{(\alpha)}-{e}_{ij}^{(\alpha-q)}+\displaystyle\sum_{s\neq
q}^{N}(w^{(\alpha-q)}_{s}-w^{(\alpha)}_{s})e_{ij}^{(s)}\Bigg{]}$ (12)
$\displaystyle=\frac{1}{w^{(\alpha)}_{q}}\Bigg{[}{e}_{ij}^{(\alpha)}-\displaystyle\sum_{s\neq
q}^{N}w^{(\alpha)}_{s}e_{ij}^{(s)}\Bigg{]}$ (13)

If we then stipulate that the weights used to estimate ${e}_{ij}^{(\alpha)}$
are related to the weights used to estimate ${e}_{ij}^{(\alpha-q)}$ by a
constant, $w^{(\alpha)}_{s}=Cw^{(\alpha-q)}_{s}$, which, given Equation 7,
implies, $C=1-w^{(\alpha)}_{q}$, we can calculate the values of the edges in
the single-sample network, $e_{ij}^{(q)}$, in terms of the two “aggregate”
networks:

$\displaystyle e_{ij}^{(q)}$
$\displaystyle=\frac{1}{w^{(\alpha)}_{q}}\Bigg{[}{e}_{ij}^{(\alpha)}-C\displaystyle\sum_{s\neq
q}^{N}w^{(\alpha-q)}_{s}e_{ij}^{(s)}\Bigg{]}$ (14)
$\displaystyle=\frac{1}{w^{(\alpha)}_{q}}\Bigg{[}{e}_{ij}^{(\alpha)}-(1-w^{(\alpha)}_{q}){e}_{ij}^{(\alpha-q)}\Bigg{]}$
(15)
$\displaystyle=\frac{1}{w^{(\alpha)}_{q}}\Bigg{[}{e}_{ij}^{(\alpha)}-{e}_{ij}^{(\alpha-q)}\Bigg{]+{e}_{ij}^{(\alpha-q)}}$
(16)

Finally, we can simplify this equation by (optionally) assuming that each
sample is given equal weight ($w^{(\alpha)}_{q}=\frac{1}{N}$):

$e_{ij}^{(q)}=N\Bigg{[}{e}_{ij}^{(\alpha)}-{e}_{ij}^{(\alpha-q)}\Bigg{]+{e}_{ij}^{(\alpha-q)}}$
(17)

##  Appendix B Application to Pearson Correlation

To begin, we remember that the Pearson correlation ($r$) between two
variables, $X$ and $Y$ can be defined as:

$r=\frac{1}{N-1}\displaystyle\sum_{i}^{N}\Bigg{(}\frac{X_{i}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y}}{S_{Y}}\Bigg{)},\text{
where }\bar{X}=\frac{1}{N}\displaystyle\sum_{i}^{N}X_{i}\text{ and
}S_{X}=\sqrt{\frac{1}{N-1}\displaystyle\sum_{i}^{N}(X_{i}-\bar{X})^{2}}$ (18)

We can then use the Pearson to calculate two “aggregate” correlation networks,
one using all samples (resulting in $r$), and the other using all but one of
the samples (resulting in $r^{\prime}$):

$r^{\prime}=\frac{1}{N-2}\displaystyle\sum_{i\neq
q}^{N}\Bigg{(}\frac{X_{i}-\bar{X^{\prime}}}{S^{\prime}_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y^{\prime}}}{S^{\prime}_{Y}}\Bigg{)},\text{
where }\bar{X^{\prime}}=\frac{1}{N-1}\displaystyle\sum_{i\neq
q}^{N}X_{i}\text{ and
}S^{\prime}_{X}=\sqrt{\frac{1}{N-2}\displaystyle\sum_{i\neq
q}^{N}(X_{i}-\bar{X^{\prime}})^{2}}$ (19)

Now, using the “LIONESS” equation for deriving single samples we see that:

$\displaystyle e_{ij}^{(q)}$
$\displaystyle=N(e_{ij}^{(r)}-e_{ij}^{(r^{\prime})})+e_{ij}^{(r^{\prime})}$
(20) $\displaystyle r_{xy}^{(q)}$
$\displaystyle=N(r_{xy}-r^{\prime}_{xy})+r^{\prime}_{xy}$ (21)
$\displaystyle=N\Bigg{[}\frac{1}{N-1}\displaystyle\sum_{i}^{N}\Bigg{(}\frac{X_{i}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y}}{S_{Y}}\Bigg{)}-\frac{1}{N-2}\displaystyle\sum_{i\neq
q}^{N}\Bigg{(}\frac{X_{i}-\bar{X^{\prime}}}{S^{\prime}_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y^{\prime}}}{S^{\prime}_{Y}}\Bigg{)}\Bigg{]}+\frac{1}{N-2}\displaystyle\sum_{i\neq
q}^{N}\Bigg{(}\frac{X_{i}-\bar{X^{\prime}}}{S^{\prime}_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y^{\prime}}}{S^{\prime}_{Y}}\Bigg{)}$
(22)
$\displaystyle=\frac{N}{N-1}\displaystyle\sum_{i}^{N}\Bigg{(}\frac{X_{i}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y}}{S_{Y}}\Bigg{)}-\frac{N-1}{N-2}\displaystyle\sum_{i\neq
q}^{N}\Bigg{(}\frac{X_{i}-\bar{X^{\prime}}}{S^{\prime}_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y^{\prime}}}{S^{\prime}_{Y}}\Bigg{)}$
(23)
$\displaystyle=\frac{N}{N-1}\Bigg{(}\frac{X_{q}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{q}-\bar{Y}}{S_{Y}}\Bigg{)}+\frac{N}{N-1}\displaystyle\sum_{i\neq
q}^{N}\Bigg{(}\frac{X_{i}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y}}{S_{Y}}\Bigg{)}-\frac{N-1}{N-2}\displaystyle\sum_{i\neq
q}^{N}\Bigg{(}\frac{X_{i}-\bar{X^{\prime}}}{S^{\prime}_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y^{\prime}}}{S^{\prime}_{Y}}\Bigg{)}$
(24)
$\displaystyle=\frac{N}{N-1}\Bigg{(}\frac{X_{q}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{q}-\bar{Y}}{S_{Y}}\Bigg{)}+\displaystyle\sum_{i\neq
q}^{N}\Bigg{[}\Bigg{(}\frac{N}{N-1}\Bigg{)}\Bigg{(}\frac{X_{i}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y}}{S_{Y}}\Bigg{)}-\Bigg{(}\frac{N-1}{N-2}\Bigg{)}\Bigg{(}\frac{X_{i}-\bar{X^{\prime}}}{S^{\prime}_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y^{\prime}}}{S^{\prime}_{Y}}\Bigg{)}\Bigg{]}$
(25)

We also note that (from equations 18 and 19 above):

$\displaystyle\bar{X}$
$\displaystyle=\frac{1}{N}\displaystyle\sum_{i}^{N}X_{i}$ (26)
$\displaystyle=\frac{1}{N}X_{q}+\frac{1}{N}\displaystyle\sum_{i\neq
q}^{N}X_{i}$ (27)
$\displaystyle=\frac{1}{N}X_{q}+\frac{N-1}{N}\bar{X^{\prime}}$ (28)
$\displaystyle=\bar{X^{\prime}}+\frac{1}{N}(X_{q}-\bar{X^{\prime}})$ (29)

And consequently,

$X_{i}-\bar{X}=X_{i}-\bar{X^{\prime}}+\frac{1}{N}(\bar{X^{\prime}}-X_{q})$
(30)

Thus when the difference between $X_{q}$ and $\bar{X^{\prime}}$ is much less
than the total number of samples being considered ($N$),
$\bar{X^{\prime}}\rightarrow\bar{X}$ and $S^{\prime}_{X}\rightarrow S_{X}$.
This is most likely for large values of $N$ in which case the above equation
can be simplified to:

$r_{xy}^{(q)}\rightarrow\frac{N}{N-1}\Bigg{(}\frac{X_{q}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{q}-\bar{Y}}{S_{Y}}\Bigg{)}$
(31)

We can now see that, for Pearson-correlation derived “single-sample” networks,
the linearity assumption holds, given the above limits on $N$ and the
relationship between $X_{q}$ and $\bar{X^{\prime}}$. Specifically, if we use
equation 31 and average over all “single-sample” Pearson-derived networks, we
see:

$\frac{1}{N}\displaystyle\sum_{i}^{N}r_{xy}^{(i)}=\frac{1}{N-1}\displaystyle\sum_{i}^{N}\Bigg{(}\frac{X_{i}-\bar{X}}{S_{X}}\Bigg{)}\Bigg{(}\frac{Y_{i}-\bar{Y}}{S_{Y}}\Bigg{)}$
(32)

which is equal to $r$ (equation 18).