各位解螺旋的小伙伴們大家好，我是解螺旋的先鋒宇，今天我來(lái)給大家介紹一個(gè)強(qiáng)大的R包-GDCRNAtools包！這個(gè)包簡(jiǎn)直就是將差異分析，富集分析，火山圖繪制，生存分析以及ceRNA構(gòu)建強(qiáng)力整合，自從用上這個(gè)包，TCGA數(shù)據(jù)分析，可視化以及生信里面ceRNA的構(gòu)建簡(jiǎn)直省心得不能再省心了，既然這個(gè)包如此強(qiáng)大，相信大家已經(jīng)無(wú)法按耐住學(xué)習(xí)的沖動(dòng)，那就隨著小編的腳步一起來(lái)看看這個(gè)包到底強(qiáng)大在何處~
包的安裝
首先我們來(lái)安裝這個(gè)包，小編R的版本是4.0.2，所以我就直接使用bioconductor的安裝方式安裝了，然后library這個(gè)包
if (!requireNamespace('BiocManager', quietly = TRUE)) install.packages('BiocManager')BiocManager::install('GDCRNATools')library(GDCRNATools)差異分析
GDCRNAtools包的gdcDEAnalysis提供了三種做差異分析的方法，分別是limma，edgeR以及DESeq2,這三種方法可以說(shuō)是現(xiàn)在做差異分析最常用的方法，一個(gè)包都納入了，學(xué)會(huì)了這個(gè)包就直接相當(dāng)于學(xué)會(huì)三個(gè)包
首先我們來(lái)構(gòu)建一個(gè)表達(dá)矩陣和分組文件，這里小編只是為了節(jié)省跑代碼的時(shí)間，所以沒(méi)有使用TCGA的數(shù)據(jù)，可以看到我們構(gòu)建的數(shù)據(jù)是和TCGA下載的counts文件是類似的，所以小伙伴們可以直接下載TCGA的counts數(shù)據(jù)就可以進(jìn)行差異分析了
genes <- c('ENSG00000000938','ENSG00000000971','ENSG00000001036','ENSG00000001084','ENSG00000001167','ENSG00000001460')samples <- c('TCGA-2F-A9KO-01', 'TCGA-2F-A9KP-01','TCGA-2F-A9KQ-01', 'TCGA-2F-A9KR-11','TCGA-2F-A9KT-11', 'TCGA-2F-A9KW-11')metaMatrix <- data.frame(sample_type=rep(c('Tumor','Normal'),each=3),sample=samples)rnaMatrix <- matrix(c(6092,11652,5426,4383,3334,2656,8436,2547,7943,3741,6302,13976,1506,6467,5324,3651,1566,2780,834,4623,10275,5639,6183,4548,24702,43,1987,269,3322,2410,2815,2089,3804,230,883,5415), 6,6)rownames(rnaMatrix) <- genescolnames(rnaMatrix) <- samples
接下來(lái)我們就使用gdcDEAnalysis來(lái)進(jìn)行差異分析
DEGAll <- gdcDEAnalysis(counts = rnaMatrix, group = metaMatrix$sample_type, comparison = 'Tumor-Normal', method = 'limma')看到這里做生信分析的小伙伴肯定想著要是我早點(diǎn)遇到這個(gè)包，我就不用limma，DESeq2，edgeR要分別用三套代碼來(lái)跑了。
火山圖繪制
差異分析做完，這個(gè)時(shí)候來(lái)一個(gè)火山圖就是非常適合了，以前我們繪制火山圖還需要增加額外的一步，那就是根據(jù)logfoldchange的正負(fù)值得到基因表達(dá)是up還是down，但是在GDCRNAtools包中我們可以直接使用gdcVolcanoPlot函數(shù)一步搞定，只需要把差異分析的數(shù)據(jù)框直接放進(jìn)去就可以出圖。
為了讓大家能夠看公眾號(hào)就能看得很清楚明白，這里我們也可以構(gòu)建一個(gè)差異分析后的數(shù)據(jù)框：
genes <- c('ENSG00000231806','ENSG00000261211','ENSG00000260920','ENSG00000228594','ENSG00000125170','ENSG00000179909','ENSG00000280012','ENSG00000134612','ENSG00000213071')symbol <- c('PCAT7','AL031123.2','AL031985.3','FNDC10','DOK4','ZNF154','RPL23AP61','FOLH1B','LPAL2')group <- rep(c('long_non_coding','protein_coding','pseudogene'), each=3)logFC <- c(2.8,2.3,-1.1,1.9,-1.2,-1.6,1.5,2.1,-1.1)FDR <- rep(c(0.1,0.00001,0.0002), each=3)deg <- data.frame(symbol, group, logFC, FDR)rownames(deg) <- genes
構(gòu)建好后直接進(jìn)行火山圖的繪制
gdcVolcanoPlot(deg.all=deg)
當(dāng)然這里完全可以使用第一步得到的差異分析結(jié)果直接放進(jìn)這個(gè)函數(shù)，就可以得到火山圖
富集分析
接下來(lái)我們進(jìn)行差異分析，相信很多剛開(kāi)始做生信分析的小伙伴都會(huì)遇到一個(gè)難題，那就是ID轉(zhuǎn)換，很多時(shí)候我們想把TCGA差異的基因拿去看看下游能夠富集到什么通路，結(jié)果像clusterprofiler包都是不接受TCGA的基因編號(hào)的，需要進(jìn)行轉(zhuǎn)換，但是我們現(xiàn)在大可不必，我們可以使用GDCRNAtools包的gdcEnrichAnalysis函數(shù)直接進(jìn)行GO和KEGG,DO富集分析。
首先為了簡(jiǎn)化，我們來(lái)取幾個(gè)TCGA的基因來(lái)演示，大家自己在做的時(shí)候可以直接把相應(yīng)的gene那一列拿進(jìn)來(lái)就可以
deg <- c('ENSG00000000938','ENSG00000000971','ENSG00000001036','ENSG00000001084','ENSG00000001167','ENSG00000001460')
然后進(jìn)行富集分析(這里是GO,KEGG,DO富集分析一起做的)
enrichOutput <- gdcEnrichAnalysis(gene=deg, simplify=TRUE)但是如果我們要指定進(jìn)行KEGG進(jìn)行分析怎么辦呢，我們可以直接使用把上面的結(jié)果輸入到GDCRNAtools包中的gdcEnrichPlot函數(shù)：
gdcEnrichPlot(enrichment, type = 'bar', category = 'KEGG',num.terms = 10, bar.color = 'red')
當(dāng)然我們也可以在type參數(shù)下面設(shè)置其它的選項(xiàng)，比如GO以及GO的三個(gè)子類。
ceRNA構(gòu)建小工具
很多時(shí)候我們想構(gòu)建ceRNA的時(shí)候，很多小伙伴可能就是想著去在線數(shù)據(jù)庫(kù)搜索，但是檢索出來(lái)了很多l(xiāng)ncRNA，miRNA，mRNA的時(shí)候，應(yīng)該如何進(jìn)一步篩選呢，這個(gè)時(shí)候我們可以把對(duì)應(yīng)的TCGA數(shù)據(jù)提取出來(lái)，然后看它們表達(dá)的相關(guān)性情況。
比如在這里我們繼續(xù)使用構(gòu)建一個(gè)類似TCGA數(shù)據(jù)結(jié)構(gòu)的數(shù)據(jù)框
genes <- c('ENSG00000000938','ENSG00000000971','ENSG00000001036', 'ENSG00000001084','ENSG00000001167','ENSG00000001460')samples <- c('TCGA-2F-A9KO-01', 'TCGA-2F-A9KP-01', 'TCGA-2F-A9KQ-01', 'TCGA-2F-A9KR-01', 'TCGA-2F-A9KT-01', 'TCGA-2F-A9KW-01')metaMatrix <- data.frame(sample_type=rep('PrimaryTumor',6), sample=samples)rnaExpr <- matrix(c(2.7,7.0,4.9,6.9,4.6,2.5, 0.5,2.5,5.7,6.5,4.9,3.8, 2.1,2.9,5.9,5.7,4.5,3.5, 2.7,5.9,4.5,5.8,5.2,3.0, 2.5,2.2,5.3,4.4,4.4,2.9, 2.4,3.8,6.2,3.8,3.8,4.2),6,6)rownames(rnaExpr) <- genescolnames(rnaExpr) <- samples然后我們使用shinyCorPlot函數(shù)，就可以自己選擇看哪些基因和哪些基因之間在TCGA的表達(dá)相關(guān)性
生存曲線繪制
很多時(shí)候我們畫(huà)某個(gè)基因高低表達(dá)對(duì)應(yīng)的生存曲線的時(shí)候，首先想到的是要構(gòu)建一個(gè)生存函數(shù)，然后在ggsurplot函數(shù)里面繪制，然后還要加上p值等等，但是GDCRNAtools包的gdcKMPlot函數(shù)把這些都整合了在一起，讓我們一個(gè)函數(shù)直接搞定，對(duì)于生信還是起步的階段的小伙伴簡(jiǎn)直太貼心啦~
首先我們也使用上面的數(shù)據(jù)，只不過(guò)在metaMatrix里面加入生存數(shù)據(jù)
genes <- c('ENSG00000000938','ENSG00000000971','ENSG00000001036','ENSG00000001084','ENSG00000001167','ENSG00000001460')samples <- c('TCGA-2F-A9KO-01', 'TCGA-2F-A9KP-01','TCGA-2F-A9KQ-01', 'TCGA-2F-A9KR-01','TCGA-2F-A9KT-01', 'TCGA-2F-A9KW-01')metaMatrix <- data.frame(sample_type=rep('PrimaryTumor',6),sample=samples,days_to_death=seq(100,600,100))rnaExpr <- matrix(c(2.7,7.0,4.9,6.9,4.6,2.5,0.5,2.5,5.7,6.5,4.9,3.8,2.1,2.9,5.9,5.7,4.5,3.5,2.7,5.9,4.5,5.8,5.2,3.0,2.5,2.2,5.3,4.4,4.4,2.9,2.4,3.8,6.2,3.8,3.8,4.2),6,6)rownames(rnaExpr) <- genescolnames(rnaExpr) <- samples
然后我們使用gdcKMPlot函數(shù)進(jìn)行繪圖
gdcKMPlot(gene='ENSG00000001167', rna.expr=rnaExpr, metadata=metaMatrix, sep='median')
這樣一張樸素簡(jiǎn)約但是又拿得出手的生存曲線就繪制好了！
看到這里相信大家已經(jīng)摩拳擦掌，準(zhǔn)備在自己的TCGA數(shù)據(jù)里面躍躍欲試?yán)玻蔷托袆?dòng)起來(lái)，看能不能分析出一點(diǎn)有意思的東西，最后歡迎大家繼續(xù)關(guān)注挑圈聯(lián)靠公眾號(hào)，小編會(huì)在以后的推文中給大家介紹更多生產(chǎn)力包，從而大大提高大家玩生信的效率！
歡迎大家關(guān)注解螺旋生信頻道-挑圈聯(lián)靠公號(hào)~
聯(lián)系客服
