截止2022年4月6日，GEO数据库共包含61,987个高通量测序数据以及65,021个芯片数据，为生信分析、临床研究以及实验设计提供了强有力的数据支撑，今天，主要为大家带来的是GEO数据库的差异分析操作。

图片

一、分析思路选择

图片

二、FAQ

2.1 如何判断双通道芯片（例：GSE111471）

1. 双通道芯片的sample title一般为A vs B

图片

2. 双通道芯片的GSM页面会有两个channel（关于什么是GSM以及GSM在哪，可以查看四文搞定GEO数据库转录组差异分析之简介）

图片

2.2 如何判断单通道芯片（例：GSE182065）

不是双通道芯片，就是单通道芯片 图片

2.3 为什么双通道芯片表达量求平均就可以了

如下图所示，双通道芯片原理为两个样本同时测（单通道为一个样一个表达量），其表达量本就是Cy5/Cy3，也就是实验组与对照组的比值（取log2），因此，一般将重复样本之间测量的value求平均后即可作为差异分析的结果。 图片 图片

2.4 为什么！为什么有了TPM/FPKM/RPKM值不能用limma，还要去做上游分析得到count值，不嫌费劲吗？ 图片

友情提示：本问题是给有一定基础的同学看的，现在不理解也可以跳过，回来看的时候记得点赞。

图片 图片

肯定有同学想用limma尝试分析RNA-seq数据了，但是我们不妨再多思考一下

图片

第一种解法：

limma帮助文档****、****DEseq2帮助文档****、****edgeR帮助文档（移动到阅读原文可观看，下不重复）

学习原始文档，你只需要把这三个帮助文档中几百页的英文单词一个一个翻译出来，把公式弄懂，你就肯定明白了！是不是特别简单

图片

第二种解法：

已知：

1. RPKM, FPKM and TPM, clearly explained（先了解基础知识）

2. 如果咱们有count矩阵，那么最优解就是edgeR和DEseq2

3. 如果是TPM/FPKM/RPKM值等，目前最优解不知道

4. count值可以转为FPKM/TPM/RPKM

求证：

表达矩阵如果为TPM/FPKM/RPKM值如何差异分析？

---

答1：如果我们有count值，且count可以转为FPKM，有没有一种可能，咱们比较一下count和FPKM差异分析结果的区别？

那么刚好我就有一个count数据（GSE140844），做出来是什么效果呢？

首先说明一下分析方法：

原始数据的count值用的是edgeR和DEseq2两种算法，

count转为FPKM后，用的是limma算法

阈值：logFC>0，*****P*****<0.05

图片

可以看出，FPKM值做limma的结果：

1. 筛选出的差异基因少

2. 和主流算法的交集少的可怜

3. 甚至还有与DEseq2/edgeR调控方向相反的结果出现

答2：如果我们有count值和FPKM值，且count和FPKM均可以转为TPM，有没有一种可能，咱们比较一下TPM和count****差异分析结果的区别

那么刚好我就有一个数据GSE162734****（count+FPKM都有），做出来是什么效果呢？

首先说明一下分析方法：

原始数据的count值用的是edgeR和DEseq2两种算法，

count/FPKM转为TPM后，用的是limma算法

阈值：logFC>0，*****P*****<0.05

可以看出，TPM值做limma的结果：

1. 筛选出的差异基因与count分析结果类似

2. 交集数量>90%

3. 交集方向基本相同

总结：

三、GEO转录组数据处理

3.1 array差异分析

3.1.1 双通道array差异分析（GSE33812）

需要用到两个文件：

两个文件：GSE的series matrix文件，GPL的annotation或full table文件，两个文件下载位置可参考四文搞定GEO数据库转录组差异分析之简介

一句话介绍我们要干什么：

如下图所示，用GPL4133中ID-gene的对应关系来注释GSE38812的ID列，并且对GSE38812中进行差异分析，得到的结果即为差异结果。

图片

准备的文件

图片

下载GPL4133的annot文件或者full table 皆可

将GPL中包含#和!的行都删除，并只保留ID列和Gene symbol列（如上图表格中GPL4133所示），保存到ann.txt中，GSE文件不用动

# 现在已默认已下载了series matrix和ann.txt文件到D盘的tem文件夹中

3.1.2 单通道array差异分析

准备文件：

图片

注：GPL文件中仅需保留ID列和Gene Symbol列，当然如果能完整理解以下代码的，也不用这样严格，之后再删也行

图片 图片

# if (!require("BiocManager", quietly = TRUE))

注意：以上例子假设基因表达谱矩阵除ID_REF列以外，前16列（连续）为实验组，后16列（连续）为对照组

如果你不想解压series matrix.gz的话，可以在R中修改probeMatrix变量的列的顺序（选择性运行）

probeMatrix <- probeMatrix[,c(18:33,2:17)] 

最终效果图：

3.1.3 High throughput sequencing

1. 所有的非count、非TPM标准化格式的数据，果断从fasq开始重新比对、定量（上游分析）

2. 本文并未介绍上游分析相关内容，有需要可以点赞、点击在看催更，点赞100+更新buxianggeng

3. TPM的数据使用单通道array差异分析中的流程即可

4. 本小节仅涉及count格式数据的DEseq2算法的分析流程

5. GEO转录组数据分析并不包含ID转换的相关教程，有需要可以点赞、点击在看催更，点赞100+更新buxianggeng

# BiocManager::install("DESeq2")

最终效果图：

注意：以上例子假设基因表达谱矩阵除genes列以外，前3列（连续）为实验组，后3列（连续）为对照组

图片

四、说在最后

1. 本期待填的坑：

图片

2. 坑都随缘填，如果有志之士想要交流的可后台留言，看见问题会回复

3. 从下期开始总结常见的生物信息分析结果图及其生物学意义

4. 感谢观看到最后