1 Installation

You can install the current release version of NetActivity from Bioconductor with:

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("NetActivity")

You can install the development version of NetActivity from GitHub with:

# install.packages("devtools")
devtools::install_github("yocra3/NetActivity")
devtools::install_github("yocra3/NetActivityData")

2 Introduction

The coordination of multiple genes is required to perform specific biological functions, so grouping gene expression in gene sets provides better biological insights than studying individual genes. We have developed a framework, NetActivityTrain, to encode individual gene expression measurements into gene set scores. This framework is implemented in Nextflow and is based on sparsely-connected autoencoders. In our framework, each gene set is represented by a neuron in the hidden layer, which is connected to only the gene set genes from the input layer. Nonetheless, to ensure a better representation of the gene sets, all gene sets are connected to all genes in the output layer.

Our framework, NetActivityTrain, is implemented in Nextflow, which is useful to train the models but not to apply the models to new data. To overcome this issue, we have implemented the package NetActivity. NetActivity enables to easily apply the gene set representations obtained with NetActivityTrain to new datasets. We have included two models containing GO Biological processes and KEGG pathways trained with GTEx or TCGA data in NetActivityData.

3 Pre-processing

Computation of gene set scores is performed with computeGeneSetScores. The computeGeneSetScores functions requires a SummarizedExperiment with the gene expression data standardized at gene level. The function prepareSummarizedExperiment facilitates the process.

NetActivity can be applied to RNAseq or microarray data. We will exemplify how to process each dataset independently.

3.1 RNAseq

For this vignette, we will show how to process RNAseq using the airway dataset:

library(NetActivity)
library(DESeq2)
library(airway)
data(airway)

The airway dataset contains 63677 genes and 8 samples. Gene expression data is represented as counts. NetActivity functions require normalized data to work. We suggest using Variance Stabilizing Transformation from DESeq2 to normalize gene expression data:

ddsSE <- DESeqDataSet(airway, design = ~ cell + dex)
vst <- varianceStabilizingTransformation(ddsSE)

Once we have normalized the data, we can proceed with the pre-processing. prepareSummarizedExperiment requires two arguments: a SummarizedExperiment and a model trained from NetActivityTrain. The current version of NetActivityData includes a model with GO-BP (Biological Processes) terms and KEGG pathways trained with GTEx or with TCGA. In this vignette, we will use the model trained in GTEx:

out <- prepareSummarizedExperiment(vst, "gtex_gokegg")
## Warning in prepareSummarizedExperiment(vst, "gtex_gokegg"): 50 genes present in
## the model not found in input data. The expression of all samples will be set to
## 0.
out
## class: SummarizedExperiment 
## dim: 63727 8 
## metadata(0):
## assays(1): ''
## rownames(63727): ENSG00000000003 ENSG00000000005 ... ENSG00000278637
##   ENSG00000280789
## rowData names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(10): SampleName cell ... BioSample sizeFactor

prepareSummarizedExperiment performs two steps. First, it checks whether all genes in the model are present in the input SummarizedExperiment. Missing genes are added as a row of 0s. Second, the function standardizes gene expression data by gene. We can compare the genes values before and after prepareSummarizedExperiment to see the effect of normalization:

library(tidyverse)
rbind(assay(vst[1:5, ]) %>%
        data.frame() %>%
        mutate(Gene = rownames(.)) %>%
        gather(Sample, Expression, 1:8) %>%
        mutate(Step = "Before"),
    assay(out[1:5, ]) %>%
        data.frame() %>%
        mutate(Gene = rownames(.)) %>%
        gather(Sample, Expression, 1:8) %>%
        mutate(Step = "After")) %>%
    mutate(Step = factor(Step, levels = c("Before", "After"))) %>%
    ggplot(aes(x = Gene, y = Expression, col = Gene)) +
        geom_boxplot() +
        theme_bw() +
        facet_grid(~ Step, scales = "free") +
        theme(axis.ticks = element_blank(), axis.text.x = element_blank())        
