GEO芯片如果超过了两组，也可以一次搞定差异分析

昨天我们讲到，GEO芯片只有两组，如何做差异分析，以及如何避免差异分析的结果完全相反的情况。
GEO芯片分析中的大坑，差异基因完全相反！

留言中有朋友提到，可以用makeContrasts函数限定，对于这种情况，我只能用赵飞的话来回复一下：裁判请回到裁判席！

实际上，这位朋友说的是limma差异矩阵构建的第二种方法，在很多时候被认为是首选方法，只是有点复杂而已。
我们先看如果用昨天的方法，一套差异分析是什么样子的

group <- c(rep("con",3),rep("treat",3)) 
group <- factor(group,levels = c("con","treat"),ordered = F)
design <- model.matrix(~group)
colnames(design) <- levels(group)
design
fit <- lmFit(exprSet,design)
fit2 <- eBayes(fit)
allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf) 

如果用今天的方法是什么样子的呢？

## 分组
group <- c(rep("con",3),rep("treat",3)) 
## 因子化
group <- factor(group)
## 分组矩阵
design <- model.matrix(~0 + group)
## 分组矩阵命名
colnames(design) <- levels(group)
design
## 比较矩阵，其中treat - con根据自己的组别更改
contrast.matrix <- makeContrasts(treat - con, levels=design)
contrast.matrix
## 线性拟合
fit <- lmFit(exprSet, design)
fit2 <- contrasts.fit(fit, contrast.matrix) 
## 贝叶斯检验
fit2 <- eBayes(fit2)
## 提取差异结果，注意这里的coef是1
allDiff=topTable(fit2,adjust='fdr',coef=1,number=Inf) 

多了几个函数，我讲一下区别，
第一个是分组矩阵

## 分组矩阵
design <- model.matrix(~0 + group)

这里多了个0,一开始我是拼了命的要理解，加0和不加0有什么区别，原理是什么，最终健明告诉我，这个叫做参数，是作者自定义的，用来区分两种不同的情况，如果你使用作用的包，就按照他的要求来就行了。
第二个是多了个比较矩阵

## 比较矩阵，其中treat - con根据自己的组别更改
contrast.matrix <- makeContrasts(treat - con, levels=design)

这个表达式中，我们能修改的只有treat - con, 根据组别不同，可以变成cancer - normal, 正是有了这个函数在限定比较的顺序，所以这种方法是不会出错的。
第三个是多了一个contrasts.fit函数

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

但是，不需要我们作任何修改。
总结一下，

整个上面的代码，需要我们修改的仅仅是group分组信息，以及makeContrasts中的比较信息，其他的都不需要修改。

其中，唯一不知道的是exprSet，这实际上是表达矩阵，格式是这个样子的，行名为基因名称，列名是样本名称。

如果这种方法就这么大神通，我肯定会选第一个方案，因为明显简单一点。但是，第二种方法，可不得了，他可以同时计算多组的差异分析,比如以下有三组。

group <- c(rep("con",2),rep("drugA",2),rep("drugB",2)) 

如果想比较drugA-con, drugB-con,drugA-drugB, 我们只要用上面的方法计算三次就可以了。但是有点麻烦，不够优雅，第二种方法只要makeContrasts中增加需要比较的组别，就可以一次搞定。

group <- c(rep("con",2),rep("drugA",2),rep("drugB",2)) 
group <- factor(group)
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
design
contrast.matrix <- makeContrasts(drugA - con, 
                                 drugB - con, 
                                 drugA - drugB, 
                                 levels=design)
contrast.matrix
fit <- lmFit(exprSet, design)
fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2)

重点就是看makeContrasts函数，里面可以自定义任何组别之间的比较。
最后如何获取结果呢，只要改变coef的值即可，1，2，3分别对应drugA-con, drugB-con,drugA-drugB。

allDiff1=topTable(fit2,adjust='fdr',coef=1,number=Inf) 
allDiff2=topTable(fit2,adjust='fdr',coef=2,number=Inf) 
allDiff3=topTable(fit2,adjust='fdr',coef=3,number=Inf)

这个1，2，3实际上指的是，contrast.matrix中的三列

所以coef可以用列名来代替，获得的结果一样

allDiff1=topTable(fit2,adjust='fdr',coef="drugA - con",number=Inf) 
allDiff2=topTable(fit2,adjust='fdr',coef="drugB - con",number=Inf) 
allDiff3=topTable(fit2,adjust='fdr',coef="drugA - drugB",number=Inf)

好了，limma差异分析两种方法都讲了，我们明天讲一下GEO芯片中配对样本的差异分析。
我是果子，明天见。这将是一个长达6个月的连续写作。