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

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

細胞通訊分析可以給我們一些細胞類群之間相互調控/交流的資訊，這種細胞之間的調控主要是通過受配體結合，傳遞訊號來實現的。不同的分化、疾病過程，可能存在特異的細胞通訊關係，因此闡明這些通訊關係至關重要。

CellPhoneDB配有詳實的受配體資料庫，其整合了此前的公共資料庫，還會手動矯正，以得到更加準確的受配體註釋。此外，針對受配體有多個亞基的情況，也進行了註釋。下面這張圖顯示了CellPhoneDB配有的資料庫包含多少種分泌蛋白和膜蛋白、蛋白質複合物、受配體關係，以及它們來源於什麼資料庫。

1. CellPhoneDB推斷細胞通訊的原理

在給定表達矩陣和細胞註釋之後，對於gene1-gene2這個互作關係，計算某一個clusterA裡面gene1的表達均值，計算另一個clusterB中gene2的表達均值，二者的均值為MEAN；在隨機更換細胞的label之後，依據新的標籤，計算“clusterA”裡面gene1的表達均值，"clusterB"中gene2的表達均值，再求一個平均值mean，這樣的過程重複多次，就可以得到一個mean的分佈，即null distribution。MEAN在這個分佈中所在的位置以及更極端的位置，構成的佔比，就是p值（p值的定義）。所以CellPhoneDB推測兩種細胞型別之間顯著富集的受配體關係，本質上還是基於一個細胞型別裡面的受體表達量，以及另一種細胞型別裡面的配體表達量。此外，如果某種關係無處不在（在所有細胞型別之間都很明顯），則找不出來。

此外還有幾個需要注意的地方：

  • 大樣本時會下采樣，只分析1/3的細胞
  • 多個亞基時考慮表達低的那一個亞基
  • 表達佔比達到一定閾值的基因才會被分析，預設是10%

2. 如何展示結果

這是原文獻給的視覺化例子，這裡有兩個地方需要注意：

  1. 右邊的熱圖表示細胞型別兩兩之間的相互作用的數量，我們可以看到沿著對角線，左右是對稱的，也就是A-B與B-A的互作數目是一樣的，為什麼會這樣？
  2. 左邊是具體受配體對，細胞對的互作氣泡圖，點的大小表示顯著水平，顏色則是The means of the average expression level of interacting molecule 1 in cluster 1 and interacting molecule 2 in cluster 2 注意到了嗎，說的是interacting molecule 1/2，而沒有說哪一個是受體哪一個是配體。

原因都和CellPhoneDB內建的gene-gene互作關係列表有關。CellPhoneDB區分不了受體還是配體，對於gene1-gene2，可以是gene1配體gene2受體，也可以是gene1受體gene2配體（如下圖）。我個人覺得也是由於這個原因，右邊那個熱圖為了說起來方便，才把不管做受體還是做配體的關係都算作是兩種細胞的互作關係，因此A-B和B-A在熱圖中的數值是一樣的（不然橫縱座標寫個interacting molecule，看到的人自然會問，這個分子是受體還是配體呢，加一起就省事了——都包含）。

這一點，github有提到：

也是這個原因，我看到文章如果用了CellPhoneDB的話，會留意它的圖，如果是用有向圖表示細胞群兩兩之間的關係數量，我會想這樣做合不合適（當然是不合適的）

3. 實際分析

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

3.1 輸入檔案的格式

註釋檔案
一共兩列，Cell列cell_type列，有列名；.csv, .txt字尾都行

表達檔案
normalize之後的矩陣，一般簡單相除normalize一下就行；.csv, .txt字尾都行

3.2 執行

軟體的安裝這裡就不講了，建立一個conda環境，pip install下載安裝就可以了

執行CellPhoneDB的主程式碼很簡單：

source /home/huangsiyuan/miniconda3/bin/activate cpdb

file_count=/home/huangsiyuan/cpdb/test_normat.txt
file_anno=/home/huangsiyuan/cpdb/test_anno.txt
outdir=/home/huangsiyuan/cpdb/test

if [ ! -d ${outdir} ]; then
  mkdir ${outdir}
fi

cellphonedb method statistical_analysis  \
            --counts-data hgnc_symbol \
            --output-path ${outdir} \
            --threshold 0.01 \ #Percentage of cells expressing the specific ligand or receptor
            --threads 10 \
            ${file_anno} ${file_count}

source /home/huangsiyuan/miniconda3/bin/deactivate cpdb

#如果細胞數太多，可以新增下采樣引數，預設只分析1/3的細胞
#--subsampling
#--subsampling-log true #對於沒有log轉化的資料，還要加這個引數

這一步之後在test資料夾裡面會生成4個檔案

deconvoluted.txt
means.txt
pvalues.txt
significant_means.txt

其中，

  • means.txt行是受配體pair，列是細胞pair，值為受體、配體在相應的cluster中表達均值的平均數；
  • pvalues.txt格式與means.txt類似，值為p值；
  • significant_means.txt格式和內容都與means.txt類似，不過僅保留了p值小於0.05的平均數。

4. 結果的視覺化

在這一步中，我一般只用到上述的means.txt和pvalues.txt檔案
我們還是先仿照文獻原文，畫出那兩張圖

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

pvalues=read.table("./test/pvalues.txt",header = T,sep = "\t",stringsAsFactors = F)
pvalues=pvalues[,12:dim(pvalues)[2]] #此時不關注前11列
statdf=as.data.frame(colSums(pvalues < 0.05)) #統計在某一種細胞pair的情況之下，顯著的受配體pair的數目；閾值可以自己選
colnames(statdf)=c("number")

#排在前面的分子定義為indexa；排在後面的分子定義為indexb
statdf$indexb=str_replace(rownames(statdf),"^.*\\.","")
statdf$indexa=str_replace(rownames(statdf),"\\..*$","")
#設定合適的細胞型別的順序
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=number))+geom_tile(color="white")+
  scale_fill_gradientn(colours = c("#4393C3","#ffdbba","#B2182B"),limits=c(0,20))+
  scale_x_discrete("cluster 1 produces molecule 1")+
  scale_y_discrete("cluster 2 produces molecule 2")+
  theme_minimal()+
  theme(
    axis.text.x.bottom = element_text(hjust = 1, vjust = NULL, angle = 45),
    panel.grid = element_blank()
  )
ggsave(filename = "interaction.num.1.pdf",device = "pdf",width = 12,height = 10,units = c("cm"))

這裡與文獻中圖不一致的地方是，我這個圖並不是關於對角線對稱的，因為我沒有將A-B，B-A的互作關係求和

舉個例子
在CellPhoneDB輸出的結果中，經統計，A-B有10個顯著的互作關係，B-A有20個顯著的互作關係【①】。然而A-B的互作其實包含A做配體8次,A做受體2次，B-A的互作其實包含B做配體19次,B做受體1次，所以嚴格來講，A和B兩種細胞互作，A做配體9次，B做配體21次【②】，這些資訊是CellPhoneDB給不了的。當然互作關係還是共計30次【③】。

換言之，文獻中對稱的圖給的資訊③，我上面那個圖給的資訊①，資訊②是不知道的（如果肉眼一個一個去看CellPhoneDB資料庫中gene1-gene2哪個是受體哪個是配體，還是可以統計出來的）。

因本文篇幅較長，餘下的視覺化部分將在下一篇展示，敬請期待~

參考文獻

[1] Efremova M, Vento-Tormo M, Teichmann S A, et al. CellPhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand–receptor complexes[J]. Nature protocols, 2020, 15(4): 1484-1506.

因水平有限，有錯誤的地方，歡迎批評指正！

相關文章