Limma求差异基因构建矩阵的两种方式

limma是我学习的第一个R包，当时只会改一下别人的流程，遇到报错就蒙圈了。 而两种矩阵的设计也是自己的心结，一直搞不懂。现在来梳理一下。

使用的数据来自于 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39461

选取前12个数据集，两种细胞系，分别敲减。3vs3和3vs3的数据 

前面的处理如下(下载的是cel文件，12个)：

1. library(affy)

2. celFiles <- list.celfiles(path = "data", full.names=TRUE)

3. data.affy <- ReadAffy(filenames = celFiles)

4. data.rma <- rma(data.affy)

5. expr.rma <- exprs(data.rma)

下面分别求abl和LNCaP细胞中敲减后的差异基因，在limma中可以一次获得

1. library(limma)

2. group <- c(rep(c("abl_con","abl_si"),3),rep(c("LNCaP_con","LNCaP_si"),3))

3. f <- factor(group,levels=c("abl_con","abl_si","LNCaP_con","LNCaP_si"))

4. # 分组矩阵

5. design <- model.matrix(~0+f)

6. colnames(design) <- levels(f)

7. fit <- lmFit(expr.rma, design)

8. # 比较矩阵

9. contrast.matrix <- makeContrasts(abl_si-abl_con,

10.                                 LNCaP_si-LNCaP_con,

11.                                 levels=design)

12. fit2 <- contrasts.fit(fit, contrast.matrix)

13. fit2 <- eBayes(fit2)

如果需要abl的差异基因,coef=1，对应makeContrasts中的ablsi-ablcon

1. topTable(fit2, coef=1,number = Inf,p.value=0.05, lfc=1.5)

如果需要LNCaP的差异基因，coef=2,对应makeContrasts中的LNCaPsi-LNCaPcon

1. topTable(fit2, coef=2,number = Inf,p.value=0.05, lfc=1.5)

最后的结果是，什么差异基因都没有，因为有批次效应，所以我们也讲了如何去除批次效应。

假如我们把abl和LNCaP细胞分开，分别求差异基因，也是可以的，这时候可以选择更加简单的方式：

以LNCaP为例,这时候情况太简单，只有处理和对照，所以连比较矩阵都可以不要的！

1. design <- c(0,1,0,1)

2. model <- model.matrix(~design)

3. fit_limma <- lmFit(data, design = model)

4. fit_limma <- eBayes(fit_limma)

5. tophits_LNCaP <- topTable(fit_limma, number = Inf, coef = 2,p.value = 0.05, lfc = 1.5)

其中coef必须等于2， 也可以是个字符串，在这里他是model矩阵设计中最后一列的名字，实际上就是design。 所以，也可以写成

1. tophits_LNCaP <- topTable(fit_limma, number = Inf, coef = "design",p.value =