一站式分析R包来了，承包了生信各种分析！太全能了！

一文学会“多面手”GDCRNAtools包用法

各位解螺旋的小伙伴们大家好，我是解螺旋的先锋宇，今天我来给大家介绍一个强大的R包-GDCRNAtools包！这个包简直就是将差异分析，富集分析，火山图绘制，生存分析以及ceRNA构建强力整合，自从用上这个包，TCGA数据分析，可视化以及生信里面ceRNA的构建简直省心得不能再省心了，既然这个包如此强大，相信大家已经无法按耐住学习的冲动，那就随着小编的脚步一起来看看这个包到底强大在何处~

包的安装

首先我们来安装这个包，小编R的版本是4.0.2，所以我就直接使用bioconductor的安装方式安装了，然后library这个包

if (!requireNamespace('BiocManager', quietly = TRUE)) install.packages('BiocManager')BiocManager::install('GDCRNATools')library(GDCRNATools)

差异分析

GDCRNAtools包的gdcDEAnalysis提供了三种做差异分析的方法，分别是limma，edgeR以及DESeq2,这三种方法可以说是现在做差异分析最常用的方法，一个包都纳入了，学会了这个包就直接相当于学会三个包

首先我们来构建一个表达矩阵和分组文件，这里小编只是为了节省跑代码的时间，所以没有使用TCGA的数据，可以看到我们构建的数据是和TCGA下载的counts文件是类似的，所以小伙伴们可以直接下载TCGA的counts数据就可以进行差异分析了

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要分别用三套代码来跑了。

火山图绘制

差异分析做完，这个时候来一个火山图就是非常适合了，以前我们绘制火山图还需要增加额外的一步，那就是根据logfoldchange的正负值得到基因表达是up还是down，但是在GDCRNAtools包中我们可以直接使用gdcVolcanoPlot函数一步搞定，只需要把差异分析的数据框直接放进去就可以出图。

为了让大家能够看公众号就能看得很清楚明白，这里我们也可以构建一个差异分析后的数据框：

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)

当然这里完全可以使用第一步得到的差异分析结果直接放进这个函数，就可以得到火山图

富集分析

接下来我们进行差异分析，相信很多刚开始做生信分析的小伙伴都会遇到一个难题，那就是ID转换，很多时候我们想把TCGA差异的基因拿去看看下游能够富集到什么通路，结果像clusterprofiler包都是不接受TCGA的基因编号的，需要进行转换，但是我们现在大可不必，我们可以使用GDCRNAtools包的gdcEnrichAnalysis函数直接进行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函数：

gdcEnrichPlot(enrichment, type = 'bar', category = 'KEGG',num.terms = 10, bar.color = 'red')

当然我们也可以在type参数下面设置其它的选项，比如GO以及GO的三个子类。

ceRNA构建小工具

很多时候我们想构建ceRNA的时候，很多小伙伴可能就是想着去在线数据库搜索，但是检索出来了很多lncRNA，miRNA，mRNA的时候，应该如何进一步筛选呢，这个时候我们可以把对应的TCGA数据提取出来，然后看它们表达的相关性情况。

比如在这里我们继续使用构建一个类似TCGA数据结构的数据框

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函数，就可以自己选择看哪些基因和哪些基因之间在TCGA的表达相关性

生存曲线绘制

很多时候我们画某个基因高低表达对应的生存曲线的时候，首先想到的是要构建一个生存函数，然后在ggsurplot函数里面绘制，然后还要加上p值等等，但是GDCRNAtools包的gdcKMPlot函数把这些都整合了在一起，让我们一个函数直接搞定，对于生信还是起步的阶段的小伙伴简直太贴心啦~

首先我们也使用上面的数据，只不过在metaMatrix里面加入生存数据

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函数进行绘图

gdcKMPlot(gene='ENSG00000001167', rna.expr=rnaExpr, metadata=metaMatrix, sep='median')

这样一张朴素简约但是又拿得出手的生存曲线就绘制好了！

看到这里相信大家已经摩拳擦掌，准备在自己的TCGA数据里面跃跃欲试啦，那就行动起来，看能不能分析出一点有意思的东西，最后欢迎大家继续关注挑圈联靠公众号，小编会在以后的推文中给大家介绍更多生产力包，从而大大提高大家玩生信的效率！

欢迎大家关注解螺旋生信频道-挑圈联靠公号~