RegEnrich: an R package for gene regulator enrichment analysis

In RegEnrich: Gene regulator enrichment analysis

---

knitr::opts_chunk$set(collapse = TRUE, 
                      comment = "#>", 
                      fig.width = 5, 
                      fig.height = 4.5, 
                      fig.align = "center",
                      echo=TRUE, 
                      warning=FALSE, 
                      message=TRUE, 
                      tidy.opts=list(width.cutoff=80), 
                      tidy=FALSE)
rm(list = ls())
gc(reset = TRUE)
options(max.print = 200, width = 110)

Introduction

This package is a pipeline to identify the key gene regulators in a biological process, for example in cell differentiation and in cell development after stimulation. Given gene expression data and sample information, there are four major steps in this pipeline: (1) differential expression analysis; (2) regulator-target network inference; (3) enrichment analysis; and (4) regulators scoring and ranking.

In this tutorial, we are showing you how to perform RegEnrich analysis by starting with a quick example, followed by detailed explanation in each step and three case studies.

---

A quick example

To illustrate how to use RegEnrich pipline, here we simply show the basic steps, assuming you have had all input data and the parameters ready.

Including RegEnrich library

library(RegEnrich)

Initializing RegenrichSet object

object = RegenrichSet(expr = expressionMatrix, # expression data (matrix)
                      colData = sampleInfo, # sample information (data frame)
                      reg = regulators, # regulators
                      method = "limma", # differentila expression analysis method
                      design = designMatrix, # desing model matrix
                      contrast = contrast, # contrast
                      networkConstruction = "COEN", # network inference method
                      enrichTest = "FET") # enrichment analysis method

Runing four major steps and obtaining results

# Differential expression analysis
object = regenrich_diffExpr(object)
# Regulator-target network inference
object = regenrich_network(object)
# Enrichment analysis
object = regenrich_enrich(object)
# Regulator scoring and ranking
object = regenrich_rankScore(object)

# Obtaining results
res = results_score(object)

The code can be even simpler if pipe (%>%) is used.

# Perform 4 steps in RegEnrich analysis
object = regenrich_diffExpr(object) %>% 
  regenrich_network() %>% 
  regenrich_enrich() %>% 
  regenrich_rankScore()

res = results_score(object)

---

RegenrichSet object initialization

The RegenrichSet object

RegEnrich package is programmed in a style of S4 object system. The most fundamental class in RegEnrich is RegenrichSet, and majority functions are working with this class in the following analysis steps. So the RegenrichSet object must be initialized prior to RegEnrich analysis. The RegenrichSet object can be easily initialized by RegenrichSet function, which require several input data along with other parameters.

Input data

There are three fundamental input data for RegEnrich pipline.

Expression data

The first one is expression data (expr), which is a table (matrix) with m rows (genes or proteins) and n columns (samples). Here the expression data can be of genes, measured by microarray or RNA sequencing, and can be of proteins, measured by mass spectrometry, etc.

Here we downloaded the gene expression data (RNA-sequencing data in FPKM format) that is from [@bouquet2016longitudinal] in GEO database (GSE63085). And only 52 samples were included in Lyme_GSE63085 dataset and can be loaded in the following way.

data(Lyme_GSE63085)
FPKM = Lyme_GSE63085$FPKM

Althouth here RNA-sequencing data is reprensented in FPKM (Fragments Per Kilobase of transcript per Million) format, we do recomment raw read count format. To further work on this FPKM data, we convert the FPKM data (plus 1) into logarithm to the base 2.

log2FPKM = log2(FPKM + 1)
print(log2FPKM[1:6, 1:5])

Sample information

The second input data is sample information (colData), which is also a table (data.frame), showing which samples (rows) belonging to which groups or having what features (columns). Here we use the sample information for the 52 samples in Lyme_GSE63085 dataset.

sampleInfo = Lyme_GSE63085$sampleInfo
head(sampleInfo)

Regulators

The third input is the regulators in the studied organisms. If the organism is human, RegEnrich by default provided transcrpiton factors and co-factors in humans as the regulators which are obtained from [@han2015trrust; @marbach2016tissue; and @liu2015regnetwork].

data(TFs)
head(TFs)

You can define your own regulators in RegEnrich. For example, for the studies in other organisms such as mouse, yeast, drosophila and arabidopsis thaliana, you can use the transcription (co-)factors in the following links, and cite corresponding paper:

Mouse, Yeast, Drosophila, and Arabidopsis thaliana

But one thing to keep in mind, the type of names or IDs of regulators should be the same as those of genes in the expression data. For example, if ENSEMBL gene ID is used in the expression data, then the regulators should be represented by ENSEMBL gene ID as well.

