前面我们处理了表达矩阵以及临床信息，下面开始做差异分析；
1.png
首先，我们看一下文章里是怎么做的：
2.png
这里让我比较疑惑的是取|logFC|>1的话，并不能得到文章中580个DEG的结果，而是有上千个差异基因，为了减少DEG的数量，这里我就取|logFC|>1.5，大家有兴趣的话也可以试下|logFC|>2。
3.png
一 差异分析
这里大家注意输入输出就OK了，在做其它数据的时候可以灵活运用；
#导入上一步的数据
load('./Rdata/exp_group.Rdata')
#加载limma包做差异分析
library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(data,design)
fit=eBayes(fit)
options(digits = 4) #设置全局的数字有效位数为4
deg = topTable(fit,coef=2,adjust='BH', n=Inf)
#设定上下调基因阈值，这里根据自己的需要更改FC的值；
deg$g=ifelse(deg$P.Value>0.05,'stable',
ifelse( deg$logFC >1.5,'up',
ifelse( deg$logFC < -1.5,'down','stable') )
)
#统计上下调基因数量
table(deg$g)
4.png
这里我们可以发现取|log2FC|>1.5有1000多个DEG；
#取出差异表达基因
cg = rownames(deg)[deg$g!="stable"]
deg1=deg[cg,]
#保存至excel表中，后面需与WGCNA结果取交集
write.table(as.data.frame(deg1), file="./result/DEGs_FC1.5.xls", sep="\t", row.names=T)
二 画图
下面的代码只需改辅助线那里的数值，这里取|log2FC|=1.5，取|log2FC|=2的话只需改为2和-2即可；
#画火山图
library(ggplot2)
p <- ggplot(
# 数据、映射、颜色
deg, aes(x = logFC, y = -log10(P.Value), colour=g)) +
geom_point(alpha=0.4, size=3.5) +
scale_color_manual(values=c("#c8e09f","#bdbdbf", "#ec5141"))+
# 辅助线
geom_vline(xintercept=c(-1.5,1.5),lty=4,col="black",lwd=0.8) +
geom_hline(yintercept = -log10(0.05),lty=4,col="black",lwd=0.8) +
# 坐标轴
labs(x="log2(fold change)",
y="-log10 (p-value)")+
ggtitle("GSEGSE140797")+
theme_bw()+
# 图例
theme(plot.title = element_text(hjust = 0.5),
legend.position="right",
legend.title = element_blank())
p
5.png
大家如果想用R语言画更漂亮的图，我这有一本《R语言数据可视化之美：专业图表绘制指南（增强版）》电子版书籍，里面有各种各样的图以及代码，有需要的朋友可以后台留言给我；
往期文章复现：
网友评论