Introduction to SeqNet
In SeqNet: Generate RNA-Seq Data from Gene-Gene Association Networks

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 5
)

library(SeqNet)
set.seed(12345)

Overview

SeqNet generates random gene-gene association networks and simulates RNA-seq count data from them, as described in Grimes and Datta (2021). A network is built out of overlapping modules that represent pathways, giving it the topological properties (hub genes, community structure) that are characteristic of real gene regulatory networks. Once a network exists, SeqNet can:

assign connection strengths so it defines a valid Gaussian graphical model,
simulate RNA-seq counts whose gene-gene correlations follow that network, either non-parametrically (matching a reference RNA-seq dataset) or parametrically (zero-inflated negative binomial), and
perturb a network to create a second, differential network that can be used for benchmarking differential co-expression or differential network methods.

This vignette walks through that workflow.

Building a random network

random_network() creates a network of p nodes made up of a number of overlapping modules:

nw <- random_network(p = 100, n_modules = 5)
nw

The printed summary reports the number of nodes, edges, and modules, along with global network characteristics (e.g. average degree, clustering coefficient). The network can be visualized directly:

g <- plot_network(nw)

plot_network() returns the layout it used (g), which can be reused so that the same network is always drawn with nodes in the same position. Modules can be highlighted on top of an existing layout:

plot_modules(nw, g)

Assigning connection strengths

A freshly created network only specifies which genes are connected, not how strongly. gen_partial_correlations() assigns edge weights so that the network's association matrix is a valid (positive-definite) partial correlation matrix (i.e. a Gaussian graphical model):

nw <- gen_partial_correlations(nw)
is_weighted(nw)

heatmap_network() visualizes the resulting association matrix:

heatmap_network(nw)

Simulating RNA-seq data

Two functions simulate expression data whose correlation structure follows the network. gen_rnaseq() uses a Gaussian copula: it first draws multivariate-normal data based on the network's partial correlations, then transforms each gene's marginal distribution to match a reference RNA-seq dataset (via the inverse CDF). If no reference is supplied, SeqNet uses a bundled reference dataset that is a subset of the TCGA breast invasive carcinoma cohort:

x <- gen_rnaseq(n = 20, network = nw, verbose = FALSE)$x
dim(x)

Alternatively, gen_zinb() simulates directly from a zero-inflated negative binomial distribution fit to each gene, rather than resampling from the empirical reference distribution:

x_zinb <- gen_zinb(n = 20, network = nw, verbose = FALSE)$x
dim(x_zinb)

Both approaches preserve the correlation structure implied by nw; they differ in how each gene's marginal (univariate) distribution is generated.

Differential networks

perturb_network() creates a modified copy of a network by rewiring connections around one or more hub genes (and, optionally, additional random genes). This simulates the kind of localized rewiring seen between, for example, healthy and diseased tissue:

nw_diff <- perturb_network(nw, n_hubs = 1, n_nodes = 5)
plot_network_diff(nw, nw_diff, g)

The differential network plot colors edges that are unique to each network, making it easy to see where the two networks disagree.

Comparing expression of a gene pair

When comparing simulated (or real) expression data across multiple groups, plot_gene_pair() plots the relationship between two genes, optionally faceted or colored by group:

x1 <- gen_rnaseq(n = 20, network = nw, verbose = FALSE)$x
x2 <- gen_rnaseq(n = 20, network = nw_diff, verbose = FALSE)$x
genes <- colnames(x1)
plot_gene_pair(list(network_1 = x1, network_2 = x2), genes[1], genes[2])

Learning more

Each function's help page (e.g. ?random_network, ?gen_rnaseq, ?perturb_network) documents additional arguments for controlling module size and overlap, network size, and simulation parameters. See citation("SeqNet") for how to cite the package, and Grimes and Datta (2021) for the full methodology.