单细胞分析实录(19): 基于CellPhoneDB的细胞通讯分析及可视化 (下篇)

在上一篇帖子中，我介绍了CellPhoneDB的原理、实际操作，以及一些值得注意的地方。这一篇继续细胞通讯分析的可视化。

公众号后台回复20210723获取本次演示的测试数据，以及主要的可视化代码。

所有的数据和结果文件均已打包，下载后直接就能跑下面的代码画图。

下面的代码可以绘制对称热图

（如果你不清楚为啥热图要沿着对角线对称，可以看一下之前的推文）

library(tidyverse)
library(RColorBrewer)
library(scales)

pvalues=read.table("./test/pvalues.txt",header = T,sep = "	",stringsAsFactors = F)
pvalues=pvalues[,12:dim(pvalues)[2]]
statdf=as.data.frame(colSums(pvalues < 0.05))
colnames(statdf)=c("number")

statdf$indexb=str_replace(rownames(statdf),"^.*\.","")
statdf$indexa=str_replace(rownames(statdf),"\..*$","")
statdf$total_number=0

for (i in 1:dim(statdf)[1]) {
  tmp_indexb=statdf[i,"indexb"]
  tmp_indexa=statdf[i,"indexa"]
  if (tmp_indexa == tmp_indexb) {
    statdf[i,"total_number"] = statdf[i,"number"]
  } else {
    statdf[i,"total_number"] = statdf[statdf$indexb==tmp_indexb & statdf$indexa==tmp_indexa,"number"]+
      statdf[statdf$indexa==tmp_indexb & statdf$indexb==tmp_indexa,"number"]
  }
}

rankname=sort(unique(statdf$indexa)) 
statdf$indexa=factor(statdf$indexa,levels = rankname)
statdf$indexb=factor(statdf$indexb,levels = rankname)

statdf%>%ggplot(aes(x=indexa,y=indexb,fill=total_number))+geom_tile(color="white")+
  scale_fill_gradientn(colours = c("#4393C3","#ffdbba","#B2182B"),limits=c(0,20))+
  scale_x_discrete("cluster 1")+
  scale_y_discrete("cluster 2")+
  theme_minimal()+
  theme(
    axis.text.x.bottom = element_text(hjust = 1, vjust = NULL, angle = 45),
    panel.grid = element_blank()
  )
ggsave(filename = "interaction.num.2.pdf",device = "pdf",width = 12,height = 10,units = c("cm"))

还可以用网络图表示互作关系的数量

代码如下

library(tidyverse)
library(RColorBrewer)
library(scales)
library(igraph)

pvalues=read.table("./test/pvalues.txt",header = T,sep = "	",stringsAsFactors = F)
pvalues=pvalues[,12:dim(pvalues)[2]]
statdf=as.data.frame(colSums(pvalues < 0.05))
colnames(statdf)=c("number")

statdf$indexb=str_replace(rownames(statdf),"^.*\.","")
statdf$indexa=str_replace(rownames(statdf),"\..*$","")
rankname=sort(unique(statdf$indexa)) 

A=c()
B=c()
C=c()
remaining=rankname
for (i in rankname[-6]) {
  remaining=setdiff(remaining,i)
  for (j in remaining) {
    count=statdf[statdf$indexa == i & statdf$indexb == j,"number"]+
      statdf[statdf$indexb == i & statdf$indexa == j,"number"]
    A=append(A,i)
    B=append(B,j)
    C=append(C,count)
  }
}

statdf2=data.frame(indexa=A,indexb=B,number=C)
statdf2=statdf2 %>% rbind(statdf[statdf$indexa==statdf$indexb,c("indexa","indexb","number")])
statdf2=statdf2[statdf2$number > 0,] #过滤掉值为0的观测

#设置节点和连线的颜色
color1=c("#8DD3C7", "#FDB462", "#B3DE69", "#FCCDE5", "#D9D9D9", "#BC80BD")
names(color1)=rankname
color2=colorRampPalette(brewer.pal(9, "Reds")[3:7])(20) #将颜色分成多少份，取决于互作关系数目的最大值
names(color2)=1:20 #每一份颜色用对应的数字命名

#做网络图
##下面的四行代码相对固定
net <- graph_from_data_frame(statdf2[,c("indexa","indexb","number")])
edge.start <- igraph::ends(net, es=igraph::E(net), names=FALSE)
group <-  cluster_optimal(net)
coords <- layout_in_circle(net, order = order(membership(group)))

