前面我們處理了表達(dá)矩陣以及臨床信息，下面開(kāi)始做差異分析；
1.png
首先，我們看一下文章里是怎么做的：
2.png
這里讓我比較疑惑的是取|logFC|>1的話，并不能得到文章中580個(gè)DEG的結(jié)果，而是有上千個(gè)差異基因，為了減少DEG的數(shù)量，這里我就取|logFC|>1.5，大家有興趣的話也可以試下|logFC|>2。
3.png
一 差異分析
這里大家注意輸入輸出就OK了，在做其它數(shù)據(jù)的時(shí)候可以靈活運(yùn)用；
#導(dǎo)入上一步的數(shù)據(jù)
load('./Rdata/exp_group.Rdata')
#加載limma包做差異分析
library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(data,design)
fit=eBayes(fit)
options(digits = 4) #設(shè)置全局的數(shù)字有效位數(shù)為4
deg = topTable(fit,coef=2,adjust='BH', n=Inf)
#設(shè)定上下調(diào)基因閾值，這里根據(jù)自己的需要更改FC的值；
deg$g=ifelse(deg$P.Value>0.05,'stable',
ifelse( deg$logFC >1.5,'up',
ifelse( deg$logFC < -1.5,'down','stable') )
)
#統(tǒng)計(jì)上下調(diào)基因數(shù)量
table(deg$g)
4.png
這里我們可以發(fā)現(xiàn)取|log2FC|>1.5有1000多個(gè)DEG；
#取出差異表達(dá)基因
cg = rownames(deg)[deg$g!="stable"]
deg1=deg[cg,]
#保存至excel表中，后面需與WGCNA結(jié)果取交集
write.table(as.data.frame(deg1), file="./result/DEGs_FC1.5.xls", sep="\t", row.names=T)
二 畫(huà)圖
下面的代碼只需改輔助線那里的數(shù)值，這里取|log2FC|=1.5，取|log2FC|=2的話只需改為2和-2即可；
#畫(huà)火山圖
library(ggplot2)
p <- ggplot(
# 數(shù)據(jù)、映射、顏色
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) +
# 坐標(biāo)軸
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語(yǔ)言畫(huà)更漂亮的圖，我這有一本《R語(yǔ)言數(shù)據(jù)可視化之美：專業(yè)圖表繪制指南（增強(qiáng)版）》電子版書(shū)籍，里面有各種各樣的圖以及代碼，有需要的朋友可以后臺(tái)留言給我；
往期文章復(fù)現(xiàn)：
