之前有小伙伴咨询一幅差异基因展示的图,我随口就说这类似于火山图,就是点图,ggplot就可以完成,不晓得他最后做出来没有[图片上传失败...(image-b8ad3b-1670206965508)]
。对了,之前Bulk RNA转录组系列关于DESeq2的帖子,封面图就是这样的。术语叫做MAplot,展示数据分布情况的可视化!最近在文献中看到这样的差异基因展示方式,所以复现一下。这里我们从counts矩阵的转录组差异分析开始做这个图,很简单,原理就是散点图,按照火山图的做法筛选差异基因即可。重试差异分析:
setwd('D:/KS项目/公众号文章/MAplot')
data <- read.csv('counts.csv', header = T, row.names = 1)
genecounts <- data[apply(data, 1, sum) > 0, ]
SampleAnno <- data.frame(condition = factor(rep(c('neg', 'pos'), each = 3),
levels = c('neg', 'pos')))
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(DESeq2)
library(ggrepel)
dds <- DESeqDataSetFromMatrix(countData=genecounts, colData=SampleAnno, design= ~condition)
dds1 <- dds
dds1 <- DESeq(dds1)
res <- results(dds1, contrast = c('condition', 'pos', 'neg'))
res <- as.data.frame(res)
write.csv(res, file = 'DEGs.csv')#差异分析结果
构建数据:
df <- as.data.frame(counts(dds1, normalized=TRUE))
df_merge = merge(res,df,by ='row.names')
rownames(df_merge) = df_merge$Row.names
df_merge = select(df_merge,-1)
#挑选需要展示的差异基因
genes = c("CLSTN2",
"MME",
"SYNPO",
"ENC1",
"PDGFRB",
"TAGLN3",
"ZFHX4",
"THY1",
"CYP26A1",
"RGS4",
"CDH10",
"JAM2",
"LMX1B",
"NTRK1",
"ST8SIA1",
"FUT9",
"DOK5",
"CABP7",
"APOL6",
"CHI3L1")
添加label:
df_merge$label = ifelse(rownames(df_merge) %in% genes,rownames(df_merge),"")
df_merge$condition = ifelse(df_merge$log2FoldChange >= 1 & df_merge$padj < 0.01,"UP",
ifelse(df_merge$log2FoldChange <= -1 & df_merge$padj < 0.01,"DOWN","none"))
df_merge$condition = ifelse(rownames(df_merge) %in% genes,"label",df_merge$condition)
ggplot2作图:
ggplot(df_merge,aes(x=baseMean,y=log2FoldChange)) +
geom_point(aes(color=condition,size = condition,alpha = condition),
show.legend = F) +
geom_hline(yintercept = c(-1,1),lty=2,lwd = 1) +
theme_bw(base_size = 12) +
ggtitle("Neg vs Pos") +
scale_color_manual(values = c("cadetblue","black","grey","darkorange")) +
scale_size_manual(values = c(1.5,2,1.5,1.5)) +
scale_alpha_manual(values = c(0.6,1,0.6,0.6)) +
scale_x_log10() +
theme(panel.grid=element_blank()) +
scale_y_continuous(limits = c(-10,8),breaks = c(-10,-5,0,4,8)) +
xlab("Mean Normalized Counts") +
ylab("log2FoldChange")+
geom_text_repel(aes(label=label), color="black",fontface="italic",
size=4, segment.size=0.5,hjust=0,
nudge_x=1, nudge_y = 1)
效果和火山图差不多,就是一个展现形式而已,萝卜青菜各有所爱!示例数据及代码已上传群文件!更多精彩请至我的公众号---KS科研分享与服务!
网友评论