Chord Diagram弦狀圖-用來(lái)畫(huà)富集分析的結(jié)果
前面看的文獻(xiàn)中，有一個(gè)很漂亮的圈圖，https://mp.weixin.qq.com/s/NIqSTvJ916ZGk_dwsRlH5w，就是kegg和go分析后的基因和通路之間的關(guān)系圖，用cnetplot函數(shù)畫(huà)出來(lái)的圖也是一樣的含義，不過(guò)這個(gè)圈圖看起來(lái)就真的很好看，就想畫(huà)，主要其實(shí)顏色是可以選擇的，用https://www.sioe.cn/yingyong/yanse-rgb-16/中的顏色就可以選擇，所以可以畫(huà)出來(lái)更好看的圖。不過(guò)后面在數(shù)據(jù)變換的時(shí)候，就還是沒(méi)轉(zhuǎn)過(guò)來(lái)可以直接用最后面畫(huà)圖用的數(shù)據(jù)，就還在思考如何一步步構(gòu)建包中的中間數(shù)據(jù)從而得到足后的數(shù)據(jù)，后來(lái)老大告訴我直接構(gòu)建后面的數(shù)據(jù)就可以了。
弦狀圖這種類(lèi)型的圖表可視化了實(shí)體之間的相互關(guān)系。實(shí)體之間的連接用于顯示它們共享某些相同的東西。節(jié)點(diǎn)沿圓排列，通過(guò)使用圓弧將點(diǎn)之間的關(guān)系連接在一起。值被分配給每個(gè)連接，它由每個(gè)弧的大小成比例地表示。顏色可以用來(lái)把數(shù)據(jù)分組到不同的類(lèi)別，這有助于進(jìn)行比較和區(qū)分組。不過(guò)當(dāng)有太多的連接顯示時(shí)，可能就會(huì)出現(xiàn)混亂。
這使得弦圖非常適合于比較數(shù)據(jù)集內(nèi)或不同數(shù)據(jù)組之間的相似性。
使用R語(yǔ)言中的繪圖包，比如 circlize, RCircos, CIRCUS, and OmicCircos等來(lái)繪制。
一張顯示人口流動(dòng)的圖如下，來(lái)自http://download.gsb.bund.de/BIB/global_flow/
GOplot包可用于生物數(shù)據(jù)的可視化。但是要注意該包不能用于執(zhí)行這些分析，只能把分析結(jié)果進(jìn)行可視化。
這個(gè)包的用法參考：http://wencke.github.io/。主要是可以把富集分析后的結(jié)果進(jìn)行圈圖展示。鏈接中的講解代碼很長(zhǎng)，提煉后的代碼如下
david <- EC$david
genelist <- EC$genelist
circ <- circle_dat(EC$david, EC$genelist)
chord <- chord_dat(data = circ, genes = EC$genes, process = EC$process)
GOChord(chord, space = 0.02, gene.order = 'logFC', gene.space = 0.25, gene.size = 5)
畫(huà)出的圖如下
EC是這個(gè)GOplot包中的實(shí)例數(shù)據(jù)，為了能用這個(gè)包中的作圖函數(shù)，就先了解這個(gè)包中的數(shù)據(jù)。那么david和genelist的數(shù)據(jù)長(zhǎng)下面這樣。
生成的中間數(shù)據(jù)chord長(zhǎng)下面這樣。
最后用GOChord函數(shù)畫(huà)圖的數(shù)據(jù)chord的數(shù)據(jù)格式如下圖所示。這是我們要把自己的數(shù)據(jù)變成這樣的數(shù)據(jù)的樣子。就是下面的chord數(shù)據(jù)格式。變成下面這樣就可以畫(huà)圈圖了。
上面最后畫(huà)圖的輸入數(shù)據(jù)函數(shù)就是chord，最后的logFC是從前面的差異分析的結(jié)果得來(lái)的，也就是genelist中的logFC。然后chord數(shù)據(jù)除了logFC以外，前面的數(shù)據(jù)就是如果某個(gè)基因在通路中，就是1，否則就是0。這中數(shù)據(jù)就是最前面代碼中的chord_dat函數(shù)實(shí)現(xiàn)的，但是問(wèn)題是這個(gè)chord 數(shù)據(jù)是從circ得到的，而circ是從david和genelist得到的，david數(shù)據(jù)是富集分析的結(jié)果，genelist是差異分析的結(jié)果，現(xiàn)在我沒(méi)有差異分析的結(jié)果，我只是有前面比如通過(guò)cox回歸篩選得到的基因名，怎么辦呢？
問(wèn)題解決的辦法就是，直接制作數(shù)據(jù)格式到chord這樣的數(shù)據(jù)，把logFC變?yōu)?就可以了。
- 下面的代碼是從前面cox回歸篩選到的基因開(kāi)始的
load(file = 'sur_logKM_cox.Rdata')
cox=rownames(cox_results[cox_results[,4]<0.01,])
log=names(log_rank_p[log_rank_p<0.01])
surGenes=intersect(cox,log)
上面得到的surGenes如下：就是一堆字符串啦
但是用老大出品的這個(gè)AnnoProbe進(jìn)行ID轉(zhuǎn)換，得到了如下圖所示非常方便后續(xù)進(jìn)行操作的表格
library(AnnoProbe) #用老大出品的這個(gè)AnnoProbe進(jìn)行ID轉(zhuǎn)換，
surGenes=annoGene(surGenes,ID_type = 'SYMBOL','human')
tail(sort(table(surGenes$biotypes)))
head(surGenes)
接下來(lái)富集分析
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
df <- bitr(unique(surGenes$SYMBOL), fromType = "SYMBOL",
toType = c( "ENTREZID"),
OrgDb = org.Hs.eg.db)#富集分析要求數(shù)據(jù)數(shù)據(jù)的是ENTREZID,要從前面得到的SYMBOL用bitr函數(shù)得到ENTREZID
head(df)
gene_up=df$ENTREZID
#下面是為了畫(huà)圖做的插曲
#將bp的結(jié)果保存
bp_go <- enrichGO(gene_up,
OrgDb = "org.Hs.eg.db",
ont="BP",
pvalueCutoff = 0.001,
qvalueCutoff = 0.001)
enrichgo=DOSE::setReadable(bp_go, OrgDb='org.Hs.eg.db',keyType='ENTREZID')#setReadable函數(shù)直接進(jìn)行ID轉(zhuǎn)換了，從ENTREZID再得到symbol
mm <- enrichgo@result[1:5,c(2,7,8)]
得到的mm數(shù)據(jù)如下（mm是我名字的縮寫(xiě)，哈哈）
tmp=do.call(rbind,
apply(mm, 1,function(x){
data.frame(go=x[1],
gene=strsplit(x[3],'/')[[1]])
})
)
得到的tmp數(shù)據(jù)格式如下
接下來(lái)就是將這個(gè)長(zhǎng)形數(shù)據(jù)用reshape包中的dcast函數(shù)變?yōu)閷捫螖?shù)據(jù)，如下
library(reshape2)
tmp2=dcast(tmp,go~gene)
tmp2[is.na(tmp2)]=0
rownames(tmp2)=tmp2[,1]
tmp2=tmp2[,-1]
tmp2=t(tmp2)
tmp2[tmp2!=0]=1
tmp2=as.data.frame(tmp2)
tmp2$logFC=0
cg=rownames(tmp2)
tmp2=apply(tmp2,2,as.numeric)
rownames(tmp2)=cg
得到的tmp2數(shù)據(jù)如下，這就是已經(jīng)做好了和前面的實(shí)例數(shù)據(jù)chord長(zhǎng)得一摸一樣的數(shù)據(jù)了。
接下來(lái)畫(huà)圖
GOChord(tmp2, space = 0.02, gene.order = 'logFC', gene.space = 0.25, gene.size = 5)
區(qū)別是我們沒(méi)有用變化倍數(shù)進(jìn)行聚類(lèi)，所以沒(méi)有l(wèi)ogFC的漸變的變化，不過(guò)已經(jīng)超級(jí)好看了！
如果想有倍數(shù)變化的結(jié)果，那么就是一開(kāi)始就是有差異分析的結(jié)果，就可以下面的代碼，得到的就是本文第二張圖片了。
# Create the plot
chord <- chord_dat(data = circ, genes = EC$genes, process = EC$process)
GOChord(chord, space = 0.02, gene.order = 'logFC', gene.space = 0.25, gene.size = 5)
參考的代碼全套在http://wencke.github.io/，也可以搜索到翻譯過(guò)來(lái)的中文版的教程。
