## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE) ## ----eval = FALSE------------------------------------------------------------- # if (!require("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # # BiocManager::install("NetActivity") ## ----eval = FALSE------------------------------------------------------------- # # install.packages("devtools") # devtools::install_github("yocra3/NetActivity") # devtools::install_github("yocra3/NetActivityData") ## ----warning = FALSE, message = FALSE----------------------------------------- library(NetActivity) library(DESeq2) library(airway) data(airway) ## ----------------------------------------------------------------------------- ddsSE <- DESeqDataSet(airway, design = ~ cell + dex) vst <- varianceStabilizingTransformation(ddsSE) ## ----------------------------------------------------------------------------- out <- prepareSummarizedExperiment(vst, "gtex_gokegg") out ## ----warning = FALSE, message = FALSE----------------------------------------- library(tidyverse) ## ----plot, fig.cap = "Standardization. Effect of standardization on gene expression values. Before represent the gene expression values passed to prepareSummarizedExperiment. After are the values obtained after prepareSummarizedExperiment."---- 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()) ## ----------------------------------------------------------------------------- assay(out["ENSG00000278637", ]) ## ----message = FALSE, warning = FALSE----------------------------------------- library(limma) library(Fletcher2013a) data(Exp1) Exp1 ## ----------------------------------------------------------------------------- SE_fletcher <- SummarizedExperiment(exprs(Exp1), colData = pData(Exp1), rowData = fData(Exp1)) ## ----warning=FALSE, message = FALSE------------------------------------------- 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') SE_fletcher <- SE_fletcher[!is.na(rownames(SE_fletcher)), ] SE_fletcher <- SE_fletcher[!duplicated(rownames(SE_fletcher)), ] SE_fletcher ## ----------------------------------------------------------------------------- out_array <- prepareSummarizedExperiment(SE_fletcher, "gtex_gokegg") out_array ## ----------------------------------------------------------------------------- scores <- computeGeneSetScores(out_array, "gtex_gokegg") scores ## ----------------------------------------------------------------------------- rowData(scores) ## ----------------------------------------------------------------------------- mod <- model.matrix(~ Treatment + Time, colData(scores)) fit <- lmFit(assay(scores), mod) %>% eBayes() topTab <- topTable(fit, coef = 2:4, n = Inf) head(topTab) ## ----------------------------------------------------------------------------- topTab$GeneSetName <- rowData(scores)[rownames(topTab), "Term"] head(topTab) ## ----plotScores, fig.cap = "GO:1990440 activity scores. GO:1990440 presented the most significant difference due to treatment."---- 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") ## ----plotWeight, fig.cap = "Weights of GO:1990440 gene set. The figure represents the weights used for computing the GO:1990440 gene set score. Weights are in absolute value to enable an easier comparison of their magnitude. Positive weights are shown in blue and negative in red.", fig.wide = TRUE---- weights <- rowData(scores)["GO:1990440", ]$Weights_SYMBOL[[1]] data.frame(weight = weights, gene = names(weights)) %>% mutate(Direction = ifelse(weight > 0, "Positive", "Negative")) %>% ggplot(aes(x = gene, y = abs(weight), fill = Direction)) + geom_bar(stat = "identity") + theme_bw() + ylab("Weight") + xlab("Gene") ## ----------------------------------------------------------------------------- sessionInfo()