Standardization. Effect of standardization on gene expression values. Before represent the gene expression values passed to prepareSummarizedExperiment. After are the values obtained after prepareSummarizedExperiment.

Figure 1: Standardization
Effect of standardization on gene expression values. Before represent the gene expression values passed to prepareSummarizedExperiment. After are the values obtained after prepareSummarizedExperiment.

In Figure 1, we can see that the original gene expression values are normalized values with different medians and dispersions. After standardization, all gene values are centered to a mean of 0 and have more similar distributions.

Next, we can check the values of a gene not present originally in the data:

assay(out["ENSG00000278637", ]) 
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000278637          0          0          0          0          0
##                 SRR1039517 SRR1039520 SRR1039521
## ENSG00000278637          0          0          0

Notice that functions included in NetActivity require that the SummarizedExperiment and the model shared the same annotation. In this example, both objects were annotated using ENSEMBL annotation, so no additional steps were required. In case the genes in the SummarizedExperiment use a different annotation, see how to convert them in the next section. Notice that the current models included in NetActivityData contain the genes using ENSEMBL annotation.

3.2 Microarray data

For this vignette, we will show how to process microarray data using the Fletcher2013a dataset:

library(limma)
library(Fletcher2013a)
data(Exp1)
Exp1
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 47231 features, 46 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: MA1 MA2 ... MI6 (46 total)
##   varLabels: Sample Time ... Target (10 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: ILMN_1343291 ILMN_1343295 ... ILMN_3311190 (47231
##     total)
##   fvarLabels: PROBEID ENTREZ SYMBOL
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation: illuminaHumanv4

Exp1 is an ExpressionSet containing the information for 47231 genes and 46 samples. Genes are named with Affymetrix ids, but the SYMBOL information is available.

The first step is converting the ExpressionSet to a SummarizedExperiment:

SE_fletcher <- SummarizedExperiment(exprs(Exp1), colData = pData(Exp1), rowData = fData(Exp1))

Second, we will set the rownames to ENSEMBL. For this, we will first map probe ids to SYMBOL and then to ENSEMBL id. Notice that mapping probes to SYMBOL and then to ENSEMBL might result in probes without annotation or duplicated probes with the same ENSEMBL id. For the sake of simplicity, in this vignette we will remove probes without SYMBOL or with duplicated SYMBOLs based on order. Nonetheless, the users are encouraged to apply other criteria in order to get the final set.

library(AnnotationDbi)
library(org.Hs.eg.db)
rownames(SE_fletcher) <- rowData(SE_fletcher)$SYMBOL
SE_fletcher <- SE_fletcher[!is.na(rownames(SE_fletcher)), ]
SE_fletcher <- SE_fletcher[!duplicated(rownames(SE_fletcher)), ]

rownames(SE_fletcher) <- mapIds(org.Hs.eg.db,
    keys = rownames(SE_fletcher),
    column = 'ENSEMBL',
    keytype = 'SYMBOL')
## 'select()' returned 1:many mapping between keys and columns
SE_fletcher <- SE_fletcher[!is.na(rownames(SE_fletcher)), ]
SE_fletcher <- SE_fletcher[!duplicated(rownames(SE_fletcher)), ]
SE_fletcher
## class: SummarizedExperiment 
## dim: 18059 46 
## metadata(0):
## assays(1): ''
## rownames(18059): ENSG00000156508 ENSG00000111640 ... ENSG00000221493
##   ENSG00000221545
## rowData names(3): PROBEID ENTREZ SYMBOL
## colnames(46): MA1 MA2 ... MI4 MI6
## colData names(10): Sample Time ... TreatmentGroups Target

Now, we can pass the resulting SummarizedExperiment to prepareSummarizedExperiment:

out_array <- prepareSummarizedExperiment(SE_fletcher, "gtex_gokegg")
## Warning in prepareSummarizedExperiment(SE_fletcher, "gtex_gokegg"): 793 genes
## present in the model not found in input data. The expression of all samples
## will be set to 0.
out_array
## class: SummarizedExperiment 
## dim: 18852 46 
## metadata(0):
## assays(1): ''
## rownames(18852): ENSG00000156508 ENSG00000111640 ... ENSG00000278637
##   ENSG00000280789
## rowData names(0):
## colnames(46): MA1 MA2 ... MI4 MI6
## colData names(10): Sample Time ... TreatmentGroups Target

