參考鏈接：
什么是細(xì)胞通訊
不同細(xì)胞類型和組織中的細(xì)胞之間相互作用(CCI)坐榆，主要是通過(guò)各種分子包括離子燎斩、代謝物棵癣、整聯(lián)蛋白钠署、受體憋肖、連接蛋白蒋得、結(jié)構(gòu)蛋白乙墙、配體和細(xì)胞外基質(zhì)的分泌蛋白來(lái)實(shí)現(xiàn)的。一些分子如激素咸产、生長(zhǎng)因子、趨化因子仲闽、細(xì)胞因子和神經(jīng)遞質(zhì)等配體介導(dǎo)細(xì)胞之間的信息交流(CCC)脑溢。 人們可以從基因表達(dá)中推斷出CCI和CCC，尤其是在單細(xì)胞轉(zhuǎn)錄組測(cè)序之后赖欣，以RNA為代表的細(xì)胞信息之間的相互作用更能夠?qū)ζ溥M(jìn)行解析
主要是流程就是根據(jù)數(shù)據(jù)庫(kù)記錄的參與細(xì)胞間通訊的相互作用蛋白列表屑彻，篩選表達(dá)矩陣后，看配體和受體表達(dá)量根據(jù)某個(gè)評(píng)分函數(shù)的交流得分顶吮，多個(gè)配體和受體聚合成為細(xì)胞與細(xì)胞之間的總體交互狀態(tài)社牲，最后通過(guò)網(wǎng)絡(luò)可視化展現(xiàn)出來(lái)
細(xì)胞通訊的方式
- 自分泌
自分泌信號(hào)轉(zhuǎn)導(dǎo)是指細(xì)胞內(nèi)通訊，細(xì)胞分泌配體云矫，這些配體用于通過(guò)同源受體誘導(dǎo)同一細(xì)胞上表達(dá)的那些分子的細(xì)胞應(yīng)答旁分泌- 旁分泌
旁分泌細(xì)胞間的通訊不需要細(xì)胞間的接觸，而是取決于信號(hào)分子在分泌后從一個(gè)細(xì)胞擴(kuò)散到另一個(gè)細(xì)胞- 近分泌
即依賴于接觸的細(xì)胞間通訊依賴于間隙連接或膜納米管等其他結(jié)構(gòu)汗菜，使信號(hào)分子直接在細(xì)胞之間傳遞让禀，而不會(huì)分泌到細(xì)胞外- 內(nèi)分泌
內(nèi)分泌細(xì)胞間的通訊代表細(xì)胞間的通訊，信號(hào)分子被分泌并通過(guò)諸如血漿的細(xì)胞外液傳播很長(zhǎng)一段距離
分析原理
在給定表達(dá)矩陣和細(xì)胞注釋之后陨界，對(duì)于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录淡，這樣的過(guò)程重復(fù)多次捌木，就可以得到一個(gè)mean的分布，即null distribution嫉戚。MEAN在這個(gè)分布中所在的位置以及更極端的位置刨裆，構(gòu)成的占比澈圈，就是p值（p值的定義）。所以CellPhoneDB推測(cè)兩種細(xì)胞類型之間顯著富集的受配體關(guān)系帆啃，本質(zhì)上還是基于一個(gè)細(xì)胞類型里面的受體表達(dá)量瞬女，以及另一種細(xì)胞類型里面的配體表達(dá)量。此外努潘，如果某種關(guān)系無(wú)處不在（在所有細(xì)胞類型之間都很明顯）诽偷，則找不出來(lái)
代碼實(shí)戰(zhàn)
rm(list = ls()) # 清除所有對(duì)象
library(Matrix)
# 將工作空間設(shè)置到此處,后面用到的文件都在此空間中（根據(jù)自己需求設(shè)置）
setwd("D:/personal_file/Bioinfomatics Study/data")
getwd()
# [1] "D:/personal_file/Bioinfomatics Study/data"
# 讀入數(shù)據(jù)并制作表達(dá)矩陣
matrix_data <- Matrix::readMM("matrix.txt")
print(object.size(matrix_data), unit="MB")
new_matrix_data <- Matrix::as.matrix(matrix_data)
col <- read.table("barcodes.txt")
col <- as.character(col$V1)
row <- read.table("features.txt")
row <- as.character(row$V1)
row <- toupper(row)
dimnames(new_matrix_data)=list(row,col)
exp <- as.data.frame(new_matrix_data) # 轉(zhuǎn)化為數(shù)據(jù)框便于操作
print(object.size(exp), unit="GB")
rm(new_matrix_data,matrix_data)
gc()
讀入數(shù)據(jù)并整理，最終得到的是exp表達(dá)矩陣慈俯，表達(dá)矩陣的行代表基因渤刃，列代表細(xì)胞樣本。因?yàn)樽x入后就是小數(shù)點(diǎn)形式的贴膘，表示該表達(dá)矩陣已normalization了
# 篩選出配受體基因并整理exp
LRgene <- read.table("geneOfpairs.txt", col.names = "Gene")
exp <- exp[rownames(exp) %in% LRgene$Gene,]
# 計(jì)算NK卖子、Tumor和MDSC類型細(xì)胞的基因表達(dá)水平
ann <- read.table("mdsc_tumor_nk_celltype.txt", head=T, row.names = 1)
MDSC_c <- ann$cell[ann$celltype=="MDSC"] #MDSC類型細(xì)胞
NK_c <- ann$cell[ann$celltype=="NK"]
Tumor_c <- ann$cell[ann$celltype=="Tumor"]
MDSC_Gene_Mean <- rowMeans(exp[MDSC_c]) #計(jì)算MDSC類型細(xì)胞的基因表達(dá)水平
NK_Gene_Mean <- rowMeans(exp[NK_c])
Tumor_Gene_Mean <- rowMeans(exp[Tumor_c])
geneMean <- cbind(row.names(exp), MDSC_Gene_Mean, NK_Gene_Mean, Tumor_Gene_Mean)
rownames(geneMean) <- seq(1, length(row.names(geneMean)))
colnames(geneMean) <- c("gene", "MDSC", "NK", "Tumor")
geneMean <- as.data.frame(geneMean)
geneOfpairs.txt文件代表的是候補(bǔ)的配受體基因，讀入它隨后用于篩選exp表達(dá)矩陣刑峡。隨后我們計(jì)算NK洋闽、Tumor、MDSC各自的基因的表達(dá)水平突梦，并封裝成geneMean數(shù)據(jù)框形式
if(F) {
cal0 <- function(data){ # 計(jì)算表達(dá)相應(yīng)基因的比例
x <- data[data!=0]
p <- length(x) / length(data)
return(p)
} #計(jì)算向量中非零的比率
calGenePro <- function(data){
rowname <- rownames(data) #數(shù)據(jù)的行名
ret <- apply(data, 1, cal0)
names(ret) <- rowname
return(ret)
}
MDSC_Gene_Pro <- calGenePro(exp[MDSC_c]) #MDSC類型細(xì)胞中表達(dá)相應(yīng)基因的比例
NK_Gene_Pro <- calGenePro(exp[NK_c])
Tumor_Gene_Pro <- calGenePro(exp[Tumor_c])
geneMean$MDSC <- round(MDSC_Gene_Mean*(ifelse(MDSC_Gene_Pro >= 0.1, 1, 0)), 3)
geneMean$NK <- round(NK_Gene_Mean*(ifelse(NK_Gene_Pro >= 0.1, 1, 0)), 3)
geneMean$Tumor <- round(Tumor_Gene_Mean*(ifelse(Tumor_Gene_Pro >= 0.1, 1, 0)), 3)
}
write.table(geneMean, "./output/gene_means.txt", sep = "\t", quote = F, row.names = F)
接下來(lái)是計(jì)算基因在各個(gè)分群中表達(dá)的比例诫舅，因?yàn)镻DCD1在3個(gè)分群中表達(dá)比例都很少，小于10%（即基因在某個(gè)分群中表達(dá)比例小于10%宫患，那么此基因應(yīng)當(dāng)在該分群中被剔除）刊懈，而后面我們要計(jì)算CD274-PDCD1在我們數(shù)據(jù)中的interaction強(qiáng)弱并畫圖，就沒(méi)篩選了基因（如要運(yùn)行娃闲，將F改為T）虚汛，geneMean最終形式為：
# 計(jì)算各個(gè)配受體對(duì)的平均表達(dá)量
LRpairs <- read.table("L-Rpairs.txt", header = T)
noLR <- LRgene$Gene[!(LRgene$Gene %in% geneMean$gene)]
LRpairs <- LRpairs[!(LRpairs$Ligand %in% noLR),] # 去除用不到的LRpairs
LRpairs <- LRpairs[!(LRpairs$Receptor %in% noLR),]
typePairs <- c("Tumor-NK","Tumor-MDSC", "NK-Tumor", "NK-MDSC", "MDSC-NK", "MDSC-Tumor")
temp <- c(0, 0) # 臨時(shí)存儲(chǔ)配受體表達(dá)量
names(temp) <- paste("celltype", 1:2, sep = "")
means <- matrix(0, nrow = length(rownames(LRpairs)), ncol = 7,
dimnames = list(1:length(rownames(LRpairs)), c("LR", typePairs)))
for (p in 1:length(rownames(LRpairs))) { # 計(jì)算平均表達(dá)量
pair <- as.character(LRpairs[p,])
n <- 2
for (typePair in typePairs) {
typePair <- strsplit(typePair, split = "-")[[1]]
for (i in 1:length(pair)) {
index <- grep(pair[i], geneMean$gene, fixed = TRUE)
temp[i] <- geneMean[index, typePair[i]]
}
temp <- as.numeric(temp)
expLevel <- ifelse(!(temp[1] & temp[2]), 0, (temp[1]+temp[2])/2)
means[p, n] <- round(expLevel, 3)
n <- n + 1
}
means[p, 1] <- paste(pair[1], pair[2], sep = "-")
}
means <- as.data.frame(means)
write.table(means, "./output/means.txt", sep = "\t", quote = F, row.names = F)
這部分代碼即計(jì)算候補(bǔ)的配受體對(duì)的表達(dá)量，如果配受體任何一方表達(dá)量為0皇帮，那么最終的表達(dá)量為0卷哩，否則為兩者的均值。最后我們會(huì)得出means數(shù)據(jù)框属拾，其內(nèi)容及形式為：
# 計(jì)算各個(gè)配受體的p-values
p_values <- matrix(0, nrow = 129, ncol = 6,
dimnames = list(1:length(rownames(LRpairs)), typePairs))
for (i in 1:1000) {
ann_Random <- cbind(sample(ann$celltype, length(ann$celltype), replace = F), ann$cell) # shuffle
ann_Random <- as.data.frame(ann_Random)
colnames(ann_Random) <- c("celltype", "cell")
MDSC_c_Random <- ann_Random$cell[ann_Random$celltype=="MDSC"] #MDSC類型細(xì)胞
NK_c_Random <- ann_Random$cell[ann_Random$celltype=="NK"]
Tumor_c_Random <- ann_Random$cell[ann_Random$celltype=="Tumor"]
MDSC_Gene_Mean_Random <- round(rowMeans(exp[MDSC_c_Random]), 3) #計(jì)算MDSC類型細(xì)胞的基因表達(dá)水平
NK_Gene_Mean_Random <- round(rowMeans(exp[NK_c_Random]), 3)
Tumor_Gene_Mean_Random <- round(rowMeans(exp[Tumor_c_Random]), 3)
geneMean_Random <- cbind(row.names(exp), MDSC_Gene_Mean_Random, NK_Gene_Mean_Random, Tumor_Gene_Mean_Random)
rownames(geneMean_Random) <- seq(1, length(row.names(geneMean_Random)))
colnames(geneMean_Random) <- c("gene", "MDSC", "NK", "Tumor")
geneMean_Random <- as.data.frame(geneMean_Random)
means_Random <- matrix(0, nrow = length(rownames(LRpairs)), ncol = 7,
dimnames = list(1:length(rownames(LRpairs)), c("LR", typePairs)))
for (p in 1:length(rownames(LRpairs))) { # 計(jì)算平均表達(dá)量
pair <- as.character(LRpairs[p,])
n <- 2
for (typePair in typePairs) {
typePair <- strsplit(typePair, split = "-")[[1]]
for (i in 1:length(pair)) {
index <- grep(pair[i], geneMean_Random$gene, fixed = TRUE)
temp[i] <- geneMean_Random[index, typePair[i]]
}
temp <- as.numeric(temp)
expLevel <- ifelse(!(temp[1] & temp[2]), 0, (temp[1]+temp[2])/2)
means_Random[p, n] <- round(expLevel, 3)
n <- n + 1
}
means_Random[p, 1] <- paste(pair[1], pair[2], sep = "-")
}
means_Random <- as.data.frame(means_Random)
TEMP <- means_Random[2:7] > means[2:7]
p_values <- p_values + TEMP
}
p_values <- p_values/1000 # 計(jì)算P值
p_values <- as.data.frame(p_values)
p_values$`Tumor-NK`[means$`Tumor-NK`==0] <- 1.0
p_values$`Tumor-MDSC`[means$`Tumor-MDSC`==0] <- 1.0
p_values$`NK-Tumor`[means$`NK-Tumor`==0] <- 1.0
p_values$`NK-MDSC`[means$`NK-MDSC`==0] <- 1.0
p_values$`MDSC-NK`[means$`MDSC-NK`==0] <- 1.0
p_values$`MDSC-Tumor`[means$`MDSC-Tumor`==0] <- 1.0
p_values$LR <- means$LR
p_values <- p_values[c(7, 1, 2, 3, 4, 5, 6)]
write.table(p_values, "./output/p-values.txt", sep = "\t", quote = F, row.names = F)
這是細(xì)胞通訊分析的關(guān)鍵部分将谊。每次將細(xì)胞的標(biāo)簽打亂，隨后得到了means_Random渐白，即打亂后計(jì)算得到的配受體對(duì)的表達(dá)量
P值的定義：當(dāng)原假設(shè)為真時(shí)尊浓，比所得到的樣本觀察結(jié)果更極端的結(jié)果出現(xiàn)的概率。將means_Random和means（原結(jié)果）進(jìn)行比較纯衍，如上所示眠砾，得到一個(gè)判定矩陣TEMP，里面是1（大于）和0（不大于）。重復(fù)1000次褒颈，相加得到的p_values就是比樣本觀察結(jié)果更極端的結(jié)果出現(xiàn)的次數(shù)柒巫，除以1000得到概率，符合P值的定義
而第二個(gè)紅框的代碼谷丸，是將means中的0值的P值替換為1堡掏，因?yàn)?代表不表達(dá)，所以沒(méi)有必要考慮這些P值刨疼，故都設(shè)為1泉唁。p_values內(nèi)容及形式為：
# 制作significant_means.txt 閾值為0.05
significant_means <- means
significant_means$`Tumor-NK` <- as.numeric(p_values$`Tumor-NK` <=0.05) * as.numeric(significant_means$`Tumor-NK`)
significant_means$`Tumor-MDSC` <- as.numeric(p_values$`Tumor-MDSC` <=0.05) * as.numeric(significant_means$`Tumor-MDSC`)
significant_means$`NK-Tumor` <- as.numeric(p_values$`NK-Tumor` <=0.05) * as.numeric(significant_means$`NK-Tumor`)
significant_means$`NK-MDSC` <- as.numeric(p_values$`NK-MDSC` <=0.05) * as.numeric(significant_means$`NK-MDSC`)
significant_means$`MDSC-NK` <- as.numeric(p_values$`MDSC-NK` <=0.05) * as.numeric(significant_means$`MDSC-NK`)
significant_means$`MDSC-Tumor` <- as.numeric(p_values$`MDSC-Tumor` <=0.05) * as.numeric(significant_means$`MDSC-Tumor`)
write.table(significant_means, "./output/significant_means.txt", sep = "\t", quote = F, row.names = F)
這段代碼就是挑選顯著的配受體，閾值為0.05--默認(rèn)顯著性水平a為0.05
# 繪制circos弦圖
library(circlize)
data <- matrix(0, nrow = nrow(significant_means)*(ncol(significant_means)-1), ncol = 3,
dimnames = list(rep(significant_means$LR, each = ncol(significant_means)-1),
c("cell1", "cell2", "relation")))
t <- 1
for (i in 1:nrow(significant_means)) { # 為data填寫數(shù)據(jù)
n <- 2
for (typePair in typePairs) {
typePair <- strsplit(typePair, split = "-")[[1]]
data[t, 1] <- typePair[1]
data[t, 2] <- typePair[2]
data[t, 3] <- significant_means[i, n]
n <- n + 1
t <- t + 1
}
}
data <- as.data.frame(data)
data$relation <- as.numeric(data$relation)
grid.col = NULL # 開始繪制
grid.col[c("MDSC", "NK", "Tumor")] <- c("red", "grey", "grey")
visible <- rep(FALSE, nrow(data))
visible[grep("CD274.PDCD1", rownames(data), fixed = TRUE)] <- TRUE
pdf(file="circlize.pdf", width=8, height=8, pointsize=8)
chordDiagram(data, grid.col = grid.col, order = c("NK", "MDSC", "Tumor"),
directional = 1, link.visible = visible,
direction.type = c("diffHeight", "arrows"),
link.arr.type = "big.arrow", diffHeight = uh(0.05, "mm"),
annotationTrack = c("name", "grid"),
link.sort = TRUE, link.decreasing = FALSE)
title(main = "CD274-PDCD1", cex.main = 1.8)
dev.off()
最后一段代碼就是利用和弦圖表示出CD274-PDCD1在該數(shù)據(jù)中的interaction強(qiáng)弱（NK揩慕，Tumor亭畜，MDSC），如下所示（可以看出interaction并不強(qiáng)）：
我們仿照了CellPhoneDB完成了細(xì)胞通訊分析迎卤。值得一提的是拴鸵，關(guān)于P值需不需要矯正以及篩選基因規(guī)則應(yīng)該具體問(wèn)題具體分析，不可一味與別人一致蜗搔，因?yàn)閿?shù)據(jù)是不同的劲藐。另外，除了和弦圖外樟凄，還可用多種形式表達(dá)細(xì)胞間的關(guān)系聘芜，有興趣的小伙伴可以都去試一下
