In Feng-Zhang/RNAseqFlow:

knitr::opts_chunk$set(
  collapse = TRUE,
  echo = TRUE,
  comment = "#>",
  warning=FALSE,
  message = FALSE
)
#rmarkdown::render("vignettes/base.Rmd",output_format = "html_document")
set.seed(123)

Differential analysis wrapping DEseq2 functions

library(RNAflow)
input <- create_count_phe()

head(input[[1]])
input[[2]]
dds <- create_DEseq(count_data=input[[1]],col_data=input[[2]],design_names = "condition+type",group_name ="condition",ref_level = "untreated")
res <- DESeq_res(dds,fold_change_line=0.58,adj_pvalue_line =0.05)
DESeq_PCA(dds) # plot PCA
volcano(res) # plot volcano
DEsig <- row.names(res)[res$regulate!="Normal"]

Enrichment analysis wrapping clusterProfiler functions

library(org.Dm.eg.db)
library(clusterProfiler)

The over-representation test based on GO terms

idType(OrgDb = org.Dm.eg.db)
go_res <- enrichGO(DEsig,keyType="ENSEMBL",OrgDb=org.Dm.eg.db)
head(go_res)

The over-representation test based on pathway from KEGG database

DEsig_entrezid <- bitr(DEsig, fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Dm.eg.db)[,2]
kegg_res <- enrichKEGG(DEsig_entrezid,keyType="ncbi-geneid",organism="dme");
head(kegg_res)

Note the keyType parameter should be "kegg" or "ncbi-geneid" human, and "ncbi-geneid" for other organism.

GSEA based on GO terms

gene_list <- res$log2FoldChange
names(gene_list)=row.names(res)
gene_list <- sort(gene_list, decreasing = TRUE)
gsea_go_res <- gseGO(gene_list,keyType="ENSEMBL",OrgDb=org.Dm.eg.db,pvalueCutoff = 0.1)
head(gsea_go_res)

if(nrow(gsea_go_res)>=1){
  gsea_go_res <- setReadable(gsea_go_res, OrgDb = org.Dm.eg.db)
  print(gseaplot(gsea_go_res,geneSetID=gsea_go_res$ID[1],title=paste("BP : ",gsea_go_res$Description[1],sep="")))
}

GSEA based on pathway from KEGG database

gsea_kegg_res <- GSEA_kegg(gene_list,type="ENSEMBL",organism = "dme",db=org.Dm.eg.db,pvalueCutoff=1)
head(gsea_kegg_res)
if(nrow(gsea_kegg_res)>1){
  print(gseaplot(gsea_kegg_res,geneSetID=gsea_kegg_res$ID[1],title=paste("KEGG : ",gsea_kegg_res$Description[1],sep="")))
}