明明是一个热图就能搞定的事情为什么要复杂到蛋壳图呢
前些天的教程：直接为CellPhoneDB创建一个独立的conda环境，以及：把Seurat对象里面表达量矩阵和细胞表型信息输出给CellPhoneDB做细胞通讯，给大家演示了如何对pbmc3k单细胞数据集做细胞通讯，并且在：CellPhoneDB的单细胞通讯结果的理解 给大家演示了细胞通讯结果的多个txt文件的含义。
并且做了一个简单的可视化，见：CellPhoneDB的单细胞通讯结果的可视化之气泡图，差不多让大家理解了所谓的细胞通讯， 就是在两个不同单细胞亚群里面，各自高表达受体配体基因对里面的一个。比如CD74_MIF 和 C5AR1_RPS19 这两个配对情况，在mono到cd4的方向比较明显，所以它们就是CD74和C5AR1在mono高表达，而MIF和RPS19在cd4的t细胞里面高表达。
然后大家最感兴趣的是每个单细胞数据集里面的多个单细胞亚群各自两两之间的受体配体基因对的数量，其实就是一个热图，以及对应的数据，值得注意的是这个 单细胞亚群各自两两之间的受体配体基因对的数量 统计文件，其实是需要使用CellPhoneDB自带的一些统计绘图函数：
conda activate cellphonedb
# 必须要保证当前路径下面有前面的步骤输出的out文件夹哦
cellphonedb plot dot_plot
cellphonedb plot heatmap_plot cellphonedb_meta.txt
# 做完后为了跟前面的区分，我把 out文件夹，修改名字为 out-by-symbol 文件夹啦
主要是为了得到 count_network.txt 文件，里面的内容就是 单细胞亚群各自两两之间的受体配体基因对的数量 统计结果。但是不少粉丝留言表示他自己的 cellphonedb plot dot_plot 代码会失败，因为服务器里面的R语言环境问题。
其实这个 count_network.txt 文件，是可以自己写代码的。如下所示：
rm(list = ls())
getwd()
mypvals <- read.table("./out-by-symbol/pvalues.txt",
header = T,sep = "\t",stringsAsFactors = F)
mymeans <- read.table("./out-by-symbol/means.txt",
header = T,sep = "\t",stringsAsFactors = F)
colnames(mypvals)[12:ncol(mypvals)]
sm = as.data.frame(
do.call(rbind,
lapply( 12:ncol(mypvals) , function(i){
return(c( strsplit(colnames(mypvals)[i],'\\.')[[1]],
sum(mypvals[,i] <0.05)))
}))
)
head(sm)
colnames(sm)=c('SOURCE' ,'TARGET' ,'count')
sm$count = as.numeric( sm$count )
sm
write.table(sm,file = 'count_network.txt',
sep = '\t',
quote = F,row.names = F)
这个时候得到 count_network.txt 文件，跟前面的 cellphonedb plot dot_plot 代码其实是一回事。
而且 这个 count_network.txt 文件里面的内容很容易可视化，代码如下所示:
library(reshape2)
sm_df =dcast(as.data.frame(sm),SOURCE~TARGET )
sm_df[is.na(sm_df)]=0
sm_df
rownames(sm_df) = sm_df[,1]
sm_df = sm_df[,-1]
p1=pheatmap::pheatmap(sm_df,display_numbers = T)
p2=pheatmap::pheatmap(log(sm_df+1),display_numbers = T)
library(patchwork)
library(cowplot)
library(ggplotify)
as.ggplot(p1) + as.ggplot(p2)
ggplot2::ggsave('heatmap.pdf',width = 10,height = 5)
很容易看到 单细胞亚群各自两两之间的受体配体基因对的数量 ，比如这个pbmc3k数据集里面的 FCGR3A_Mono 跟其它单细胞亚群的通讯数量就多很多 ：
单细胞亚群各自两两之间的受体配体基因对的数量
但是大家看文献，会发现绝大部分细胞通讯结果的展示都是类似于下面的蛋壳图:
来源于文献：《Single-cell transcriptome profiling of an adult human cell atlas of 15 major organs》，发表于2020年底，在genome biology杂志，链接是：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02210-0
上面的蛋壳图的右边其实就是 每个单细胞亚群跟其它单细胞亚群之间的连接数量，而上面的数值就是它与其它所有单细胞亚群连接数量的总和，还不如我们前面的热图展现的更加直观。
如果你确实需要把前面的单细胞亚群各自两两之间的受体配体基因对的数量热图，转换为蛋壳图，也是有成熟的代码，如下所示:
rm(list = ls())
library(psych)
library(qgraph)
library(igraph)
library(purrr) # map function
netf<- read.table("count_network.txt", header = T)
mynet <- netf
head(mynet)
net<- graph_from_data_frame(mynet)
plot(net)
allcolour=c('#a6cee3','#1f78b4','#b2df8a','#33a02c' )
karate_groups <- cluster_optimal(net)
coords <- layout_in_circle(net, order =
order(membership(karate_groups))) # 设置网络布局
E(net)$width <- E(net)$count/40 # 边点权重（粗细）
plot(net, edge.arrow.size=.1,
edge.curved=0,
vertex.color=allcolour,
vertex.frame.color="#555555",
vertex.label.color="black",
layout = coords,
vertex.label.cex=.7)
net2 <- net # 复制一份备用
for (i in 1: length(unique(mynet$SOURCE)) ){
E(net)[map(unique(mynet$SOURCE),function(x) {
get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))
})%>% unlist()]$color <- allcolour[i]
} # 这波操作谁有更好的解决方案？
plot(net, edge.arrow.size=.1,
edge.curved=0,
vertex.color=allcolour,
vertex.frame.color="#555555",
vertex.label.color="black",
layout = coords,
vertex.size = 30,
vertex.label.cex=.7)
dev.off()
length(unique(mynet$SOURCE))
# 查看需要绘制多少张图，以方便布局
# 需要人工肉眼，自己调整 mfrow=c(2,2)
par(mfrow=c(2,2), mar=c(.3,.3,.3,.3))
for (i in 1: length(unique(mynet$SOURCE)) ){
net1<-net2
E(net1)[map(unique(mynet$SOURCE),function(x) {
get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))
})%>% unlist()]$color <- allcolour[i]
plot(net1, edge.arrow.size=.1,
edge.curved=0.4,
vertex.color=allcolour,
vertex.frame.color="#555555",
vertex.label.color="black",
layout = coords,
vertex.label.cex=1,
vertex.size = 30)
}
如果是单细胞常规分析可以参考前面的例子：人人都能学会的单细胞聚类分群注释 ，我们演示了第一层次的分群。如果你对单细胞数据分析还没有基础认知，可以看基础10讲：
腾讯云开发者
