細(xì)胞通訊分析可以給我們一些細(xì)胞類群之間相互調(diào)控/交流的信息，這種細(xì)胞之間的調(diào)控主要是通過受配體結(jié)合，傳遞信號(hào)來實(shí)現(xiàn)的。不同的分化、疾病過程，可能存在特異的細(xì)胞通訊關(guān)系，因此闡明這些通訊關(guān)系至關(guān)重要。

CellPhoneDB配有詳實(shí)的受配體數(shù)據(jù)庫，其整合了此前的公共數(shù)據(jù)庫，還會(huì)手動(dòng)矯正，以得到更加準(zhǔn)確的受配體注釋。此外，針對受配體有多個(gè)亞基的情況，也進(jìn)行了注釋。下面這張圖顯示了CellPhoneDB配有的數(shù)據(jù)庫包含多少種分泌蛋白和膜蛋白、蛋白質(zhì)復(fù)合物、受配體關(guān)系，以及它們來源于什么數(shù)據(jù)庫。

1. CellPhoneDB推斷細(xì)胞通訊的原理

在給定表達(dá)矩陣和細(xì)胞注釋之后，對于gene1-gene2這個(gè)互作關(guān)系，計(jì)算某一個(gè)clusterA里面gene1的表達(dá)均值，計(jì)算另一個(gè)clusterB中g(shù)ene2的表達(dá)均值，二者的均值為MEAN；在隨機(jī)更換細(xì)胞的label之后，依據(jù)新的標(biāo)簽，計(jì)算“clusterA”里面gene1的表達(dá)均值，"clusterB"中g(shù)ene2的表達(dá)均值，再求一個(gè)平均值mean，這樣的過程重復(fù)多次，就可以得到一個(gè)mean的分布，即null distribution。MEAN在這個(gè)分布中所在的位置以及更極端的位置，構(gòu)成的占比，就是p值（p值的定義）。所以CellPhoneDB推測兩種細(xì)胞類型之間顯著富集的受配體關(guān)系，本質(zhì)上還是基于一個(gè)細(xì)胞類型里面的受體表達(dá)量，以及另一種細(xì)胞類型里面的配體表達(dá)量。此外，如果某種關(guān)系無處不在（在所有細(xì)胞類型之間都很明顯），則找不出來。

此外還有幾個(gè)需要注意的地方：

• 大樣本時(shí)會(huì)下采樣，只分析1/3的細(xì)胞
• 多個(gè)亞基時(shí)考慮表達(dá)低的那一個(gè)亞基
• 表達(dá)占比達(dá)到一定閾值的基因才會(huì)被分析，默認(rèn)是10%

2. 如何展示結(jié)果

這是原文獻(xiàn)給的可視化例子，這里有兩個(gè)地方需要注意：

1. 右邊的熱圖表示細(xì)胞類型兩兩之間的相互作用的數(shù)量，我們可以看到沿著對角線，左右是對稱的，也就是A-B與B-A的互作數(shù)目是一樣的，為什么會(huì)這樣？
2. 左邊是具體受配體對，細(xì)胞對的互作氣泡圖，點(diǎn)的大小表示顯著水平，顏色則是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，而沒有說哪一個(gè)是受體哪一個(gè)是配體。

原因都和CellPhoneDB內(nèi)置的gene-gene互作關(guān)系列表有關(guān)。CellPhoneDB區(qū)分不了受體還是配體，對于gene1-gene2，可以是gene1配體gene2受體，也可以是gene1受體gene2配體（如下圖）。我個(gè)人覺得也是由于這個(gè)原因，右邊那個(gè)熱圖為了說起來方便，才把不管做受體還是做配體的關(guān)系都算作是兩種細(xì)胞的互作關(guān)系，因此A-B和B-A在熱圖中的數(shù)值是一樣的（不然橫縱坐標(biāo)寫個(gè)interacting molecule，看到的人自然會(huì)問，這個(gè)分子是受體還是配體呢，加一起就省事了——都包含）。

這一點(diǎn)，github有提到：

也是這個(gè)原因，我看到文章如果用了CellPhoneDB的話，會(huì)留意它的圖，如果是用有向圖表示細(xì)胞群兩兩之間的關(guān)系數(shù)量，我會(huì)想這樣做合不合適（當(dāng)然是不合適的）

3. 實(shí)際分析

公眾號(hào)后臺(tái)回復(fù)20210723獲取本次演示的測試數(shù)據(jù)，以及主要的可視化代碼。

3.1 輸入文件的格式

注釋文件
一共兩列，Cell列cell_type列，有列名；.csv, .txt后綴都行

表達(dá)文件
normalize之后的矩陣，一般簡單相除normalize一下就行；.csv, .txt后綴都行

3.2 運(yùn)行

軟件的安裝這里就不講了，創(chuàng)建一個(gè)conda環(huán)境，pip install下載安裝就可以了

運(yùn)行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

#如果細(xì)胞數(shù)太多，可以添加下采樣參數(shù)，默認(rèn)只分析1/3的細(xì)胞
#--subsampling
#--subsampling-log true #對于沒有l(wèi)og轉(zhuǎn)化的數(shù)據(jù)，還要加這個(gè)參數(shù)

這一步之后在test文件夾里面會(huì)生成4個(gè)文件

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

其中，

• means.txt行是受配體pair，列是細(xì)胞pair，值為受體、配體在相應(yīng)的cluster中表達(dá)均值的平均數(shù)；
• pvalues.txt格式與means.txt類似，值為p值；
• significant_means.txt格式和內(nèi)容都與means.txt類似，不過僅保留了p值小于0.05的平均數(shù)。

4. 結(jié)果的可視化

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

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

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

#排在前面的分子定義為indexa；排在后面的分子定義為indexb
statdf$indexb=str_replace(rownames(statdf),"^.*\\.","")
statdf$indexa=str_replace(rownames(statdf),"\\..*$","")
#設(shè)置合適的細(xì)胞類型的順序
rankname=sort(unique(statdf$indexa)) 
#轉(zhuǎn)成因子類型，畫圖時(shí)，圖形將按照預(yù)先設(shè)置的順序排列
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"))

這里與文獻(xiàn)中圖不一致的地方是，我這個(gè)圖并不是關(guān)于對角線對稱的，因?yàn)槲覜]有將A-B，B-A的互作關(guān)系求和

舉個(gè)例子
在CellPhoneDB輸出的結(jié)果中，經(jīng)統(tǒng)計(jì)，A-B有10個(gè)顯著的互作關(guān)系，B-A有20個(gè)顯著的互作關(guān)系【①】。然而A-B的互作其實(shí)包含A做配體8次,A做受體2次，B-A的互作其實(shí)包含B做配體19次,B做受體1次，所以嚴(yán)格來講，A和B兩種細(xì)胞互作，A做配體9次，B做配體21次【②】，這些信息是CellPhoneDB給不了的。當(dāng)然互作關(guān)系還是共計(jì)30次【③】。

換言之，文獻(xiàn)中對稱的圖給的信息③，我上面那個(gè)圖給的信息①，信息②是不知道的（如果肉眼一個(gè)一個(gè)去看CellPhoneDB數(shù)據(jù)庫中g(shù)ene1-gene2哪個(gè)是受體哪個(gè)是配體，還是可以統(tǒng)計(jì)出來的）。

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

參考文獻(xiàn)

[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.

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