GEO芯片分析中的大坑,差异基因完全相反!
2019-06-11 本文已影响19人
9d760c7ce737
昨晚上吐下泄,今天躺在床上一整天,到了晚上才好一点,强撑着写一篇帖子,flag不能破啊。
GEO的芯片分析,可以看一下这一篇。
最有诚意的GEO数据库教程
GEO的芯片分析,一到了差异分析部分,几乎没有悬念,因为其中的步骤都是limma包的作者定义的,我们没有修改的空间,这里面常见的坑就是,
差异分析的结果会完全相反,原本高表达的会变成低表达。
不懂内在逻辑的时候,想要正确完全靠的是运气,主要看两组的命名,比如,一个3vs3的实验,我们的分组是这个样子的。
group <- c(rep("con",3),rep("treat",3))
前面三个是control,后面三个是treat。这时候,算出来的结果就是treat比上control,没有问题。因为默认情况下因子的排序是按照字母的前后顺序排列的,control在treat之前。
如果我们把两组命名为,cancer 和normal
group <- c(rep("cancer",3),rep("normal",3))
那么由于n在c的后面,所以最终的结果会变成normal比上cancer,所以的差异基因变化值都会相反。
如何 避免呢?这里要用到因子的排序。只要每次在levels里面,把真实的对照组放在前面就可以了。
group <- c(rep("con",3),rep("treat",3))
group <- factor(group,levels = c("con","treat"),ordered = F)
换一个名字也是一样的,只要把对照组放在前面即可
group <- c(rep("treat",8),rep("vector",10))
group <- factor(group,levels = c("vector","treat"),ordered = F)
设置了分组后,接下来的操作是固定的
## 1.构建比较矩阵
design <- model.matrix(~group)
## 比较矩阵命名
colnames(design) <- levels(group)
design
#2.线性模型拟合
fit <- lmFit(exprSet,design)
#3.贝叶斯检验
fit2 <- eBayes(fit)
最终只有输出差异基因的时候要强调一下,coef要选择2。
allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf)
其中,group不一定是要很整齐的前面是对照,后面是处理,样本的顺序可以是乱的。比如这一篇中的group也是可以的。
GEO的样本名称太多而且排序不规则,你们都是手动分组的么?
明天我们介绍一下limma包中的多组差异分析,三组并不需要算三次,一次就可以了。