常用的細(xì)胞通訊軟件：
- CellphoneDB：是公開的人工校正的凭疮，儲(chǔ)存受體嗤谚、配體以及兩種相互作用的數(shù)據(jù)庫(kù)滔蝉。此外穆役，還考慮了結(jié)構(gòu)組成，能夠描述異構(gòu)復(fù)合物啡专。（配體-受體+多聚體）
- iTALK：通過平均表達(dá)量方式，篩選高表達(dá)的胚體和受體博其，根據(jù)結(jié)果作圈圖套才。（配體-受體）
- CellChat：CellChat將細(xì)胞的基因表達(dá)數(shù)據(jù)作為輸入，并結(jié)合配體受體及其輔助因子的相互作用來模擬細(xì)胞間通訊慕淡。（配體-受體+多聚體+輔因子）
- NicheNet // NicheNet多樣本分析 // 目標(biāo)基因的配體和靶基因活性預(yù)測(cè)：通過將相互作用細(xì)胞的表達(dá)數(shù)據(jù)與信號(hào)和基因調(diào)控網(wǎng)絡(luò)的先驗(yàn)知識(shí)相結(jié)合來預(yù)測(cè)相互作用細(xì)胞之間的配體-靶標(biāo)聯(lián)系的方法背伴。( 配體-受體+信號(hào)通路）
附：NicheNet使用的常見問題匯總其它細(xì)胞互作軟件還包括
Celltalker，SingleCellSignalR峰髓，scTensor和SoptSC（這幾個(gè)也是基于配體-受體相互作用）
1. CellPhoneDB介紹
Nature protocol原文
CellPhoneDB是包含配體傻寂、受體及其相互作用的數(shù)據(jù)庫(kù)，可以對(duì)細(xì)胞間的通訊分子進(jìn)行全面携兵、系統(tǒng)的分析疾掰，研究不同細(xì)胞類型之間的相互交流及通訊網(wǎng)絡(luò)。
CellPhoneDB 是目前使用最廣泛的細(xì)胞通訊分析軟件徐紧，CellPhoneDB的配體-受體數(shù)據(jù)庫(kù)集成于UniProt静檬、Ensembl、PDB并级、IUPHAR等拂檩，共存儲(chǔ)978種蛋白質(zhì)：其中501種是分泌蛋白，而585種是膜蛋白嘲碧，這些蛋白質(zhì)參與1,396種相互作用稻励。CellPhoneDB數(shù)據(jù)庫(kù)還考慮了配體和受體的亞基結(jié)構(gòu)，準(zhǔn)確地表示了異聚體復(fù)合物愈涩，有466個(gè)是異聚體望抽，對(duì)于準(zhǔn)確研究多亞基復(fù)合物介導(dǎo)的細(xì)胞通訊很關(guān)鍵。CellPhoneDB有474種相互作用涉及分泌蛋白履婉，而490種相互作用涉及膜蛋白煤篙，總共有250個(gè)涉及整合素的相互作用。
數(shù)據(jù)庫(kù)鏈接：https://www.cellphonedb.org/ppi-resources毁腿。
目前支持的物種：人（也可通過人和鼠的同源基因比對(duì)應(yīng)用于鼠）
2. CellPhoneDB原理
參考：https://www.cellphonedb.org/explore-sc-rna-seq
- CellphoneDB基于一種細(xì)胞類型的受體和另一種細(xì)胞類型的配體的表達(dá)辑奈，得出兩種細(xì)胞群之間豐富的受體-配體相互作用。對(duì)于細(xì)胞群所表達(dá)的基因狸棍，計(jì)算表達(dá)該基因的細(xì)胞百分比和基因表達(dá)平均值身害。若該基因只在該群中10%及以下的細(xì)胞中表達(dá)（默認(rèn)值為0.1）味悄，則被移除草戈。
- 將所有細(xì)胞的簇標(biāo)簽隨機(jī)排列1000次（可選值），并確定簇中受體平均表達(dá)水平和相互作用簇中配體平均表達(dá)水平的平均值侍瑟。對(duì)于兩種細(xì)胞類型之間每對(duì)比較中的每一個(gè)受體-配體對(duì)唐片，這產(chǎn)生了一個(gè)零分布（null distribution丙猬，亦稱伯努利分布、兩點(diǎn)分布费韭，指的是對(duì)于隨機(jī)變量X有, 參數(shù)為p(0<p<1)茧球，如果它分別以概率p和1-p取1和0為值。EX= p,DX=p(1-p)星持。伯努利試驗(yàn)成功的次數(shù)服從伯努利分布,參數(shù)p是試驗(yàn)成功的概率）抢埋。
- 通過計(jì)算等于或高于實(shí)際平均值的平均值的比例，可以得到了一個(gè)P值督暂，表示給定受體-配體復(fù)合物細(xì)胞類型特異性的可能性揪垄。換句話說，如果觀察到的平均值在前5%逻翁，則相互作用的P值為0.05饥努。
- 根據(jù)兩種細(xì)胞類型中富集到的顯著的受體-配體對(duì)的數(shù)量，對(duì)細(xì)胞類型之間高度特異性的相互作用進(jìn)行排序八回，以便手動(dòng)篩選出生物學(xué)相關(guān)的相互作用關(guān)系酷愧。
3. CellPhoneDB下載和使用
3.1 CellPhoneDB下載和主要功能介紹
安裝：
CellPhoneDB軟件倉(cāng)庫(kù)地址：https://github.com/Teichlab/cellphonedb
PIP安裝：pip install cellphonedb（python版本需要python>3.5）
CellPhoneDB安裝完成后，可以在終端測(cè)試是否成功缠诅。CellPhoneDB主要分成method溶浴、plot、query和database4個(gè)模塊滴铅。
我們主要進(jìn)行分析的是method（分析模塊）和plot模塊（畫圖模塊）戳葵。
query是進(jìn)行數(shù)據(jù)查詢的模塊，查詢基因有哪些互作結(jié)果（get_interaction_gene結(jié)果為數(shù)據(jù)庫(kù)涉及到的基因汉匙，find_interactions_by_element可以找到特定基因的受體配體作用數(shù)據(jù)對(duì)）拱烁。
database是輸入數(shù)據(jù)庫(kù)，一般默認(rèn)可以不輸入噩翠，直接使用即可戏自。
- 四個(gè)重要的函數(shù)：
-
核心可選參數(shù)：
示例：
#（1）設(shè)置迭代和線程的數(shù)量
cellphonedb method statistical_analysis yourmetafile.txt yourcountsfile.txt --iterations=10 --threads=2
#（2）子項(xiàng)目文件夾
cellphonedb method analysis yourmetafile.txt yourcountsfile.txt --project-name=new_project
# （3）設(shè)置輸出路徑
mkdir custom_folder
cellphonedb method statistical_analysis yourmetafile.txt yourcountsfile.txt --output-path=custom_folder
#（4） 二次抽樣
cellphonedb method analysis yourmetafile.txt yourcountsfile.txt --subsampling --subsampling-log false --subsampling-num-cells 3000
3.2 數(shù)據(jù)準(zhǔn)備：注釋好的pbmc3k（數(shù)據(jù)下載和準(zhǔn)備見monocle）
pbmc3k <- readRDS("pbmc.rds")
由于細(xì)胞類型已經(jīng)注釋好，接下來準(zhǔn)備cellphonedb的文件：表達(dá)譜文件cellphonedb_count.txt和細(xì)胞分組注釋文件cellphonedb_meta.txt伤锚。
write.table(as.matrix(pbmc3k@assays$RNA@data), 'cellphonedb_count.txt', sep='\t', quote=F)
meta_data <- cbind(rownames(pbmc3k@meta.data), pbmc3k@meta.data[,'cell_type', drop=F])
meta_data <- as.matrix(meta_data)
meta_data[is.na(meta_data)] = "Unkown" # 細(xì)胞類型中不能有NA
write.table(meta_data, 'cellphonedb_meta.txt', sep='\t', quote=F, row.names=F)
3.3 終端運(yùn)行cellphonedb：
安裝好cellphonedb并加入環(huán)境變量后擅笔，使用終端運(yùn)行如下代碼：
cellphonedb method statistical_analysis cellphonedb_meta.txt cellphonedb_count.txt --counts-data=gene_name
#如果我們count的基因是基因名格式，需要添加參數(shù)--counts-data=gene_name屯援，如果行名為ensemble名稱的話猛们，可以不添加這個(gè)參數(shù)，使用默認(rèn)值即可狞洋。
終端運(yùn)行完上述代碼弯淘，在工作目錄下會(huì)得到一個(gè)out文件夾，里面有四個(gè)txt：
deconvoluted.txt：基因在亞群中的平均表達(dá)量
mean.txt：每對(duì)受體-配體的平均表達(dá)量
pvalues.txt：每對(duì)受體-配體的p值
significant_means.txt：每對(duì)受體-配體顯著性結(jié)果的平均表達(dá)量值
#cellphonedb 自己的繪圖
cellphonedb plot dot_plot
cellphonedb plot heatmap_plot cellphonedb_meta.txt
tree out
#使用tree來查看cellphonedb生成的out文件夾下的文件吉懊。如果沒有tree庐橙，可以使用brew install tree假勿。
out
├── count_network.txt ## 繪制網(wǎng)絡(luò)圖文件。注意這個(gè)文件是在cellphonedb plot中生成的态鳖，所以上面兩個(gè)plot不要漏掉了转培。
├── deconvoluted.txt
├── heatmap_count.pdf
├── heatmap_log_count.pdf
├── interaction_count.txt
├── means.txt ## 繪制點(diǎn)圖文件
├── plot.pdf
├── pvalues.txt ## 繪制點(diǎn)圖文件
└── significant_means.txt
得到3個(gè)圖：
- plot.pdf
CellPhoneDB氣泡圖：每一列為兩個(gè)細(xì)胞亞群（如：B|DC T），每一行為一對(duì)受體配體名稱（如：CD2_CD58）浆竭，顏色代表兩個(gè)亞群這兩個(gè)基因的平均表達(dá)量高低浸须，越紅表示表達(dá)越高，氣泡大小代表P值的-log10值邦泄，氣泡越大羽戒，說明其越具顯著性。
- heatmap_count.pdf和heatmap_log_count.pdf
CellPhoneDB Heatmap：下面左圖和右圖其實(shí)都是對(duì)interaction_count.txt結(jié)果的展示虎韵，左圖為亞群之間相互作用受體配體數(shù)目的熱圖易稠，右圖為將這個(gè)數(shù)目取了log10。
上述文件和圖形的解讀可以參考官網(wǎng)：https://www.cellphonedb.org/documentation
3.4 使用R對(duì)結(jié)果進(jìn)行可視化（主要是網(wǎng)絡(luò)圖和點(diǎn)圖）
pbmc='/practice/cellphonedb/out/' #outs下的文件放在這里了
library(psych)
library(qgraph)
library(igraph)
library(tidyverse)
netf<- "count_network.txt"
mynet <- read.delim(paste0(pbmc,"count_network.txt"), check.names = FALSE) #讀取數(shù)據(jù)
table(mynet$count)
mynet %>% filter(count>0) -> mynet # 過濾count為0的數(shù)據(jù)（有零會(huì)報(bào)錯(cuò)）
head(mynet)
net<- graph_from_data_frame(mynet) #構(gòu)建net對(duì)象
plot(net)
為了給這個(gè)網(wǎng)絡(luò)圖的邊點(diǎn)mapping上不同的屬性我們引入一串顏色
allcolour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB",
"#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",
"#808000","#FF00FF","#FA8072","#7B68EE","#9400D3",
"#800080","#A0522D","#D2B48C","#D2691E","#87CEEB",
"#40E0D0","#5F9EA0","#FF1493",
"#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4",
"#32CD32","#F0E68C","#FFFFE0","#EE82EE","#FF6347",
"#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")
karate_groups <- cluster_optimal(net) #統(tǒng)計(jì)每個(gè)端點(diǎn)的和
coords <- layout_in_circle(net, order =
order(membership(karate_groups))) # 設(shè)置網(wǎng)絡(luò)布局
E(net)$width <- E(net)$count/10 #根據(jù)count值設(shè)置邊的寬
plot(net, edge.arrow.size=.1,
edge.curved=0,
vertex.color=allcolour,
vertex.frame.color="#555555",
vertex.label.color="black",
layout = coords,
vertex.label.cex=.7)
但是邊的顏色和點(diǎn)的顏色還是對(duì)應(yīng)不上包蓝，我們來修改一番邊的屬性驶社。
net2 <- net # 復(fù)制一份備用
for (i in 1: length(unique(mynet$SOURCE)) ){ #配置發(fā)出端的顏色
E(net)[map(unique(mynet$SOURCE),function(x) {
get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))
})%>% unlist()]$color <- allcolour[i]
} # 這波操作誰(shuí)有更好的解決方案？
plot(net, edge.arrow.size=.1,
edge.curved=0,
vertex.color=allcolour,
vertex.frame.color="#555555",
vertex.label.color="black",
layout = coords,
vertex.label.cex=.7)
我們觀察到由于箭頭是雙向的测萎，所以兩個(gè)細(xì)胞類型之間的線條會(huì)被后來畫上去的覆蓋掉亡电，這里我們把線條做成曲線：
plot(net, edge.arrow.size=.1, #箭頭大小設(shè)置為0.1
edge.curved=0.2, # 只是調(diào)了曲率這個(gè)參數(shù)
vertex.color=allcolour,
vertex.frame.color="#555555", #圓圈顏色
vertex.label.color="black", #標(biāo)簽顏色
layout = coords, #網(wǎng)絡(luò)布局位點(diǎn)
vertex.label.cex=.7#標(biāo)簽大小)
接下來，我們來繪制第二組類型貝殼一樣的網(wǎng)絡(luò)圖硅瞧。
dev.off()
length(unique(mynet$SOURCE)) # 查看需要繪制多少?gòu)垐D份乒，以方便布局
par(mfrow=c(2,5), mar=c(.3,.3,.3,.3))
for (i in 1: length(unique(mynet$SOURCE)) ){
net1<-net2
E(net1)[map(unique(mynet$SOURCE),function(x) {
get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))
})%>% unlist()]$color <- allcolour[i]
plot(net1, edge.arrow.size=.1,
edge.curved=0.4,
vertex.color=allcolour,
vertex.frame.color="#555555",
vertex.label.color="black",
layout = coords,
vertex.label.cex=1)
}
但是，細(xì)胞間通訊的頻數(shù)（count）并沒有繪制在圖上腕唧，略顯不足或辖，這就補(bǔ)上吧。
dev.off()
length(unique(mynet$SOURCE))
par(mfrow=c(2,5), mar=c(.3,.3,.3,.3))
for (i in 1: length(unique(mynet$SOURCE)) ){
net1<-net2
E(net1)$count <- ""
E(net1)[map(unique(mynet$SOURCE),function(x) {
get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))
})%>% unlist()]$count <- E(net2)[map(unique(mynet$SOURCE),function(x) {
get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))
})%>% unlist()]$count # 故技重施
E(net1)[map(unique(mynet$SOURCE),function(x) {
get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))
})%>% unlist()]$color <- allcolour[i]
plot(net1, edge.arrow.size=.1,
edge.curved=0.4,
edge.label = E(net1)$count, # 繪制邊的權(quán)重
vertex.color=allcolour,
vertex.frame.color="#555555",
vertex.label.color="black",
layout = coords,
vertex.label.cex=1
)
}
找兩條邊驗(yàn)證一下對(duì)應(yīng)關(guān)系枣接。
mynet %>% filter(SOURCE == 'B',TARGET == 'DC')
SOURCE TARGET count
1 B DC 13
mynet %>% filter(SOURCE == 'Memory CD4 T',TARGET == 'B')
SOURCE TARGET count
1 Memory CD4 T B 15
用ggplot2 繪制點(diǎn)圖颂暇，關(guān)鍵是把兩張表合并到一起。
mypvals <- read.delim(paste0(pbmc,"pvalues.txt"), check.names = FALSE)
mymeans <- read.delim(paste0(pbmc,"means.txt"), check.names = FALSE)
# 這些基因list很有意思啊但惶，建議保存
chemokines <- grep("^CXC|CCL|CCR|CX3|XCL|XCR", mymeans$interacting_pair,value = T)
chemokines <- grep("^CXC|CCL|CCR|CX3|XCL|XCR", mymeans$interacting_pair,value = T)
th1 <- grep("IL2|IL12|IL18|IL27|IFNG|IL10|TNF$|TNF |LTA|LTB|STAT1|CCR5|CXCR3|IL12RB1|IFNGR1|TBX21|STAT4",
mymeans$interacting_pair,value = T)
th2 <- grep("IL4|IL5|IL25|IL10|IL13|AREG|STAT6|GATA3|IL4R",
mymeans$interacting_pair,value = T)
th17 <- grep("IL21|IL22|IL24|IL26|IL17A|IL17A|IL17F|IL17RA|IL10|RORC|RORA|STAT3|CCR4|CCR6|IL23RA|TGFB",
mymeans$interacting_pair,value = T)
treg <- grep("IL35|IL10|FOXP3|IL2RA|TGFB", mymeans$interacting_pair,value = T)
costimulatory <- grep("CD86|CD80|CD48|LILRB2|LILRB4|TNF|CD2|ICAM|SLAM|LT[AB]|NECTIN2|CD40|CD70|CD27|CD28|CD58|TSLP|PVR|CD44|CD55|CD[1-9]",
mymeans$interacting_pair,value = T)
coinhibitory <- grep("SIRP|CD47|ICOS|TIGIT|CTLA4|PDCD1|CD274|LAG3|HAVCR|VSIR",
mymeans$interacting_pair,value = T)
niche <- grep("CSF", mymeans$interacting_pair,value = T)
mymeans %>% dplyr::filter(interacting_pair %in% costimulatory)%>%
dplyr::select("interacting_pair",starts_with("NK"),ends_with("NK")) %>%
reshape2::melt() -> meansdf
colnames(meansdf)<- c("interacting_pair","CC","means")
mypvals %>% dplyr::filter(interacting_pair %in% costimulatory)%>%
dplyr::select("interacting_pair",starts_with("NK"),ends_with("NK"))%>%
reshape2::melt()-> pvalsdf
colnames(pvalsdf)<- c("interacting_pair","CC","pvals")
pvalsdf$joinlab<- paste0(pvalsdf$interacting_pair,"_",pvalsdf$CC)
meansdf$joinlab<- paste0(meansdf$interacting_pair,"_",meansdf$CC)
pldf <- merge(pvalsdf,meansdf,by = "joinlab")
summary((filter(pldf,means >1))$means)
pldf%>% filter(means >1) %>%
ggplot(aes(CC.x,interacting_pair.x) )+
geom_point(aes(color=means,size=-log10(pvals+0.0001)) ) +
scale_size_continuous(range = c(1,3))+
scale_color_gradient2(high="red",mid = "yellow",low ="darkblue",midpoint = 25 )+ theme_bw()+
theme(axis.text.x = element_text(angle = -45,hjust = -0.1,vjust = 0.8))
