之前有小伙伴咨詢一幅差異基因展示的圖，我隨口就說這類似于火山圖，就是點(diǎn)圖屑柔，ggplot就可以完成屡萤，不曉得他最后做出來沒有[圖片上傳失敗...(image-b8ad3b-1670206965508)]
。對(duì)了掸宛，之前Bulk RNA轉(zhuǎn)錄組系列關(guān)于DESeq2的帖子死陆，封面圖就是這樣的。術(shù)語叫做MAplot唧瘾，展示數(shù)據(jù)分布情況的可視化翔曲！圖片
最近在文獻(xiàn)中看到這樣的差異基因展示方式，所以復(fù)現(xiàn)一下劈愚。這里我們從counts矩陣的轉(zhuǎn)錄組差異分析開始做這個(gè)圖，很簡(jiǎn)單闻妓，原理就是散點(diǎn)圖菌羽，按照火山圖的做法篩選差異基因即可。重試差異分析：
setwd('D:/KS項(xiàng)目/公眾號(hào)文章/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')#差異分析結(jié)果
構(gòu)建數(shù)據(jù)：
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)
圖片
效果和火山圖差不多由缆，就是一個(gè)展現(xiàn)形式而已注祖，蘿卜青菜各有所愛！示例數(shù)據(jù)及代碼已上傳群文件均唉！更多精彩請(qǐng)至我的公眾號(hào)---KS科研分享與服務(wù)是晨！
