github-actions[bot]
commented
3 years ago GEO芯片中配对样本如何做差异分析 by 果子学生信
前面我们已经见了GEO芯片差异分析两组如何做差异分析
GEO芯片分析中的大坑，差异基因完全相反！
接着讲了，差异分析的另外一种方式，以及多组如何同时做差异分析。
GEO芯片如果超过了两组，也可以一次搞定差异分析
今天我们讲一下，配对样如何做差异分析。该部分内容取材于这个诚意GEO帖子。
最有诚意的GEO数据库教程
如果不配对，分析应该是这样的
library(limma)
## 如果不是配对的
group <- c(rep('before',18),rep('after',18))
group <- factor(group,levels = c("before","after"),ordered = F)
design <- model.matrix(~group)
fit <- lmFit(exprSet,design)
fit2 <- eBayes(fit)
最终可以通过添加p.value直接提取差异基因
allDiff1=topTable(fit2,adjust='fdr',coef=2,number=Inf ,p.value=0.05)
有6796个差异基因，部分数据如下
如果要是配对样本，这个例子中是1:18和19:36依次配对。差异分析这样做
第一种形式
pairinfo = factor(rep(1:18,2))
group <- c(rep('before',18),rep('after',18))
group <- factor(group,levels = c("before","after"),ordered = F)
design <- model.matrix(~group+pairinfo)
fit <- lmFit(exprSet,design)
fit2 <- eBayes(fit)
通过添加配对信息pairinfo到model.matrix函数中实现。
最终获取差异基因有两种方式：
allDiff2=topTable(fit2,adjust='fdr',coef=2,number=Inf ,p.value=0.05)
allDiff21=topTable(fit2,adjust='fdr',coef="groupafter",number=Inf ,p.value=0.05)
差异基因是7635，多出来的基因就是因为加了配对信息才显示出差异的，我们后面作图来展示。
第二种形式
如果采取昨天提到的第二种方法，可以这样做
pairinfo = factor(rep(1:18,2))
group <- c(rep('before',18),rep('after',18))
group <- factor(group)
design <- model.matrix(~0 + group +pairinfo)
colnames(design)[1:length(levels(group))] <- levels(group)
contrast.matrix <- makeContrasts(after - before, levels=design)
## 线性拟合
fit <- lmFit(exprSet, design)
fit <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit)
## 提取差异结果，注意这里的coef是1
allDiff3=topTable(fit2,adjust='fdr',coef=1,number=Inf ,p.value=0.05)
最终获得的结果是一样的。
那么现在的问题是，这些配对的数据如何展示，最好的办法是每个配对的样本间给出连线。
我们可以通过ggplot2来实现。
首先，创建```ggplot2```喜欢的数据格式：
pairinfo = factor(rep(1:18,2))
group <- c(rep('before',18),rep('after',18))
group <- factor(group,levels = c("before","after"),ordered = F)
data <- data.frame(pairinfo=pairinfo,group=group,t(exprSet))
使用ggplot2来作图
高表达的
library(ggplot2)
## 高表达
ggplot(data, aes(group,CYP1B1,fill=group)) +
geom_boxplot() +
geom_point(size=2, alpha=0.5) +
geom_line(aes(group=pairinfo), colour="black", linetype="11") +
xlab("") +
ylab(paste("Expression of ","CYP1B1"))+
theme_classic()+
theme(legend.position = "none")
低表达的
library(ggplot2)
ggplot(data, aes(group,CH25H,fill=group)) +
geom_boxplot() +
geom_point(size=2, alpha=0.5) +
geom_line(aes(group=pairinfo), colour="black", linetype="11") +
xlab("") +
ylab(paste("Expression of ","CH25H"))+
theme_classic()+
theme(legend.position = "none")
这些基因配不配对，都有差异，所以我们找出那个相差的1000多个墙头草基因来看一下。
genelist <- setdiff(rownames(allDiff3),rownames(allDiff1))
选前4个基因来分别普通作图和配对作图比较一下
## 第一个方式
plotlist <- lapply(genelist[1:4], function(x){
data =data.frame(data[,c("pairinfo","group")],gene=data[,x])
ggplot(data, aes(group,gene,col=group)) +
geom_jitter(size=2, alpha=0.5) +
xlab("") +
ylab(paste("Expression of ",x))+
theme_classic()+
theme(legend.position = "none")
})
library(cowplot)
p1 <- plot_grid(plotlist=plotlist, ncol=4,labels = LETTERS[1:4])
## 第二种方式
plotlist <- lapply(genelist[1:4], function(x){
data =data.frame(data[,c("pairinfo","group")],gene=data[,x])
ggplot(data, aes(group,gene,fill=group)) +
geom_boxplot() +
geom_point(size=2, alpha=0.5) +
geom_line(aes(group=pairinfo), colour="black", linetype="11") +
xlab("") +
ylab(paste("Expression of ",x))+
theme_classic()+
theme(legend.position = "none")
})
library(cowplot)
p2 <- plot_grid(plotlist=plotlist, ncol=4,labels = LETTERS[1:4])
plot_grid(p1,p2,ncol=1)
很明显，差异分析的配对作图效果要好很多啊。
好了，明天我们要扩展一下这个贴子中，对因子的level进行排序的那个技能，到底还有什么实用的功能。
GEO芯片分析中的大坑，差异基因完全相反！
这是果子开启的长达6个月连续写作计划的第4天，
我们明天见。
https://mp.weixin.qq.com/s/tz0CsJDumvnzY8WLqvoT-A