4 Computing the scores

Once we have pre-processed the data, we are ready to compute the gene set scores. We will continue with the output of the microarray data (out_array).

computeGeneSetScores requires two arguments: a SummarizedExperiment and a model trained from NetActivityTrain:

scores <- computeGeneSetScores(out_array, "gtex_gokegg")
scores
## class: SummarizedExperiment 
## dim: 1518 46 
## metadata(0):
## assays(1): ''
## rownames(1518): GO:0006734 GO:0008340 ... path:hsa04966 path:hsa04977
## rowData names(4): Term GeneSet Weights Weights_SYMBOL
## colnames(46): MA1 MA2 ... MI4 MI6
## colData names(10): Sample Time ... TreatmentGroups Target

computeGeneSetScores returns a SummarizedExperiment with the gene set scores for all the gene sets present in the model. The colData of the input SummarizedExperiment is preserved, so this object can be directly used for downstream analyses. The rowData contains annotation of the gene sets, which can be useful for downstream analyses:

rowData(scores)
## DataFrame with 1518 rows and 4 columns
##                                 Term       GeneSet
##                          <character>   <character>
## GO:0006734    NADH metabolic process    GO:0006734
## GO:0008340    determination of adu..    GO:0008340
## GO:0014854    response to inactivity    GO:0014854
## GO:0019081         viral translation    GO:0019081
## GO:0019835                 cytolysis    GO:0019835
## ...                              ...           ...
## path:hsa04614 Renin-angiotensin sy.. path:hsa04614
## path:hsa04744      Phototransduction path:hsa04744
## path:hsa04964 Proximal tubule bica.. path:hsa04964
## path:hsa04966 Collecting duct acid.. path:hsa04966
## path:hsa04977 Vitamin digestion an.. path:hsa04977
##                                                  Weights
##                                                   <list>
## GO:0006734           0.0209889, 0.0738816,-0.0206952,...
## GO:0008340           0.1628420,-0.0977037,-0.1346692,...
## GO:0014854          -0.2046938,-0.1480707,-0.0132261,...
## GO:0019081     0.002303295, 0.027539201,-0.000309411,...
## GO:0019835          -0.0330166, 0.0329533, 0.0290691,...
## ...                                                  ...
## path:hsa04614       -0.0728291,-0.0747844,-0.0602768,...
## path:hsa04744       0.14130385,0.13613251,0.00561176,...
## path:hsa04964        0.0112757, 0.0593170,-0.0692757,...
## path:hsa04966       -0.0475978, 0.1188805, 0.0976878,...
## path:hsa04977       -0.0478214, 0.0285422,-0.0713180,...
##                                           Weights_SYMBOL
##                                                   <list>
## GO:0006734           0.0209889, 0.0738816,-0.0206952,...
## GO:0008340           0.1628420,-0.0977037,-0.1346692,...
## GO:0014854          -0.2046938,-0.1480707,-0.0132261,...
## GO:0019081     0.002303295, 0.027539201,-0.000309411,...
## GO:0019835          -0.0330166, 0.0329533, 0.0290691,...
## ...                                                  ...
## path:hsa04614       -0.0728291,-0.0747844,-0.0602768,...
## path:hsa04744       0.14130385,0.13613251,0.00561176,...
## path:hsa04964        0.0112757, 0.0593170,-0.0692757,...
## path:hsa04966       -0.0475978, 0.1188805, 0.0976878,...
## path:hsa04977       -0.0478214, 0.0285422,-0.0713180,...

5 Differential gene set scores analysis

Once the gene set scores are computed, we can test for differential expression using limma. We will run a differential analysis due to treatment, adjusted for time points:

