各位解螺旋的小伙伴們大家好，我是解螺旋的先鋒宇，今天我來給大家介紹一個強大的R包-GDCRNAtools包！這個包簡直就是將差異分析，富集分析，火山圖繪制，生存分析以及ceRNA構建強力整合，自從用上這個包，TCGA數(shù)據(jù)分析，可視化以及生信里面ceRNA的構建簡直省心得不能再省心了，既然這個包如此強大，相信大家已經(jīng)無法按耐住學習的沖動，那就隨著小編的腳步一起來看看這個包到底強大在何處~
包的安裝
首先我們來安裝這個包，小編R的版本是4.0.2，所以我就直接使用bioconductor的安裝方式安裝了，然后library這個包
if (!requireNamespace('BiocManager', quietly = TRUE)) install.packages('BiocManager')BiocManager::install('GDCRNATools')library(GDCRNATools)差異分析
GDCRNAtools包的gdcDEAnalysis提供了三種做差異分析的方法，分別是limma，edgeR以及DESeq2,這三種方法可以說是現(xiàn)在做差異分析最常用的方法，一個包都納入了，學會了這個包就直接相當于學會三個包
首先我們來構建一個表達矩陣和分組文件，這里小編只是為了節(jié)省跑代碼的時間，所以沒有使用TCGA的數(shù)據(jù)，可以看到我們構建的數(shù)據(jù)是和TCGA下載的counts文件是類似的，所以小伙伴們可以直接下載TCGA的counts數(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-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
接下來我們就使用gdcDEAnalysis來進行差異分析
DEGAll <- gdcDEAnalysis(counts = rnaMatrix, group = metaMatrix$sample_type, comparison = 'Tumor-Normal', method = 'limma')看到這里做生信分析的小伙伴肯定想著要是我早點遇到這個包，我就不用limma，DESeq2，edgeR要分別用三套代碼來跑了。
火山圖繪制
差異分析做完，這個時候來一個火山圖就是非常適合了，以前我們繪制火山圖還需要增加額外的一步，那就是根據(jù)logfoldchange的正負值得到基因表達是up還是down，但是在GDCRNAtools包中我們可以直接使用gdcVolcanoPlot函數(shù)一步搞定，只需要把差異分析的數(shù)據(jù)框直接放進去就可以出圖。
為了讓大家能夠看公眾號就能看得很清楚明白，這里我們也可以構建一個差異分析后的數(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
構建好后直接進行火山圖的繪制
gdcVolcanoPlot(deg.all=deg)
當然這里完全可以使用第一步得到的差異分析結果直接放進這個函數(shù)，就可以得到火山圖
富集分析
接下來我們進行差異分析，相信很多剛開始做生信分析的小伙伴都會遇到一個難題，那就是ID轉(zhuǎn)換，很多時候我們想把TCGA差異的基因拿去看看下游能夠富集到什么通路，結果像clusterprofiler包都是不接受TCGA的基因編號的，需要進行轉(zhuǎn)換，但是我們現(xiàn)在大可不必，我們可以使用GDCRNAtools包的gdcEnrichAnalysis函數(shù)直接進行GO和KEGG,DO富集分析。
首先為了簡化，我們來取幾個TCGA的基因來演示，大家自己在做的時候可以直接把相應的gene那一列拿進來就可以
deg <- c('ENSG00000000938','ENSG00000000971','ENSG00000001036','ENSG00000001084','ENSG00000001167','ENSG00000001460')
然后進行富集分析(這里是GO,KEGG,DO富集分析一起做的)
enrichOutput <- gdcEnrichAnalysis(gene=deg, simplify=TRUE)但是如果我們要指定進行KEGG進行分析怎么辦呢，我們可以直接使用把上面的結果輸入到GDCRNAtools包中的gdcEnrichPlot函數(shù)：
gdcEnrichPlot(enrichment, type = 'bar', category = 'KEGG',num.terms = 10, bar.color = 'red')
當然我們也可以在type參數(shù)下面設置其它的選項，比如GO以及GO的三個子類。
ceRNA構建小工具
很多時候我們想構建ceRNA的時候，很多小伙伴可能就是想著去在線數(shù)據(jù)庫搜索，但是檢索出來了很多l(xiāng)ncRNA，miRNA，mRNA的時候，應該如何進一步篩選呢，這個時候我們可以把對應的TCGA數(shù)據(jù)提取出來，然后看它們表達的相關性情況。
比如在這里我們繼續(xù)使用構建一個類似TCGA數(shù)據(jù)結構的數(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的表達相關性
生存曲線繪制
很多時候我們畫某個基因高低表達對應的生存曲線的時候，首先想到的是要構建一個生存函數(shù)，然后在ggsurplot函數(shù)里面繪制，然后還要加上p值等等，但是GDCRNAtools包的gdcKMPlot函數(shù)把這些都整合了在一起，讓我們一個函數(shù)直接搞定，對于生信還是起步的階段的小伙伴簡直太貼心啦~
首先我們也使用上面的數(shù)據(jù)，只不過在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ù)進行繪圖
gdcKMPlot(gene='ENSG00000001167', rna.expr=rnaExpr, metadata=metaMatrix, sep='median')
這樣一張樸素簡約但是又拿得出手的生存曲線就繪制好了！
看到這里相信大家已經(jīng)摩拳擦掌，準備在自己的TCGA數(shù)據(jù)里面躍躍欲試啦，那就行動起來，看能不能分析出一點有意思的東西，最后歡迎大家繼續(xù)關注挑圈聯(lián)靠公眾號，小編會在以后的推文中給大家介紹更多生產(chǎn)力包，從而大大提高大家玩生信的效率！
歡迎大家關注解螺旋生信頻道-挑圈聯(lián)靠公號~
