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

2021 年 7 月 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. 如何展示結果

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

右邊的熱圖表示細胞類型兩兩之間的相互作用的數量，我們可以看到沿着對角線，左右是對稱的，也就是A-B與B-A的互作數目是一樣的，為什麼會這樣？
左邊是具體受配體對，細胞對的互作氣泡圖，點的大小表示顯著水平，顏色則是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.

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

Tags: 生物信息