Other parameters

In addition to previous 3 most fundamental input data, other parameters for RegEnrich analysis can be initialized here as well. These parameters include 3 groups.

First, parameters to perform differential expression analysis, such as method (differential test method), minMeanExpr (the threshold to remove the lowly expressed gene), design (design model matrix or formula), reduced (reduced model matrix or formula), contrast (to extract the coefficients in the model), etc. Here we consider the effect of different patients and time (week in sample information table) on gene expression. To identify the differential genes related to time, we can simply use LRT_LM method (likelihood retio test on two linear model) to perfrom differential expression analysis. So the corresponding parameters are defined by:

method = "LRT_LM"
# design and reduced formulae
design = ~ patientID + week
reduced = ~ patientID

Second, parameters to perform regulator-target network inference, such as networkConstruction (network inference method), topNetPercent (what percentage of the top edges to retain), etc. Here we use COEN method (weighted gene coexpression network) and other default parameters to inference regulator-target network.

networkConstruction = "COEN"

Third, parameters to perform enrichment analysis, such as enrichTest (enrichment method), minSize (minimum number of targets of regulator), etc. Here we use FET method (Fisher's exact test) and other default parameters to perform enrichment analysis.

enrichTest = "FET"

The detailed explaination of these parameters can be found in the help page of RegenrichSet function, which can be viewed by ?RegenrichSet. Unlike expression data and sample information data, these parameters can be re-specified in the later steps of RegEnrich analysis.

Initializing RegenrichSet object

To reduce the running time, we consider the first 2000 genes, and remove genes with mean log2FPKM <= 0.01.

object = RegenrichSet(expr = log2FPKM[1:2000, ],
                      colData = sampleInfo,
                      method = "LRT_LM", 
                      minMeanExpr = 0.01,
                      design = ~ patientID + week,
                      reduced = ~ patientID,
                      networkConstruction = "COEN", 
                      enrichTest = "FET")
print(object)

---

Differential expression analysis

In this step, the major goal is to obtain differential p-values and log2 fold changes of gene expression between different conditions. There are couple of packages being developed to perform differential expression analysis, such as DESeq2, edgeR, limma, DSS, EBSeq, and baySeq. The full tutorials of these packages have been already provided as vignettes or other documentation in these packages. Here, we provide a wraper function, regenrich_diffExpr, which allows you to choose either "Wald_DESeq2", "LRT_DESeq2", "limma", or "LRT_LM" to perform the differential expression analysis on the RegenrichSet obejct.

Use the parameters initialized in the RegenrichSet object

Since RegenrichSet object is initialized with method = "LRT_LM", regenrich_diffExpr function performs differential expression analysis using likelihood ratio test on two linear models that are specified by design formula and reduced formula.

object = regenrich_diffExpr(object)
print(object)
print(results_DEA(object))

LRT_LM method is implemented for data with complicated experiment designs, in which it is less meaningful to calculate log2 fold changes. In the current version of RegEnrich, calculating the log2 fold change between conditions is not implemented in LRT_LM method. And the log2 fold changes of all genes are set to 0 by default.

Re-specify the parameters for differential expression analysis

Here, parameters to perform differential expression analysis can be re-specified in regenrich_diffExpr function. Below, we show an example of using limma to obtain differential expressed genes, by changing parameters.

object2 = regenrich_diffExpr(object, method = "limma", coef = "week3")
print(object2)
print(results_DEA(object2))

More detailed explaination of regenrich_diffExpr function can be accessed by ?regenrich_diffExpr.

---

Regulator-target network inference

RegEnrich provides two computational methods, COEN and GRN, to infer regulator-target network based on expression data. For COEN method, the weighted co-expression network is constructed using WGCNA package, and then the regulator-target network is the robust subnetwork of weighted co-expression network, and the nodes and edges are removed if they are not connected to any regulators. With respect to GRN, it infers regulator-target network using random forest algorithm. This method was initially described in GENIE3 package, and it is modified in RegEnrich to support parallel computing and to control the model accuracy. Regulator-target network inferences using COEN and GRN methods are shown bellow.

COEN (based on WGCNA)

regenrich_network is the function to perform regulator-target network inference. Since the networkConstruction = "COEN" parameter has been set during the RegenrichSet initialization, so by default RegEnrich constructs a COEN network.

set.seed(123)
object = regenrich_network(object)
print(object)

What happens under the hood in above codes is updating object@topNetP slot, which is a Topnetwork class. RegEnrich provides a function to access this slot.

# TopNetwork class
print(results_topNet(object))

Please note that since the oganism of GSE63085 dataset is Homo sapien, the regulators used in RegEnrich by default are obtained from [@han2015trrust; @marbach2016tissue; and @liu2015regnetwork]. And we are using gene names in the expression data, so the regulators here are also represented by gene names.

# All parameters are stored in object@paramsIn slot
head(slot(object, "paramsIn")$reg)

Since network inference is generally very time-consuming, we suggest saving the RegenrichSet object with the regulator-target network in case of using it next time.

# Saving object to 'fileName.Rdata' file in '/folderName' directory
save(object, file = "/folderName/fileName.Rdata")

A more detailed explaination can be found in the help pages of regenrich_network (see ?regenrich_network) and TopNetwork-class (see ?"TopNetwork-class").

GRN (based on random forest)

Alternatively, you can build a gene regulatory network (GRN) by setting networkConstruction = "GRN" parameter in regenrich_network function. To accelarate computing, you can set number of CPU cores and random seeds using BPPARAM parameter. Here you can control the accuracy of network inferance by minR which are computed based on out-of-bag estimation in random forest. Please note that the lower minR is set, the less edges and potential less regulators are retained.

### not run
library(BiocParallel)
# on non-Windows computers (use 2 workers)
bpparam = register(MulticoreParam(workers = 2, RNGseed = 1234))
# on Windows computers (use 2 workers)
bpparam = register(SnowParam(workers = 2, RNGseed = 1234))

object3 = regenrich_network(object, networkConstruction = "GRN", 
                           BPPARAM = bpparam, minR = 0.3)
print(object3)
print(results_topNet(object3))
save(object3, file = "/folderName/fileName3.Rdata")

User defined network

It is also possible to provide a regulator-target network, which is obtained somewhere else. For example, this network can be constructed based on the relation network of transcription factors and their binding targets. Here we assigned the constructed COEN regulator-target network (a Topnetwork object) in object variable to object2 variable.

network_user = results_topNet(object)
print(class(network_user))
regenrich_network(object2) = network_user
print(object2)

It is also fine to provide a 3-column table (data.frame object) of network edges, in which the first column is regulators, the second column is targets, and the third column is edge weight (reliability).

network_user = slot(results_topNet(object), "elementset")
print(class(network_user))
regenrich_network(object2) = as.data.frame(network_user)
print(object2)

---

Enrichment analysis

RegEnrich provides two methods to perform enrichment analysis, i.e. Fisher's exact test (FET) and gene set enrichment analysis (GSEA). Both methods are implemented in regenrich_enrich function. regenrich_enrich function updates object@resEnrich slot, which is a Enrich class. RegEnrich provides results_enrich function to access this slot.

Fisher's exact test (FET)

Since the enrichTest = "FET" parameter has been set during the RegenrichSet initialization, so by default RegEnrich performs enrichment analysis using FET method.

object = regenrich_enrich(object)
print(results_enrich(object))
# enrich_FET = results_enrich(object)@allResult
enrich_FET = slot(results_enrich(object), "allResult")
head(enrich_FET[, 1:6])

Gene set enrichment analysis (GSEA)

Since the enrichTest = "FET" parameter has been set during the RegenrichSet initialization, but enrichTest = GSEA parameter can be re-specified in regenrich_enrich function to perform enrichment analysis using GSEA method. Typically, GSEA is slower than FET method, especially when the number of reg is large and the regulator-target network is also large. Reducing the number of permutation (nperm, default is 10,000) can be a good trial to have a look at preliminary results.

set.seed(123)
object2 = regenrich_enrich(object, enrichTest = "GSEA", nperm = 5000)
print(results_enrich(object))
# enrich_GSEA = results_enrich(object2)@allResult
enrich_GSEA = slot(results_enrich(object2), "allResult")
head(enrich_GSEA[, 1:6])

You can compare the order of enriched regulators obtained by FET and GSEA methods using plotOrders function.

plotOrders(enrich_FET[[1]][1:20], enrich_GSEA[[1]][1:20])

---

Regulator scoring and ranking

The RegEnrich score is a summarized information from both differential expression analysis and regulator enrichment analysis for regulators. This step of RegEnrich analysis is done by regenrich_rankScore function.

Above all, the differential expression analysis is perormed by LRT_LM method, regulator-target network is infered by COEN method, and enrichment analysis is performed by FET method, so the scores and ranking summarize the importance of regulators by considering regulatory interactions in the studied biological process.

object = regenrich_rankScore(object)
results_score(object)

The expression of regulator and its targets can be viewed using following code.

plotRegTarExpr(object, reg = "ARNTL2")

Note that the previous analysis is a tutorial showing you how to perform basic RegEnrich analysis. As you known this analysis is based on only first 2000 genes, the real key regulators should be different from the previous results.

RegEnrich can work with different types of dataset, such as microarray data, RNAseq raw read count data, and mass spectrometriy proteomic data. The following section shows you 2 case studies of using RegEnrich to work with these 2 types of datasets.

---

Case studies

Case 1: Microarray (single-channel) data

Background

This case study analyzes gene expression changes of primary human hepatocytes after 6 and 24 h treatment with interferon alpha (IFN-α). The gene expression was examined using single-channel Affymetrix Human Genome U133 Plus 2.0 Arrays. And the raw data and normalized data is available in GEO database under GSE31193 accession ID.

Reading the data

There are several ways to read the data. You can download the normalized data file, decompress it, and then read it using read.csv function.

Reading the raw data (.cel files) and then normalize it using other normalization method is also possible. As there is a simpler way to read data from GEO database, this method is not included in this vignette.

# Download all .cel files from GEO database (GSE31193 dataset) to current folder
# and then decompress it to GSE31193_RAW folder.
library(affy)
abatch = ReadAffy(celfile.path = "./GSE31193_RAW/")

# Data normalization
eset <- rma(abatch)

# Annotation
library(annotate)
library(hgu133plus2.db)

ID <- featureNames(eset)
Symbol <- getSYMBOL(ID, "hgu133plus2.db")
fData(eset) <- data.frame(Symbol=Symbol)

The simplist way to read the data is using GEOquery library. To use this library, you must have this package installed.

if (!requireNamespace("GEOquery"))
 BiocManager::install("GEOquery")

library(GEOquery)
eset <- getGEO(GEO = "GSE31193")[[1]]

This dataset contains samples treated with IL28B, but here we are focusing on only control samples and samples after 6 and 24 h IFN-α treatment.

# Samples information
pd0 = pData(eset)
pd = pd0[, c("title", "time:ch1", "agent:ch1")]
colnames(pd) = c("title", "time", "group")
pd$time = factor(c(`n/a` = "0", `6` = "6", `24` = "24")[pd$time],
                 levels = c("0", "6", "24"))
pd$group = factor(pd$group, levels = c("none", "IFN", "IL28B"))
pData(eset) = pd

# Only samples without or after 6 and 24 h IFN-α treatment
eset_IFN = eset[, pd$group %in% c("none", "IFN")]

# Order the samples based on time of treatment
pData(eset_IFN) = pData(eset_IFN)[order(pData(eset_IFN)$time),]

# Rename samples
colnames(eset_IFN) = pData(eset_IFN)$title

# Probes information
probeGene = fData(eset_IFN)[, "Gene Symbol", drop = FALSE]

Here to simplify the analysis, if there are multiple probes matching a gene, we use only one probe with higher average expression value to represent this gene.

probeGene$meanExpr = rowMeans(exprs(eset_IFN))
probeGene = probeGene[order(probeGene$meanExpr, decreasing = TRUE),]

# Keep a single probe for a gene, and remove the probe matching no gene.
probeGene_noDu = probeGene[!duplicated(probeGene$`Gene Symbol`), ][-1,]

data = eset_IFN[rownames(probeGene_noDu), ]
rownames(data) = probeGene_noDu$`Gene Symbol`

Because the speed of network infernece is highly influenced by the number of genes, to quickly illustrate how to use RegEnrich, here we randomly take only 5,000 genes for analysis. If you would like to see the real result in the analysis, then the following data subsetting step should be discarded.

set.seed(1234)
data = data[sample(1:nrow(data), 5000), ]

RegEnrich analysis

Here we would like to know which regulators play key roles in primary human hepatocytes after 24 h treatment with IFN-α.

expressionMatrix = exprs(data) # expression data
# rownames(expressionMatrix) = probeGene_noDu$`Gene Symbol`
sampleInfo = pData(data) # sample information

design = ~time
contrast = c(0, 0, 1) # to extract the coefficient "time24"

data(TFs)
# Initializing a RegenrichSet object
object = RegenrichSet(expr = expressionMatrix, # expression data (matrix)
                      colData = sampleInfo, # sample information (data frame)
                      reg = unique(TFs$TF_name), # regulators
                      method = "limma", # differentila expression analysis method
                      design = design, # desing fomula
                      contrast = contrast, # contrast
                      networkConstruction = "COEN", # network inference method
                      enrichTest = "FET") # enrichment analysis method

# Perform RegEnrich analysis
set.seed(123)

# This procedure takes quite a bit of time.
object = regenrich_diffExpr(object) %>%
  regenrich_network() %>%
  regenrich_enrich() %>%
  regenrich_rankScore()

res = results_score(object)
res

plotRegTarExpr(object, reg = "STAT1")

Case 2: RNAseq read count data

Background

Here we show how to apply RegEnrich on the RNAseq data by analyzing Kang et al's monocyte-macrophage-IFN stimulation dataset ( GSE130567). There are multiple experiment conditions in this study. But here we would like to focus on partial samples in which monocytes were cultured with 10 ng/ml human macrophage colonystimulating factor (M-CSF) in the presence (IFN-γ-primed macrophages) or absence (resting macrophages) of 100 U/ml human IFN-γ for 48 h. RNA were extracted and reverse transcripted followed by sequencing (50 bp, paired-end) using Illumina HiSeq 4000. Sequenced reads were mapped to reference human genome (hg19 assembly) using STAR aligner with default parameters. We will use the raw HT-seq counts for the RegEnrich analysis.

Reading the data

Since the sample information and raw read counts data are archived seperately in GEO database. First, we can read the sample information using GEOquery package.

library(GEOquery)
eset <- getGEO(GEO = "GSE130567")[[1]]
pdata = pData(eset)[, c("title", "geo_accession", "cultured in:ch1", "treatment:ch1")]
colnames(pdata) = c("title", "accession", "cultured", "treatment")
pData(eset) = pdata

# Only samples cultured with M-CSF in the presence or absence of IFN-γ
eset = eset[, pdata$treatment %in% c("NT", "IFNG-3h") & pdata$cultured == "M-CSF"]

# Sample information
sampleInfo = pData(eset)
rownames(sampleInfo) = paste0(rep(c("Resting", "IFNG"), each = 3), 1:3)
sampleInfo$treatment = factor(rep(c("Resting", "IFNG"), each = 3),
                              levels = c("Resting", "IFNG"))

Second, download read count file and decompose into a temporary folder.

tmpFolder = tempdir()
tmpFile = tempfile(pattern = "GSE130567_", tmpdir = tmpFolder, fileext = ".tar")
download.file("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE130567&format=file",
              destfile = tmpFile, mode = "wb")
untar(tmpFile, exdir = tmpFolder)
files = untar(tmpFile, list = TRUE)
filesFull = file.path(tmpFolder, files)

Then read the raw read counts in these files.

dat = list()
for (file in filesFull){
  accID = gsub(".*(GSM\\d{7}).*", "\\1", file)
  if(accID %in% sampleInfo$accession){
    zz = gzfile(file, "rt")
    zzdata = read.csv(zz, header = FALSE, sep = "\t", skip = 4, row.names = 1)
    close(zz)
    zzdata = zzdata[,1, drop = FALSE] # Extract the first numeric column
    colnames(zzdata) = accID
    dat = c(dat, list(zzdata))
  }
}
edata = do.call(cbind, dat)

edata = edata[grep(".*[0-9]+$", rownames(edata)),] # remove PAR locus genes
rownames(edata) = substr(rownames(edata), 1, 15)
colnames(edata) = rownames(sampleInfo)

# Retain genes with average read counts higher than 1
edata = edata[rowMeans(edata) > 1,]

Similar to the case 1, here we randomly take only 5,000 genes to quickly illustrate how to use RegEnrich, but to see the real result from the analysis, you should neglect the following step.

set.seed(1234)
edata = edata[sample(1:nrow(edata), 5000), ]

RegEnrich analysis

expressionMatrix = as.matrix(edata) # expression data

design = ~ treatment
reduced = ~ 1

data(TFs)
# Initializing a RegenrichSet object
object = RegenrichSet(expr = expressionMatrix, # expression data (matrix)
                      colData = sampleInfo, # sample information (data frame)
                      reg = unique(TFs$TF), # regulators
                      method = "LRT_DESeq2", # differentila expression analysis method
                      design = design, # desing fomula
                      reduced = reduced, # reduced
                      networkConstruction = "COEN", # network inference method
                      enrichTest = "FET") # enrichment analysis method

# Perform RegEnrich analysis
set.seed(123)

# This procedure takes quite a bit of time.
object = regenrich_diffExpr(object) %>%
  regenrich_network() %>%
  regenrich_enrich() %>%
  regenrich_rankScore()

res = results_score(object)
res$name = TFs[res$reg, "TF_name"]
res

plotRegTarExpr(object, reg = "ENSG00000028277")

---

Session info

sessionInfo()

---

References

RegEnrich documentation built on March 7, 2021, 2 a.m.