为什么我一行代码就可以完成3个R包的RNA-seq差异分析呢

Posted on by

在教师节收到学生的提问，刷我B站74小时视频的时候看到我演示了RNA-seq差异分析只用了一行代码就完成了3大R包的全部分析，并且输出了对应的图表结果，觉得很神奇，但是B站视频并没有配套讲义和代码还有测试数据。

首先我一直使用airway数据集做测试

airway数据集在这里我就不多说了，搜索生信技能树早期教程可以看到很多介绍。

## 表达矩阵来自于R包： airway
if(F){
 library(airway)
 data(airway)
 exprSet=assay(airway)
 group_list=colData(airway)[,3]
 save(exprSet,group_list,file = 'airway_exprSet.Rdata')
}

load(file = 'airway_exprSet.Rdata')

if(T){
 colnames(exprSet)
 pheatmap::pheatmap(cor(exprSet))
 group_list
 tmp=data.frame(g=group_list)
 rownames(tmp)=colnames(exprSet)
 # 组内的样本的相似性应该是要高于组间的！
 pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
 dim(exprSet)
 exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
 dim(exprSet)

 exprSet=log(edgeR::cpm(exprSet)+1)
 dim(exprSet)
 exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
 dim(exprSet)
 M=cor(log2(exprSet+1))
 tmp=data.frame(g=group_list)
 rownames(tmp)=colnames(M)
 pheatmap::pheatmap(M,annotation_col = tmp)
 pheatmap::pheatmap(M,annotation_col = tmp,filename = 'cor.png')

library(pheatmap)
 pheatmap(scale(cor(log2(exprSet+1))))

}

很明显可以看到， 组内的样本的相似性应该是要高于组间的！

而且为了显示这个规律，我还做了一个统计学技巧展示，当然了，很多人非常的不用心，所以把视频听10遍也看不懂，get不到我的点，需要批评！

使用我包装好的函数即可

可以看到，下面的代码非常简洁，因为仅仅是使用了 run_DEG_RNAseq 函数，就根据表达矩阵和分组信息，完成了全部的分析！

rm(list = ls())
options(stringsAsFactors = F)
load(file = 'airway_exprSet.Rdata')
group_list
group_list=relevel(group_list,ref = 'untrt')
source('run_DEG_RNA-seq.R')
run_DEG_RNAseq(exprSet,group_list,
 g1="untrt",g2="trt",
 pro='airway')

这就是大家看视频后提的问题，为什么这么神奇呢？下面的图表是如何自动出来的呢？

因为这个 run_DEG_RNAseq 函数的代码非常长，这里我就不贴在公众号了哈，大家可以在我的GitHub的GEO项目找到它！

GEO传奇代码

一不留神，这个GEO项目就成为了点赞数最多的，直接孵化出12篇数据挖掘类SCI文章，至于间接的那些就不计其数了，因为大家都是偷偷的使用，也不告诉我，甚至某些别有用心者还不告诉身边的人，要一个人独享这些代码。

https://github.com/jmzeng1314/GEO

https://github.com/jmzeng1314/GEO/tree/master/airway_RNAseq

既然是多个R包，结果该如何取舍呢？

这个时候是没有标准答案的，因为每个R包都非常热门，引用量都是好几千，你选择哪个都符合市场规律，不过，我这里有一个代码，对3个结果根据阈值筛选交集。

https://github.com/jmzeng1314/GEO/tree/master/airway_RNAseq

差异基因后是不是也可以批量GO/KEGG数据库注释呢？

当然是啊，都会写代码了，还有什么是不能为所欲为的呢？

同样的，代码也是在GitHub，需要你仔细理解，不过我有一个小小的要求，请不要把我的代码雪藏，或者刻意隐瞒。

https://github.com/jmzeng1314/GEO/tree/master/airway_RNAseq

值得一提的是这里面的一行代码是需要格外注意的哦：

group_list=relevel(group_list,ref = 'untrt')
​

Comments are closed.