前天，我们已经给大家分享了GEO芯片测序limma差异分析非配对的万能代码（点击查看下载）。昨天，我们已经给大家分享了GEO芯片测序limma差异分析配对的万能代码（点击查看下载）。那我们从GEO下载的或者我们自己测序的高通量测序（RNAseq）矩阵（Raw count。从下图我们可以看到公司给我们的一般是Raw count，或者FPKM/RPKM count，这个矩阵我们在发文章的时候一般需要上传到GEO或者TCGA）：![]()
因此，我们在挖掘GEO数据的时候，要么是Raw count矩阵，要么是FPKM/RPKM count矩阵。Raw count矩阵一般用Deseq2包进行差异分析。FPKM/RPKM count矩阵一般用limma包差异分析。今天我们给大家分享的是Raw count的Deseq2包差异分析万能代码：“Raw count的Deseq2包差异分析+热图+火山图+GO分析+KEGG分析+GSEA分析全套万能代码（附代码交流群）”
🔽①长按下方二维码关注🔽
![]()
②对话框输入关键词：全套代码
②对话框输入关键词：全套代码
②对话框输入关键词：全套代码
是全套的万能代码，自动判断同基因名处理，你只要知道下载的数据有几个样本就行。差异分析+热图+火山图+GO分析+KEGG分析+GSEA分析，万能代码一键搞定。大家可以进群，群内不答疑，可代码共享。下面开始操作：DESeq2是一个比较常用的转录组分析R包，包的使用非常简单，与之前的limma包不一样，DESeq2需要的数据是Row counts矩阵，这点非常重要。所以不管你自己的测序数据，还是GEO下载的数据，你需要明确是不是Row counts矩阵。比如我从GEO下载GSE169758的测序数据（我这里只是举例子），如下操作，打开生物医学之家网站（swyxzj.com）进入GEO官网：![]()
检索GSE号，并且明确知道是测序数据，不是芯片数据：
![]()
![]()
下载、解压，整理成下面（前面三个对照组，后面3个实验组）：
![]()
将raw count矩阵和代码放在同一个文件夹下（RNAseq_Deseq2）：![]()
### 一、DEseq2差异分析+标准化PCA图### 二、ENSG转Symbol----公众号：科研部，参考链接：https://mp.weixin.qq.com/s/l3jgf_GhglZ5pwUNhc1x1g### 三、输出差异基因### 四、作图，差异表达图，热图，火山图### 五、差异基因集的GO分析### 六、差异基因集的KEGG分析### 七、 所有基因集的GSEA富集分析
![]()
如上代码，是3个对照组和3个实验组，所以你如果你的数据是7个，你就修改为7就行。之后，你会生成一个PCA图（感觉分组分的不好，这个数据有点问题）：
![]()
同时，生成DEseq2差异分析结果文件：
![]()
那我们根据PCA图可以看到有个别样本不行，错综复杂。意思就是这个数据不太好。那我们再换一个我自己的数据重新差异分析看看，数据如下：
![]()
读入数据，DEseq2差异分析（唯一要改的就是7个对照组，几个实验组）：
![]()
![]()
分的很开，结果可靠。差异分析结果文件如下（但是是ENSG的，我们后续做GO、kegg、gsea都是需要基因symbol，所以下一步，我们要转换）：![]()
ENSG转symbol，请仔细看这篇推文（点击查看）。我们测序之后，比对到基因组文件的时候是用的哪个版本的GTF文件，你就需要下载该版本的GTF文件，用ENSGmap进去，得到基因Symbol，那我这里是自己测序的数据，用的是gencode.v39.文件，所以就下载1_gencode.v39.annotation.gtf.gz（怎么下载，请点击查看），放在这个文件夹下：![]()
运行代码（你需要自己下载gtf文件，修改8：21的21，这个代码意思就是说取根据8-21列取每行的中位数，后面好比对到gtf文件时，同名的symbol处理，21的意思是7+n，n是样本数，我是14个样本，所以7+14=21）：![]()
就得到了symbol的DEseq2分析结果文件（有symbol了）：![]()
![]()
我这里取的阈值是logFC为1，矫正后的p值小于0.05，也就是变化倍数大于2，并且有统计学意义的基因作为差异基因，用于后续的分析，打开如下：![]()
我们仔细看看后会发现，我们第二列里面有些也是ENSG号。
![]()
![]()
发现有些ENSG基因，它的名字就是ENSG号，数量不多，不影响的。
首先就是做热图，我们不可能将所有的差异基因都展示出来，所以这里筛选出padj < 0.05和abs(log2FoldChange) >2的基因拿出来做热图展示，就是说明有部分基因变化非常明显：![]()
接着我们画火山图，说明很多基因确实发生了变化（蓝色和红色的基因有很多）：
![]()
我们要先加载差异基因文件，然后将其转化为entrezID：
![]()
![]()
这样我们后面用entrezID做GO分析的时候就会报错：错误: near "7006": syntax error。如下：![]()
这个时候，我们就可以随机选取该基因的任何一个entrezID都行。
![]()
![]()
这里是将P小于0.05的结果保留下来（万能代码，这一步，你什么都不需要改，后面的KEGG分析也是一样，不需要改）：
![]()
![]()
先输入差异基因，随后kegg富集分析，并保存P小于0.05的结果：
![]()
再选择P小于0.05的前10作图，并将KEGGmap到的通路下载下来：![]()
![]()
![]()
是用所有基因来做GSEA富集分析，不是用差异基因哦，并且按照logFC排序进行GSEA富集分析，运行到下面代码：
![]()
![]()
![]()
每一张都保存起来就行了，就等于一篇论文（你的数据有多新，你就能发多大的文章）：
![]()
![]()
全套代码一共是300多行，自己多练练就会了，不会写不要紧，看得懂就行了，每一句代码都有注释。代码和演示数据我都上传了：![]()
如下免费获取：
“Raw count的Deseq2包差异分析+热图+火山图+GO分析+KEGG分析+GSEA分析全套万能代码（附代码交流群）”🔽①长按下方二维码关注🔽
![]()
②对话框输入关键词：全套代码
②对话框输入关键词：全套代码
②对话框输入关键词：全套代码
