这是一篇发表在《Circulation》的文章，图乍一看确实是一个雷达图，而且我们之前也讲过一种雷达图的做法（复现Immunity文章图表：分组雷达图展示富集结果）、（雷达图展示GSEA分析结果）。雷达图有很多种做法。这篇文章中的图我们先看一下它表达的意思：首先是图展示的是3各组的，可以理解为基因或者蛋白以及产物的表达高低，然后圈外有点，表示的是三个组两两比较的显著性！这个图还是挺有意思的，所以我们复现一下。虽然花里胡哨！
reference：Myocardial Metabolomics of Human Heart Failure With Preserved Ejection Fraction
我一开始的想法是使用雷达图相关的方式进行，但是后来发现不太好实现，也就没有进行其他的尝试了，不想太浪费时间，我觉得应该会有方法。后来我一下子想到的是利用极坐标进行绘制，那么之前图里面的线条其实就是折线图而已，显著性我们使用点来添加就可以了，ggplot实现没有压力，这么一想，那么事情就很好办了。首先我用单细胞的数据构建了一个作图数据。我们的数据里面刚好也有三个分组！至于其他的数据，自行整理成为相似的格式即可：首先我们做三组比较的差异基因！
library(Seurat)
library(dplyr)
library(ggplot2)
SMC <- subset(uterus, celltype=='Smooth muscle cells')
rm(uterus)
table(SMC$orig.ident)
# AEH EEC HC
# 959 605 575
#找找差异基因构建数据
AEH_vs_HC <- FindMarkers(SMC,
ident.1 = "AEH",
ident.2 = "HC",
group.by = "orig.ident",
logfc.threshold = 0.25,
min.pct = 0.25,
only.pos = T)
AEH_vs_HC <- AEH_vs_HC[which(AEH_vs_HC$p_val_adj<0.05),]
EEC_vs_HC <- FindMarkers(SMC,
ident.1 = "EEC",
ident.2 = "HC",
group.by = "orig.ident",
logfc.threshold = 0.25,
min.pct = 0.25,
only.pos = T)
EEC_vs_HC <- EEC_vs_HC[which(EEC_vs_HC$p_val_adj<0.05),]
EEC_vs_AEH <- FindMarkers(SMC,
ident.1 = "EEC",
ident.2 = "AEH",
group.by = "orig.ident",
logfc.threshold = 0.25,
min.pct = 0.25,
only.pos = T)
EEC_vs_AEH <- EEC_vs_AEH[which(EEC_vs_AEH$p_val_adj<0.05),]
然后选择需要展示的基因，这里其实是构建三组显著性的文件：
#选择的基因
genes <- c("HLA-A","HLA-C","HLA-B","IFITM3","ADARB2",
"TIMP1","SCGB1D2","CYTOR","CFH","IGHM","C2","CD63","CFI")
AEH_vs_HC_sig <- AEH_vs_HC[genes,] %>%
mutate(name = genes, sig=ifelse(is.na(p_val),"no_sig","sig"),group="AEH_vs_HC")%>%
select(name, sig, group)
EEC_vs_HC_sig <- EEC_vs_HC[genes,]%>%
mutate(name = genes, sig=ifelse(is.na(p_val),"no_sig","sig"),group="EEC_vs_HC")%>%
select(name, sig, group)
EEC_vs_AEH_sig <- EEC_vs_AEH[genes,]%>%
mutate(name = genes, sig=ifelse(is.na(p_val),"no_sig","sig"),group="EEC_vs_AEH")%>%
select(name, sig, group)
接下来获取基因在三个组的表达量：
#表达量数据
DefaultAssay(SMC) <- 'RNA'
avg_exp <- AverageExpression(SMC,group.by = 'orig.ident',features = genes)[[1]]
avg_exp <- avg_exp[genes,]
library(reshape2)
data_exp <- melt(avg_exp,
measure.vars = 1:3,#选择需要转化的列
value.name = 'exp')#聚合值的新列名
sig_gene_group <- rbind(AEH_vs_HC_sig,EEC_vs_HC_sig,EEC_vs_AEH_sig)
data_exp <- cbind(data_exp,sig_gene_group)
colnames(data_exp)[6] <- 'sig_group'
最后作图：
ggplot(data_exp, aes(x=gene,y=exp,group=group,colour=group)) +
geom_line(lwd=1.5) +
geom_point(data=subset(AEH_vs_HC_sig, sig=='sig'),aes(x=name, y=max(data_exp$exp)+1),size=3)+
geom_point(data=subset(EEC_vs_HC_sig, sig=='sig'),aes(x=name, y=max(data_exp$exp)+5),size=3)+
geom_point(data=subset(EEC_vs_AEH_sig, sig=='sig'),aes(x=name, y=max(data_exp$exp)+9),size=3)+
coord_polar(start=0,direction=1, clip = "on")+
xlab("") +
ylab("")+
theme(axis.ticks = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(color = 'black',size = 10),
panel.background = element_blank(),
panel.grid.major = element_line(colour = "grey80",size = 0.8),
panel.grid.minor = element_line(colour = "grey80",size = 0.8))+
scale_color_manual(values = c( "#B84D64","#3D6AAA","#EE7072","tan1","#394D9B","red"))
我们可以发现，这个图基本上全部复现了，只需要自己调整一下legend就可以了。此外还有一个问题就是我的线条上没有置信区间，这是因为这里我用的单细胞数据的平均表达量，所以只有一个数据，故而没有。如果你的数据有重复，那就就可以添加置信区间，效果如下：
网友评论