CBNplot
Chapter 1 CBNplot: Bayesian network plot for enrichment analysis results
The R package to plot Bayesian network inferred from expression data based on the enrichment analysis results including clusterProfiler or ReactomePA results (Wu et al. 2021; Yu and He 2016). It makes use of libraries including clusterProfiler, ReactomePA, bnlearn, graphite and depmap (Killian and Gatto 2021; Scutari 2010; Sales et al. 2012). The description of functions and several use cases are depicted in this book using GSE133624, which contains RNA-Seq data of bladder cancer and adjacent normal bladder tissues (Chen et al. 2019).
1.1 The preprocessing and DEG identification of GSE133624
In this book, the data from GSE133624 was used for the demonstrative purpose, and the parameters can be changed. First, DESeq2
is used to identify DEGs (Love, Huber, and Anders 2014).
library(DESeq2)
## Load dataset and make metadata
= read.table("GSE133624_reads-count-all-sample.txt", header=1, row.names=1)
counts = sapply(colnames(counts), function (x) substring(x,1,1))
meta = data.frame(meta)
meta colnames(meta) = c("Condition")
<- DESeqDataSetFromMatrix(countData = counts,
dds colData = meta,
design= ~ Condition)
## Prefiltering
<- rowSums(counts(dds) < 10) > dim(meta)[1]*0.9
filt <- dds[!filt,]
dds
## Perform DESeq2()
= DESeq(dds)
dds = results(dds, pAdjustMethod = "bonferroni")
res
## apply variance stabilizing transformation
= vst(dds, blind=FALSE)
v = assay(v) vsted
## Plot PCA of VST values
::plotPCA(v, intgroup="Condition")+
DESeq2theme_bw()
## Define the input genes, and use clusterProfiler::bitr to convert the ID.
= subset(res, padj<0.05)
sig = clusterProfiler::bitr(rownames(sig), fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Hs.eg.db)$ENTREZID
cand.entrez
## Perform enrichment analysis (ORA)
= ReactomePA::enrichPathway(gene = cand.entrez)
pway = clusterProfiler::enrichGO(cand.entrez, ont = "MF", OrgDb = org.Hs.eg.db)
pwayGO
## Convert to SYMBOL
= setReadable(pway, OrgDb=org.Hs.eg.db)
pway = setReadable(pwayGO, OrgDb=org.Hs.eg.db)
pwayGO
## Store the similarity
= enrichplot::pairwise_termsim(pway)
pway
## Define including samples
= rownames(subset(meta, Condition=="T")) incSample
Additionally, for the gene set enrichment analysis, log2 fold changes are obtained.
= clusterProfiler::bitr(rownames(res), fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Hs.eg.db)
allEntrez $ENSEMBL <- rownames(res)
res<- merge(data.frame(res), allEntrez, by="ENSEMBL")
lfc <- lfc[order(lfc$log2FoldChange, decreasing=TRUE),]
lfc <- lfc$log2FoldChange
geneList names(geneList) <- lfc$ENTREZID
<- ReactomePA::gsePathway(geneList) pwayGSE
The mean gene count in the pathway, which is used in the inference.
<- subset(pway@result, p.adjust<0.05)
sigpway paste(mean(sigpway$Count), sd(sigpway$Count))
FALSE [1] "29.6976744186047 22.0886288379294"