CBNplot：可视化富集通路中基因的调控关系

说在前面

相信大家在平时做富集分析时都会有这样的一个需求：如果能知道感兴趣的某条通路中各基因的调控关系，那么就能准确识别出hub基因；或者说找到我们感兴趣的基因在这条通路中的上下游调控关系，从而就可以进行后续的实验验证。很多情况下只有想象中是完美的，但是只要感想就会有实现的机会，对于上面说的这个想法就在今年被实现了。

想必国内的生信小伙伴都或多或少的听闻过Y叔的大名，Y叔开发的一系列生信分析软件可谓撑起了国内生信圈的半边天。而今天Immugent介绍的这个软件也是最近由Y叔和京都大学的Yasushi Okuno一同开发的CBNplot，相应的文章发表在Bioinformatics杂志上，篇名为 CBNplot: Bayesian network plots for enrichment analysis。

image.png

关于CBNplot的介绍，生信宝库会以三篇推文并且以代码实操的形式分别介绍其主要功能，下面开始介绍第一部分的用法。

---

代码展示

首先我们先从GEO上下载一个示例数据，算出差异基因后再做富集分析。

library(DESeq2)

## Load dataset and make metadata
counts = read.table("GSE133624_reads-count-all-sample.txt", header=1, row.names=1)
meta = sapply(colnames(counts), function (x) substring(x,1,1))
meta = data.frame(meta)
colnames(meta) = c("Condition")

dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = meta,
                              design= ~ Condition)
## Prefiltering
filt <- rowSums(counts(dds) < 10) > dim(meta)[1]*0.9
dds <- dds[!filt,]

## Perform DESeq2()
dds = DESeq(dds)
res = results(dds, pAdjustMethod = "bonferroni")

## apply variance stabilizing transformation
v = vst(dds, blind=FALSE)
vsted = assay(v)
## Plot PCA of VST values
DESeq2::plotPCA(v, intgroup="Condition")+
    theme_bw()

image.png

## Define the input genes, and use clusterProfiler::bitr to convert the ID.
sig = subset(res, padj<0.05)
cand.entrez = clusterProfiler::bitr(rownames(sig), fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Hs.eg.db)$ENTREZID

## Perform enrichment analysis (ORA)
pway = ReactomePA::enrichPathway(gene = cand.entrez)
pwayGO = clusterProfiler::enrichGO(cand.entrez, ont = "BP", OrgDb = org.Hs.eg.db)

## Convert to SYMBOL
pway = setReadable(pway, OrgDb=org.Hs.eg.db)
pwayGO = setReadable(pwayGO, OrgDb=org.Hs.eg.db)

## Store the similarity
pway = enrichplot::pairwise_termsim(pway)

## Define including samples
incSample = rownames(subset(meta, Condition=="T"))

allEntrez = clusterProfiler::bitr(rownames(res), fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Hs.eg.db)
res$ENSEMBL <- rownames(res)
lfc <- merge(data.frame(res), allEntrez, by="ENSEMBL")
lfc <- lfc[order(lfc$log2FoldChange, decreasing=TRUE),]
geneList <- lfc$log2FoldChange
names(geneList) <- lfc$ENTREZID

pwayGSE <- ReactomePA::gsePathway(geneList)
sigpway <- subset(pway@result, p.adjust<0.05)
paste(mean(sigpway$Count), sd(sigpway$Count))

基于富集分析的结果我们就可以使用CBNplot对我们感兴趣的通路进行展示了。

barplot(pway, showCategory = 15)
#使用bngeneplot函数绘图
bngeneplot(results = pway, exp = vsted, pathNum = 17)
#Change the label for the better readability.
bngeneplot(results = pway, exp = vsted, pathNum = 17, labelSize=7, shadowText=TRUE)

# Show the confidence of direction
bngeneplot(results = pway,
           exp = vsted,
           expSample = incSample,
           pathNum = 13, R = 50, showDir = T,
           convertSymbol = T,
           expRow = "ENSEMBL",
           strThresh = 0.7)

image.png image.png image.png

可以通过参数compareRef=TRUE并指定pathDb，可以将基因之间的关系与参考网络进行比较。默认情况下，两个有向网络的交集以重叠边的数量表示。

library(parallel)
cl = makeCluster(4)
bngeneplot(results = pway,
           exp = vsted,
           expSample = incSample,
           pathNum = 13, R = 30, compareRef = T,
           convertSymbol = T, pathDb = "reactome",
           expRow = "ENSEMBL", cl = cl)
           
bngeneplot(results = pway,
           exp = vsted,
           expSample = incSample,
           pathNum = 15, R = 10, compareRef = T,
           convertSymbol = T, pathDb = "reactome", compareRefType = "difference",
           expRow = "ENSEMBL")

image.png image.png

还可以添加一个barplot来描述边缘的强度和方向(概率)，指定strength plot =TRUE和nStrength。

bngeneplot(results = pway,
           exp = vsted,
           expSample = incSample,
           pathNum = 15, R = 10, compareRef = T,
           convertSymbol = T, pathDb = "reactome", compareRefType = "intersection",
           expRow = "ENSEMBL", sizeDep = T, dep = dep, strengthPlot = T, nStrength = 10)
           
cl = makeCluster(8)
bngeneplot(results = pway,
           exp = vsted,
           expSample = incSample,
           pathNum = c(15, 16), R = 10,
           convertSymbol = T,
           expRow = "ENSEMBL")                    

image.png image.png

---

展望

在本期推文中，小编从GEO数据库上下载了示例数据病，并后续进行了差异分析和富集分析，随后演示了如何利用CBNplot来展示感兴趣通路中的基因之间的调控关系。但是这种调控关系只是CBNplot基于基因在各样本之间的表达水平进行的预测，并不能代表实际存在的调控关系。在实际应用中，还需要根据CHIPseq，ATAC等实验数据进一步证实某两个基因之间有之间的相互作用。无论如何，预测的结果可能不是很完美但总归比没有好；基于此，我们还可以根据相关的生物学知识和文献检索先建立几个假说，最后再使用实验进行验证，

好啦，本期推文到这就结束啦，在下期的推文中，Immugent将会介绍如何使用CBNplot在通路水平进行展示。