mod <- model.matrix(~ Treatment + Time, colData(scores))
fit <- lmFit(assay(scores), mod) %>% eBayes()
topTab <- topTable(fit, coef = 2:4, n = Inf)
head(topTab)
##            TreatmentE2FGF10 TreatmentE2FGF10PD TreatmentUT       AveExpr
## GO:1990440        1.0255687         0.03220728  -0.9531159 -2.665365e-15
## GO:0045943        0.5426360         0.03697507  -0.8617941 -4.529891e-16
## GO:0000470       -0.3737068        -0.12955486   0.8604381 -1.766401e-16
## GO:2000232       -0.3266975        -0.12581617   0.7801178 -2.485934e-16
## GO:0000083        0.2643651        -0.10021547  -1.0062766 -3.378940e-16
## GO:1902570        0.3195350         0.15336997  -0.5789914 -1.704554e-17
##                   F      P.Value    adj.P.Val
## GO:1990440 64.16935 1.524191e-17 2.313722e-14
## GO:0045943 52.20593 9.428650e-16 7.156345e-13
## GO:0000470 42.27305 5.142105e-14 2.601905e-11
## GO:2000232 31.03742 1.145857e-11 4.348526e-09
## GO:0000083 30.54285 1.494377e-11 4.536928e-09
## GO:1902570 27.47854 8.245091e-11 2.086008e-08

For easing the interpretation, we can add the gene set full name to the results table:

topTab$GeneSetName <- rowData(scores)[rownames(topTab), "Term"]
head(topTab)
##            TreatmentE2FGF10 TreatmentE2FGF10PD TreatmentUT       AveExpr
## GO:1990440        1.0255687         0.03220728  -0.9531159 -2.665365e-15
## GO:0045943        0.5426360         0.03697507  -0.8617941 -4.529891e-16
## GO:0000470       -0.3737068        -0.12955486   0.8604381 -1.766401e-16
## GO:2000232       -0.3266975        -0.12581617   0.7801178 -2.485934e-16
## GO:0000083        0.2643651        -0.10021547  -1.0062766 -3.378940e-16
## GO:1902570        0.3195350         0.15336997  -0.5789914 -1.704554e-17
##                   F      P.Value    adj.P.Val
## GO:1990440 64.16935 1.524191e-17 2.313722e-14
## GO:0045943 52.20593 9.428650e-16 7.156345e-13
## GO:0000470 42.27305 5.142105e-14 2.601905e-11
## GO:2000232 31.03742 1.145857e-11 4.348526e-09
## GO:0000083 30.54285 1.494377e-11 4.536928e-09
## GO:1902570 27.47854 8.245091e-11 2.086008e-08
##                                                                                                                 GeneSetName
## GO:1990440 positive regulation of transcription from RNA polymerase II promoter in response to endoplasmic reticulum stress
## GO:0045943                                                         positive regulation of transcription by RNA polymerase I
## GO:0000470                                                                                           maturation of LSU-rRNA
## GO:2000232                                                                                    regulation of rRNA processing
## GO:0000083                                    regulation of transcription involved in G1/S transition of mitotic cell cycle
## GO:1902570                                                                                protein localization to nucleolus

Many gene sets present high differences due to treatment. We will further explore the top gene set (GO:1990440 or positive regulation of transcription from RNA polymerase II promoter in response to endoplasmic reticulum stress).

Once we have our candidate gene set, we recommend taking a look to the distribution of gene set scores based on our variable of interest:

data.frame(Expression = as.vector(assay(scores["GO:1990440", ])),
    Treatment = scores$Treatment) %>%
    ggplot(aes(x = Treatment, y = Expression, col = Treatment)) +
        geom_boxplot() +
        theme_bw() +
    ylab("NetActivity scores")
GO:1990440 activity scores. GO:1990440 presented the most significant difference due to treatment.

Figure 2: GO:1990440 activity scores
GO:1990440 presented the most significant difference due to treatment.

We can clearly see differences in GO:1990440 activation scores between the different treatment groups in Figure 2. UT samples have the lowest gene set activity scores while E2FGF10 have the highest activation scores.

Nonetheless, notice that the sign of the gene set score is arbitrary and cannot be direc