gsean.R
In gsean: Gene Set Enrichment Analysis with Networks

## ---- fig.align='center', message=FALSE, warning=FALSE, eval=FALSE------------
#  library(gsean)
#  library(TCGAbiolinks)
#  library(SummarizedExperiment)
#  # TCGA LUAD
#  query <- GDCquery(project = "TCGA-LUAD",
#                    data.category = "Gene expression",
#                    data.type = "Gene expression quantification",
#                    platform = "Illumina HiSeq",
#                    file.type  = "normalized_results",
#                    experimental.strategy = "RNA-Seq",
#                    legacy = TRUE)
#  GDCdownload(query, method = "api")
#  invisible(capture.output(data <- GDCprepare(query)))
#  exprs.LUAD <- assay(data)
#  # remove duplicated gene names
#  exprs.LUAD <- exprs.LUAD[-which(duplicated(rownames(exprs.LUAD))),]
#  # list of genes
#  recur.mut.gene <- c("KRAS", "TP53", "STK11", "RBM10", "SPATA31C1", "KRTAP4-11",
#                      "DCAF8L2", "AGAP6", "KEAP1", "SETD2", "ZNF679", "FSCB",
#                      "BRAF", "ZNF770", "U2AF1", "SMARCA4", "HRNR", "EGFR")
#  
#  # KEGG_hsa
#  load(system.file("data", "KEGG_hsa.rda", package = "gsean"))
#  
#  # GSEA
#  set.seed(1)
#  result.GSEA <- gsean(KEGG_hsa, recur.mut.gene, exprs.LUAD, threshold = 0.7)
#  invisible(capture.output(p <- GSEA.barplot(result.GSEA, category = 'pathway',
#                                             score = 'NES', pvalue = 'padj',
#                                             sort = 'padj', top = 20)))
#  p <- GSEA.barplot(result.GSEA, category = 'pathway', score = 'NES',
#                    pvalue = 'padj', sort = 'padj', top = 20)
#  p + theme(plot.margin = margin(10, 10, 10, 75))

## ---- fig.align='center', message=FALSE, warning=FALSE, eval=TRUE-------------
library(gsean)
library(pasilla)
library(DESeq2)
# pasilla count data
pasCts <- system.file("extdata", "pasilla_gene_counts.tsv",
                      package = "pasilla", mustWork = TRUE)
cts <- as.matrix(read.csv(pasCts, sep="\t", row.names = "gene_id"))
condition <- factor(c(rep("untreated", 4), rep("treated", 3)))
dds <- DESeqDataSetFromMatrix(
  countData = cts,
  colData   = data.frame(condition),
  design    = ~ 0 + condition)
# filtering
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
# differentially expressed genes
dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds, contrast = list("conditiontreated", "conditionuntreated"), listValues = c(1, -1))
statistic <- res$stat
names(statistic) <- rownames(res)
exprs.pasilla <- counts(dds, normalized = TRUE)

# convert gene id
library(org.Dm.eg.db)
gene.id <- AnnotationDbi::select(org.Dm.eg.db, names(statistic), "ENTREZID", "FLYBASE")
names(statistic) <- gene.id[,2]
rownames(exprs.pasilla) <- gene.id[,2]

# GO_dme
load(system.file("data", "GO_dme.rda", package = "gsean"))

# GSEA
set.seed(1)
result.GSEA <- gsean(GO_dme, statistic, exprs.pasilla)
invisible(capture.output(p <- GSEA.barplot(result.GSEA, category = 'pathway',
                                           score = 'NES', top = 50, pvalue = 'padj',
                                           sort = 'padj', numChar = 110) + 
  theme(plot.margin = margin(10, 10, 10, 50))))
plotly::ggplotly(p)

## ----sessionInfo--------------------------------------------------------------
sessionInfo()