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

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

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

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

  1. library(tidyverse)
  2. library(RColorBrewer)
  3. library(scales)
  5. pvalues=read.table("./test/pvalues.txt",header = T,sep = "\t",stringsAsFactors = F)
  6. pvalues=pvalues[,12:dim(pvalues)[2]]
  7. statdf=as.data.frame(colSums(pvalues < 0.05))
  8. colnames(statdf)=c("number")
  10. statdf$indexb=str_replace(rownames(statdf),"^.*\\.","")
  11. statdf$indexa=str_replace(rownames(statdf),"\\..*$","")
  12. statdf$total_number=0
  14. for (i in 1:dim(statdf)[1]) {
  15.   tmp_indexb=statdf[i,"indexb"]
  16.   tmp_indexa=statdf[i,"indexa"]
  17.   if (tmp_indexa == tmp_indexb) {
  18.     statdf[i,"total_number"] = statdf[i,"number"]
  19.   } else {
  20.     statdf[i,"total_number"] = statdf[statdf$indexb==tmp_indexb & statdf$indexa==tmp_indexa,"number"]+
  21.       statdf[statdf$indexa==tmp_indexb & statdf$indexb==tmp_indexa,"number"]
  22.   }
  23. }
  25. rankname=sort(unique(statdf$indexa)) 
  26. statdf$indexa=factor(statdf$indexa,levels = rankname)
  27. statdf$indexb=factor(statdf$indexb,levels = rankname)
  29. statdf%>%ggplot(aes(x=indexa,y=indexb,fill=total_number))+geom_tile(color="white")+
  30.   scale_fill_gradientn(colours = c("#4393C3","#ffdbba","#B2182B"),limits=c(0,20))+
  31.   scale_x_discrete("cluster 1")+
  32.   scale_y_discrete("cluster 2")+
  33.   theme_minimal()+
  34.   theme(
  35.     axis.text.x.bottom = element_text(hjust = 1, vjust = NULL, angle = 45),
  36.     panel.grid = element_blank()
  37.   )
  38. ggsave(filename = "interaction.num.2.pdf",device = "pdf",width = 12,height = 10,units = c("cm"))

代码如下

  1. library(tidyverse)
  2. library(RColorBrewer)
  3. library(scales)
  4. library(igraph)
  6. pvalues=read.table("./test/pvalues.txt",header = T,sep = "\t",stringsAsFactors = F)
  7. pvalues=pvalues[,12:dim(pvalues)[2]]
  8. statdf=as.data.frame(colSums(pvalues < 0.05))
  9. colnames(statdf)=c("number")
  11. statdf$indexb=str_replace(rownames(statdf),"^.*\\.","")
  12. statdf$indexa=str_replace(rownames(statdf),"\\..*$","")
  13. rankname=sort(unique(statdf$indexa)) 
  15. A=c()
  16. B=c()
  17. C=c()
  18. remaining=rankname
  19. for (i in rankname[-6]) {
  20.   remaining=setdiff(remaining,i)
  21.   for (j in remaining) {
  22.     count=statdf[statdf$indexa == i & statdf$indexb == j,"number"]+
  23.       statdf[statdf$indexb == i & statdf$indexa == j,"number"]
  24.     A=append(A,i)
  25.     B=append(B,j)
  26.     C=append(C,count)
  27.   }
  28. }
  30. statdf2=data.frame(indexa=A,indexb=B,number=C)
  31. statdf2=statdf2 %>% rbind(statdf[statdf$indexa==statdf$indexb,c("indexa","indexb","number")])
  32. statdf2=statdf2[statdf2$number > 0,] #过滤掉值为0的观测
  34. #设置节点和连线的颜色
  35. color1=c("#8DD3C7", "#FDB462", "#B3DE69", "#FCCDE5", "#D9D9D9", "#BC80BD")
  36. names(color1)=rankname
  37. color2=colorRampPalette(brewer.pal(9, "Reds")[3:7])(20) #将颜色分成多少份，取决于互作关系数目的最大值
  38. names(color2)=1:20 #每一份颜色用对应的数字命名
  40. #做网络图
  41. ##下面的四行代码相对固定
  42. net <- graph_from_data_frame(statdf2[,c("indexa","indexb","number")])
  43. edge.start <- igraph::ends(net, es=igraph::E(net), names=FALSE)
  44. group <-  cluster_optimal(net)
  45. coords <- layout_in_circle(net, order = order(membership(group)))
  47. E(net)$width <- E(net)$number / 2 #将数值映射到连线的宽度，有时还需要微调，这里除以2就是这个目的
  48. E(net)$color <- color2[as.character(ifelse(E(net)$number > 20,20,E(net)$number))] #用前面设置好的颜色赋给连线，颜色深浅对应数值大小
  49. E(net)$label = E(net)$number #连线的标注
  50. E(net)$label.color <- "black" #连线标注的颜色
  51. V(net)$label.color <- "black" #节点标注的颜色
  52. V(net)$color <- color1[names(V(net))] #节点的填充颜色，前面已经设置了；V(net)返回节点信息
  54. #调整节点位置的线条角度
  55. ##如果没有这两行代码，节点位置的圆圈是向右的
  56. 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]))
  57. igraph::E(net)$loop.angle[which(edge.start[,2]==edge.start[,1])] <- loop.angle[edge.start[which(edge.start[,2]==edge.start[,1]),1]]
  59. #pdf("interaction.num.3.pdf",width = 6,height = 6)
  60. plot(net,
  61.      edge.arrow.size = 0, #连线不带箭头
  62.      edge.curved = 0, #连线不弯曲
  63.      vertex.frame.color = "black", #节点外框颜色
  64.      layout = coords,
  65.      vertex.label.cex = 1, #节点标注字体大小
  66.      vertex.size = 30) #节点大小
  67. #dev.off()

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

  1. source("CCC.bubble.R")
  2. CCC(
  3.   pfile="./test/pvalues.txt",
  4.   mfile="./test/means.txt",
  5.   #neg_log10_th= -log10(0.05),
  6.   #means_exp_log2_th=1,
  7.   #neg_log10_th2=3,
  8.   #means_exp_log2_th2=c(-4,6),
  9.   #notused.cell=c("Bcell","Gcell"),
  10.   #used.cell=c("Mcell"),
  11.   #cell.pair=c("Mcell.Scell","Mcell.NKcell","Mcell.Tcell","Scell.Mcell","NKcell.Mcell","Tcell.Mcell"),#这里是自定义的顺序，若是可选细胞对的子集，则只展示子集，若有交集则只展示交集；空值情况下，会根据可选细胞对自动排序
  12.   #gene.pair=c("MIF_TNFRSF14","FN1_aVb1 complex","EGFR_MIF")#作用同上
  13. )
  14. 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，如下：

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

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

版权声明：本文为TOP-Bio原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://www.cnblogs.com/TOP-Bio/p/15057637.html