單細胞分析實錄(19): 基於CellPhoneDB的細胞通訊分析及視覺化 (下篇)

TOP生物資訊發表於2021-07-25
視覺化

在上一篇帖子中，我介紹了CellPhoneDB的原理、實際操作，以及一些值得注意的地方。這一篇繼續細胞通訊分析的視覺化。

公眾號後臺回覆20210723獲取本次演示的測試資料，以及主要的視覺化程式碼。

所有的資料和結果檔案均已打包，下載後直接就能跑下面的程式碼畫圖。

下面的程式碼可以繪製對稱熱圖

（如果你不清楚為啥熱圖要沿著對角線對稱，可以看一下之前的推文）

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

pvalues=read.table("./test/pvalues.txt",header = T,sep = "\t",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 = "\t",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行。如果你也想用，歡迎轉發上一篇推文，截圖後發給公眾號後臺，留下郵箱，小編就會發給你哦。別怪小編套路呀，寫這兩篇帖子花了不少時間呢
因水平有限，有錯誤的地方，歡迎批評指正！

相關文章