單細胞分析實錄(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.txtpvalues.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.

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

相關文章