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

在教师节收到学生提问，刷我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项目找到它！

https://github.com/jmzeng1314/GEO/blob/master/airway_RNAseq/run_DEG_RNA-seq.R

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')

后记

本来到这里应该是要贴上我们全国巡讲的宣传，但是中秋节临时加开的广州特别班，我们没想过会有很多人报名，毕竟只有区区十天不到的报名时间，但是很快就满20人了，而且非常多迫切想学习的小伙伴找到我们，即使拒绝后仍然是发生如下的对话：

因为加人就要换大一点会议室，成本就增加几千块钱，所以只能是看有没有3个以上的人报名，至少把成本cover掉！

大家仍然是可以尝试报名，广州今年就这一场了，还有很多其他城市的粉丝嗷嗷待哺者等着我们！

号外：中秋节广州3天入门课程报名马上截止：（中秋节一起来学习！）全国巡讲第16站-广州（生信入门课加量不加价）