背景介紹
通常我們做完轉錄組學，進行一系列富集分析如GO、KEGG、GSEA等，但是這些結果無非就是富集出來術語以及對應的基因集合，或者通路及其對應的基因集合。
而現在，有一個R包，能夠基於富集分析的結果，從表達矩陣數據中通過繪製貝葉斯網絡圖，來進行推斷基因、通路的上下游分析，這樣，大大方便了我們進行對應的實驗驗證，也減少了實驗流程的主觀性，增強其客觀性和可信性。
那麼接下來，我們根據CBNplot的官方文檔，進行解釋說明該包的使用方法！
官方文檔地址：https://noriakis.github.io/software/CBNplot/
數據來源
數據來源於PubMed中的公開數據 GSE133624[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133624]，這是一份人膀胱尿路上皮癌及鄰近正常組織中基因表達的轉錄組學分析數據。
直接點擊矩陣文件即可用於示例
使用教程
1.在演示之前，首先要做的是進行GSE133624的處理和差異基因的確定。首先需要加載R包，然後讀取示例數據、數據過濾、差異分析、變量轉換、PCA分析
# 對於以下包可以使用下列方式安裝，以clusterProfiler舉例# if (!require("BiocManager", quietly = TRUE))# install.packages("BiocManager")library(BiocManager)# BiocManager::install("clusterProfiler")library(DESeq2)library(ggplot2)library(clusterProfiler)library(org.Hs.eg.db)library(ReactomePA)library(enrichplot)# CBNplot包安裝使用下列方式# install_github("noriakis/CBNplot",force = TRUE)library(CBNplot)## 首先加載數據集和做元數據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")meta$Condition <- factor(meta$Condition)dds <- DESeqDataSetFromMatrix(countData = counts, colData = meta, design= ~ Condition)## 進行數據的過濾filt <- rowSums(counts(dds) < 10) > dim(meta)[1]*0.9dds <- dds[!filt,]## 進行基因的差異分析dds = DESeq(dds)res = results(dds, pAdjustMethod = "bonferroni")## 對變量進行轉換v = vst(dds, blind=FALSE)vsted = assay(v)## 繪製PCA圖plotPCA(v, intgroup="Condition")+ theme_bw()
2.通過padj篩選結果，轉換基因ID進行富集分析，將geneID映射到基因Symbol，定義包括的樣本
## 定義輸入基因，使用clusterProfiler轉換基因ID.sig = subset(res, padj<0.05)cand.entrez = bitr(rownames(sig), fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Hs.eg.db)$ENTREZID## 進行富集分析pway = enrichPathway(gene = cand.entrez)pwayGO = enrichGO(cand.entrez, ont = "BP", OrgDb = org.Hs.eg.db)## 轉換為SYMBOLpway = setReadable(pway, OrgDb=org.Hs.eg.db)pwayGO = setReadable(pwayGO, OrgDb=org.Hs.eg.db)## 存儲相似性pway = pairwise_termsim(pway)## 定義包括的樣本incSample = rownames(subset(meta, Condition=="T"))## 對於基因集富集分析，需要獲取log2FC數值# 進行基因ID的轉換allEntrez = bitr(rownames(res), fromType="ENSEMBL", toType="ENTREZID", OrgDb=org.Hs.eg.db)# 構造基因ID-Log2FC列表res$ENSEMBL <- rownames(res)lfc <- merge(data.frame(res), allEntrez, by="ENSEMBL")lfc <- lfc[order(lfc$log2FoldChange, decreasing=TRUE),]geneList <- lfc$log2FoldChangenames(geneList) <- lfc$ENTREZID#進行ges富集分析pwayGSE <- gsePathway(geneList)# 通路中的平均基因數量，用於推斷sigpway <- subset(pway@result, p.adjust<0.05)paste(mean(sigpway$Count), sd(sigpway$Count))
3.畫個圖試試，該圖指的是富集出來的第17條通路上基因之間的調控關係或上下游關係
bngeneplot(results = pway, exp = vsted, pathNum = 17)
4.下來我們有空再解讀解讀其他圖譜的含義！
參考文獻
留言列表