没想到！生信分组还有这个大坑！你被坑过吗？！

晨曦碎碎念系列之

limma进行差异分析的策略——分组的讨论

Hi，大家好，我是晨曦，这期依旧是晨曦碎碎念哦~，今天我们聚焦的问题也是一个稍微有些难度的问题——差异分析的分组问题

这里我们主要针对的就是limma包

因为即可以进行芯片数据的分析，也可以进行测序数据的分析

所以晨曦也是经常使用这个R包，这个R包其实本质上来说，我们需要提供三类信息

表达矩阵（exprSet）

分组矩阵 (design) ：就是将表达矩阵的列（各个样本）分成几组（例如最简单的case - control，或者一些时间序列的样本day0, day1, day2 ...）【通过model.matrix()得到】

比较矩阵（contrast）：意思就是如何指定函数去进行组间比较【通过makeContrasts()得到】

好了，到这里开始犯难了，分组信息的提取又分为两种形式，以前晨曦也只会机械的调用代码，但是最近小伙伴有了这方面的困惑，所以本篇推文应小伙伴们的期望诞生

晨曦碎碎念的本质其实就是通过一篇推文解决一个问题，像是茶话会一样，所以叫“碎碎念”

我们今天就趁着这个机会把limma包分组给彻底弄明白！

晨曦碎碎念系列传送门

想白嫖单细胞生信文章？这五大源头数据库，是你发文章的源泉！高频预警！你一定要收藏！

盘活国自然的新思路！你研究的热点真的是热点吗？大数据帮你定位！

好家伙！90%以上审稿人都会问到的问题，今天帮你解决！就是这么齐齐整整！

二个组 & 三个组

针对两个组，无外乎就是对照组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)#查看数据是否经过log2 dat=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)

好啦，我们这一期的推文就给大家介绍到这里，如果各位小伙伴有感兴趣的点，也欢迎各位小伙伴在评论区留言，说不定下一个晨曦碎碎念就是你感兴趣的内容哦~

我是晨曦，我们下期再见~

晨曦单细胞笔记系列传送门

1. 首次揭秘！不做实验也能发10+SCI，CNS级别空间转录组套路全解析（附超详细代码！）

2. 过关神助！99%审稿人必问，多数据集联合分析，你注意到这点了吗？

3. 太猛了！万字长文单细胞分析全流程讲解，看完就能发文章！建议收藏！（附代码）

4. 秀儿！10+生信分析最大的难点在这里！30多种方法怎么选？今天帮你解决！

5. 图好看易上手！没有比它更适合小白入手的单细胞分析了！老实讲，这操作很sao！

6. 毕业救星！这个R包在高分文章常见，实用！好学！

7. 我就不信了，生信分析你能绕开这个问题！今天一次性帮你解决！