主讲老师第二十五课：差异基因分析
PCA 分析以后，我们开始进入差异分析（DEGs.zip），选择差异表达基因。
首先，还是需要分析数据的导入与预处理，导入 .RData 文件，同时把一些低表达基因给过滤了。
注意，这一步筛选作用十分重要，必须进行，不然后面会有报错。在正式的介绍之前，先给大家介绍一篇关于 RNA-seq 分析的文献，完整阐述了 RNA 测序的整体分析流程，题目为 RNA-Seq workflow:
gene-level exploratory analysis and differential expression
1.deseq2
1.1 R 包的安装与读取
来看下如何使用 DESeq2 进行差异分析。
1.2 差异分析
首先，我们来设置一下分组信息。
同时，将分组设置为因子形式，注意转换，不然到时候很有可能 high 和 low 相反，引起所有结果都相反，到此为止，差异分析的准备工作基本就完成了，后面，就是基于 deseq2包的常规分析流程将 count 文件和分组信息进行设置后，R 包内部会进行一些的标准化操作，最终分析得到差异信息，由于这部分时间比较长，我提前把这部分的运行结果进行了保存，随后，我们来把差异分析的结果进行相应的提取关于每一项内容的具体含义，我们也可以来查看一下。
当然，也可以整体看一下分析的结果
接着，将其按校正后的 p 值进行排列，并转换成为数据框格式，并进行保存。
随后，根据 ｜logFC｜ >2 和 padj <0.05 的筛选标准，寻找显著差异表达基因。
最终，得到显著上调基因 1885 个，显著下调基因 1404 个，最后，将结果输出保存。
2.edgeR
2.1 R 包安装与读取
接着，我们来看下如何使用 edgeR 来进行测序数据的差异分析
2.2 差异分析
首先，与前面类似，对分组信息进行相应的设置，获得 DGEList接着，使用 edgeR 的 calcNormFactors 函数默认 TMM 算法对 DGEList 标准化然后，对实验设计矩阵(Design matrix)进行设置，这一步类似于 DESeq2 中的 design 参数。
1.估计离散值（Dispersion）
2.通过 glmFit() 进行拟合
当然也可以用 glmQLFit（quasi-likelihood (QL)拟合）用于解释生物学和技术性导致的基因特异性变异。
3.差异表达检验，因为前面使用的是 glmFit，那么对应的就是 glmLRT
4.提取结果
这样，整个分析部分就结束了，同样的，用相同的筛选标准，我们来挑选显著差异表达基因。
最终，可以看到，现在上调基因为 2776 个，显著下调基因为 1339 个，同时，保存分析结果。
3.limma
关于 limma 包的差异分析，其最初是为芯片的差异分析而设计的，后来对 RNA 测序数据也同样可以进行差异分析，一般是自己读取 RNA-seq 的基因的 reads 的 counts 数进行分析。
3.1 R 包安装与读取
3.2 差异分析
首先，通过对组别进行设计，然后根据精确权重法（precision weights）VOOM 算法，对输入数据进行标准化处理。
接下来的使用，和 limma 包基本就比较类似的了最终，通过 topTable() 函数提取所有基因的差异分析结果，接下来就是对差异分析进行注释，选出显著上调和显著下调对基因。
保存差异分析的结果
到此，整个差异分析基本就完成了，虽然整个过程看似繁琐，但是，实际上自己准备好输入数据后，设计分组信息，其他的所有代码都不需要进行修改，即可完成差异分析。另外，关于三个 R 包分析得到的结果不一样，这个是十分正常的，虽然三者要求的输入数据都是count 数据，但是其内部的标准化流程是不一样的，因此，得到的结果有所差异在所难免
limma进行差异分析的策略——分组的讨论
其实针对limma包，因为即可以进行芯片数据的分析，也可以进行测序数据的分析，所以主讲老师也是经常使用这个R包，这个R包其实本质上来说，我们需要提供三类信息表达矩阵 （exprSet）分组矩阵 (design) ：就是将表达矩阵的列（各个样本）分成几组（例如最简单的 case - control ，或者一些时间序列的样本 day0, day1, day2 ... ）【通过 model.matrix() 得到】比较矩阵（contrast）：意思就是如何指定函数去进行组间比较【通过 makeContrasts() 得到】好了，到这里开始犯难了，分组信息的提取又分为两种形式，以前主讲老师也只会机械的调用代码，但是最近小伙伴有了这方面的困惑，所以本篇推文应小伙伴们的期望诞生主讲老师碎碎念的本质其实就是通过一篇推文解决一个问题，像是茶话会一样，所以叫“碎碎念”好了，我们今天就趁着这个机会把 limma 包分组给彻底弄明白二个组&三个组针对两个组，无外乎就是对照组VS实验组，癌症VS健康等等，我们只需要通过提取临床信息中的关键词即可获得分组信息，其实两个组和三个组甚至于多个组在提取分组信息上来说是一样的，重点是后面的差异分组，或者可以说是构建分组信息有着不同的选择#准备工作library(AnnoProbe)library(GEOquery)#获取数据并处理表达矩阵gset <- geoChina("GSE12021")gset[[1]]#获得的Expressionset提取第一个组成a=gset[[1]]#通过赋值避免更改原始数据dat=exprs(a)#提取表达矩阵dim(dat)#看一下dat这个矩阵的维度boxplot(dat[,1:4],las=2)#查看数据是否经过log2dat=log2(dat)boxplot(dat[,1:4],las=2)#数据分布并不是十分均匀，所以我们需要标准化一下library(limma)dat=normalizeBetweenArrays(dat)boxplot(dat[,1:4],las=2)#再次查看数据，得到标准的原始数据#因为我们limma分析芯片数据需要满足两个条件：log2和counts#获取临床信息pd=pData(a)这里我们主要来看一下临床信息，这里我们的分组有三个，分别是正常组、OA（骨关节炎）、RA（类风湿关节炎）Ps:这个数据集各个样本的 title 中有一个单词拼写错误，所以后面提取相关分组信息的时候多加了一个循环这里的临床信息又有小小的划分，如果有直接的分组信息（例如：Normal or Cancer），我们可以用"=="来精确匹配，进而提取，但是如果遇到这种有很长信息的，我们可以使用下面的代码来进行一个分组信息的获取#获取分组信息library(stringr)library(tidyverse)table(pd$title)group_list <- as.data.frame(str_split_fixed(pd$title,",",n=2))[1]group_list <- c(group_list)[[1]]group_list <- ifelse(group_list=="Normal","Normal",ifelse(group_list=="Osteoarthrisits","OA",ifelse(group_list=="Osteoarthritis","OA","RA")))table(group_list)#Ps，这里OA的对应英文有两类，应该是上传数据的作者整理错误
然后获取分组信息后，我们就可以进行标准流程的差异分析了
#构建分组矩阵library(limma)design <- model.matrix(~0+group_list)colnames(design) <- levels(factor(group_list))rownames(design)=colnames(dat)#构建比较矩阵contrast.matrix<-makeContrasts("OA-Normal",levels = design)
#limma差异分析标准三步骤
fit <- lmFit(exprSet, design)#线性拟合fit2 <- contrasts.fit(fit, contrast.matrix) #这种分组模式下特有的，基本不需要改动fit2 <- eBayes(fit2)#贝叶斯检验allDiff=topTable(fit2,adjust='fdr',coef=1,number=Inf) #提取差异分析结果这里我们重点先来看一下我们这个设计矩阵这个矩阵很好理解，就是告诉我们样本属于哪一个分组但是在这里可能会有人犯难 model.matrix 函数这里有两种编写方式，简单来说就是加0和不加0，这里我以前是拼命想理解，什么截距列，不充当截距等等，直到有一天某一位老师的点拨才恍然大悟，这个加0和不加0统称为参数，这是作者定义的，你如果要使用这个包你就按着我的规矩来，否则，你就自己去写包，好了，瞬间敞亮，豁然开朗，土地平旷，屋舍俨然......当然我们也可以再理解一下不加0已经把第一列做成截距列，这里我们把哪一列做成截距列是由因子水平控制的，因子水平第一个作为截距列，其构建的矩阵天生存在了两种关系即明确了分组的同时又明确了对比关系，就是从第二列开始往后用coef这个参数可以提取固定列数与截距列的比，所以这样构建矩阵是不需要另外构建比较矩阵的加0那种方法，仅仅是分组而已，组之间需要如何比较，需要自己再制作差异比较矩阵，通过makeContrasts函数来控制如何比较，这时候我们的coef参数就是通过选择makeContrasts函数
里的分组来获得差异分析的结果
我个人是比较推荐加0的构建方法的，因为可以优雅的完成多组的差异分析结果，当然，如果你选择不加
0，对于两组的差异分析结果是一样的
接下来我们再来通过两组和三组的分组设计来加深一下我们对于 limma 包构建分组矩阵的理解吧：
#示例（两组）group <- c(rep("con",3),rep("treat",3)) #分两组group <- factor(group,levels = c("con","treat"),ordered = F)#转化因子且规定顺序，因为limma包的设计者默认因子顺序的第一位作为截距列，所以这里我们需要把control放在前面design <- model.matrix(~group)colnames(design) <- levels(group)design#此时构建的矩阵第一列充当截距列的就是control，所以我们后续通过coef参数指定列数的时候，都是指定列数比截距列#标准三步骤fit <- lmFit(exprSet,design)fit2 <- eBayes(fit)allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf)当然主讲老师这里更推荐的就是第二种方法，示例代码如下#示例（三组）group <- c(rep("con",2),rep("drugA",2),rep("drugB",2)) #分三组## 分组矩阵design <- model.matrix(~0 + group)## 分组矩阵命名colnames(design) <- levels(group)design## 比较矩阵，其中treat - con根据自己的组别更改（这里是实验组-对照组，这个位置限制了比对的顺序，和因子水平无关）contrast.matrix <- makeContrasts(drugA - con,drugB - con,drugA - drugB,levels=design)## 线性拟合fit <- lmFit(exprSet, design)fit2 <- contrasts.fit(fit, contrast.matrix) #这种分组模式下特有的，基本不需要改动## 贝叶斯检验fit2 <- eBayes(fit2)## 提取差异结果，注意这里的coef是1allDiff=topTable(fit2,adjust='fdr',coef=1,number=Inf)#最后如何获取结果呢，只要改变coef的值即可，1，2，3分别对应drugA-con, drugB-con,drugA-drugB#adjust="fdr"表示对p值的校正方法，包括了："none", "BH", "BY" and "holm"。
limma_DEG = na.omit(output) #如果数据中有NA，也可以直接去除当然这里我们最优秀的教程永远是 limma 的帮助文档，大家如果有疑问也可以去看帮助文档来进行学习usersguide.pdf (bioconductor.org)好啦，我们这一期的推文就给大家介绍到这里，如果各位小伙伴有感兴趣的点，也欢迎各位小伙伴在评论区留言，说不定下一个主讲老师碎碎念就是你感兴趣的内容哦~我是主讲老师，我们下期再见~
参考教程：
那些常用的limma操作 - 简书 (jianshu.com)表达量矩阵分组很复杂也可以使用limma的3大策略 - 云+社区 - 腾讯云 (tencent.com)limma 差异分析透彻讲解 - 简书 (jianshu.com)(11条消息) 用limma包进行多组差异表达分析_今天也是个妖精头子呀的博客-CSDN博客R语言limma包差异基因分析（两组或两组以上） - 知乎 (zhihu.com)用limma对芯片数据做差异分析 (bio-info-trainee.com)用limma包的voom方法来做RNA-seq 差异分析 (bio-info-trainee.com)