E(net)$width <- E(net)$number / 2 #将数值映射到连线的宽度，有时还需要微调，这里除以2就是这个目的
E(net)$color <- color2[as.character(ifelse(E(net)$number > 20,20,E(net)$number))] #用前面设置好的颜色赋给连线，颜色深浅对应数值大小
E(net)$label = E(net)$number #连线的标注
E(net)$label.color <- "black" #连线标注的颜色
V(net)$label.color <- "black" #节点标注的颜色
V(net)$color <- color1[names(V(net))] #节点的填充颜色，前面已经设置了；V(net)返回节点信息

#调整节点位置的线条角度
##如果没有这两行代码，节点位置的圆圈是向右的
loop.angle<-ifelse(coords[igraph::V(net),1]>0,-atan(coords[igraph::V(net),2]/coords[igraph::V(net),1]),pi-atan(coords[igraph::V(net),2]/coords[igraph::V(net),1]))
igraph::E(net)$loop.angle[which(edge.start[,2]==edge.start[,1])] <- loop.angle[edge.start[which(edge.start[,2]==edge.start[,1]),1]]

#pdf("interaction.num.3.pdf",width = 6,height = 6)
plot(net,
     edge.arrow.size = 0, #连线不带箭头
     edge.curved = 0, #连线不弯曲
     vertex.frame.color = "black", #节点外框颜色
     layout = coords,
     vertex.label.cex = 1, #节点标注字体大小
     vertex.size = 30) #节点大小
#dev.off()

气泡图——具体的互作关系

以上几种图，只是用来展示数量，具体的两种细胞之间的互作关系可以用如下的代码展示：

source("CCC.bubble.R")
CCC(
  pfile="./test/pvalues.txt",
  mfile="./test/means.txt",
  #neg_log10_th= -log10(0.05),
  #means_exp_log2_th=1,
  #neg_log10_th2=3,
  #means_exp_log2_th2=c(-4,6),
  #notused.cell=c("Bcell","Gcell"),
  #used.cell=c("Mcell"),
  #cell.pair=c("Mcell.Scell","Mcell.NKcell","Mcell.Tcell","Scell.Mcell","NKcell.Mcell","Tcell.Mcell"),#这里是自定义的顺序，若是可选细胞对的子集，则只展示子集，若有交集则只展示交集；空值情况下，会根据可选细胞对自动排序
  #gene.pair=c("MIF_TNFRSF14","FN1_aVb1 complex","EGFR_MIF")#作用同上
)
ggsave(filename = "interaction.detail.1.pdf",device = "pdf",width =20,height = 12,units = "cm")

参数解释：

• neg_log10_th和means_exp_log2_th两个参数用来筛选显著的互作关系；
• neg_log10_th2和means_exp_log2_th2两个参数用来限定最终气泡图的数值范围；
• notused.cell不包含的细胞类型；
• used.cell必须包含的细胞类型；
• cell.pair必须包含的细胞pair，以及它们的顺序；
• gene.pair必须包含的基因pair，以及它们的顺序。

后面四个参数在细化气泡图的时候，很有用。

我们先不加额外的参数，查看全部的互作关系

随后再细化，指定需要展示的细胞类型和gene pair，如下：

CCC(
  pfile="./test/pvalues.txt",
  mfile="./test/means.txt",
  cell.pair=c("Mcell.Scell","Mcell.NKcell","Mcell.Tcell","Scell.Mcell","NKcell.Mcell","Tcell.Mcell"),#这里是自定义的顺序，若是可选细胞对的子集，则只展示子集，若有交集则只展示交集；空值情况下，会根据可选细胞对自动排序
  gene.pair=c("MIF_TNFRSF14","FN1_aVb1 complex","EGFR_MIF")#作用同上
)
ggsave(filename = "interaction.detail.2.pdf",device = "pdf",width =16,height = 10,units = "cm")

最后那个CCC( )函数是小编写的，小编觉得还挺好用的。并不复杂，也才120行。如果你也想用，欢迎转发上一篇推文，截图后发给公众号后台，留下邮箱，小编就会发给你哦。别怪小编套路呀，写这两篇帖子花了不少时间呢
因水平有限，有错误的地方，欢迎批评指正！