参考链接：
什么是细胞通讯
不同细胞类型和组织中的细胞之间相互作用(CCI)，主要是通过各种分子包括离子、代谢物、整联蛋白、受体、连接蛋白、结构蛋白、配体和细胞外基质的分泌蛋白来实现的。一些分子如激素、生长因子、趋化因子、细胞因子和神经递质等配体介导细胞之间的信息交流(CCC)。 人们可以从基因表达中推断出CCI和CCC，尤其是在单细胞转录组测序之后，以RNA为代表的细胞信息之间的相互作用更能够对其进行解析
主要是流程就是根据数据库记录的参与细胞间通讯的相互作用蛋白列表，筛选表达矩阵后，看配体和受体表达量根据某个评分函数的交流得分，多个配体和受体聚合成为细胞与细胞之间的总体交互状态，最后通过网络可视化展现出来
细胞通讯的方式
- 自分泌
自分泌信号转导是指细胞内通讯，细胞分泌配体，这些配体用于通过同源受体诱导同一细胞上表达的那些分子的细胞应答旁分泌- 旁分泌
旁分泌细胞间的通讯不需要细胞间的接触，而是取决于信号分子在分泌后从一个细胞扩散到另一个细胞- 近分泌
即依赖于接触的细胞间通讯依赖于间隙连接或膜纳米管等其他结构，使信号分子直接在细胞之间传递，而不会分泌到细胞外- 内分泌
内分泌细胞间的通讯代表细胞间的通讯，信号分子被分泌并通过诸如血浆的细胞外液传播很长一段距离
分析原理
在给定表达矩阵和细胞注释之后，对于gene1-gene2这个互作关系，计算某个clusterA里面gene1的表达均值，计算另一个clusterB中gene2的表达均值，二者的均值为MEAN；在随机更换细胞的label之后，依据新的标签，计算“clusterA”里面gene1的表达均值，"clusterB"中gene2的表达均值，再求一个平均值mean，这样的过程重复多次，就可以得到一个mean的分布，即null distribution。MEAN在这个分布中所在的位置以及更极端的位置，构成的占比，就是p值（p值的定义）。所以CellPhoneDB推测两种细胞类型之间显著富集的受配体关系，本质上还是基于一个细胞类型里面的受体表达量，以及另一种细胞类型里面的配体表达量。此外，如果某种关系无处不在（在所有细胞类型之间都很明显），则找不出来
代码实战
rm(list = ls()) # 清除所有对象
library(Matrix)
# 将工作空间设置到此处,后面用到的文件都在此空间中（根据自己需求设置）
setwd("D:/personal_file/Bioinfomatics Study/data")
getwd()
# [1] "D:/personal_file/Bioinfomatics Study/data"
# 读入数据并制作表达矩阵
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) # 转化为数据框便于操作
print(object.size(exp), unit="GB")
rm(new_matrix_data,matrix_data)
gc()
读入数据并整理，最终得到的是exp表达矩阵，表达矩阵的行代表基因，列代表细胞样本。因为读入后就是小数点形式的，表示该表达矩阵已normalization了
# 筛选出配受体基因并整理exp
LRgene <- read.table("geneOfpairs.txt", col.names = "Gene")
exp <- exp[rownames(exp) %in% LRgene$Gene,]
# 计算NK、Tumor和MDSC类型细胞的基因表达水平
ann <- read.table("mdsc_tumor_nk_celltype.txt", head=T, row.names = 1)
MDSC_c <- ann$cell[ann$celltype=="MDSC"] #MDSC类型细胞
NK_c <- ann$cell[ann$celltype=="NK"]
Tumor_c <- ann$cell[ann$celltype=="Tumor"]
MDSC_Gene_Mean <- rowMeans(exp[MDSC_c]) #计算MDSC类型细胞的基因表达水平
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文件代表的是候补的配受体基因，读入它随后用于筛选exp表达矩阵。随后我们计算NK、Tumor、MDSC各自的基因的表达水平，并封装成geneMean数据框形式
if(F) {
cal0 <- function(data){ # 计算表达相应基因的比例
x <- data[data!=0]
p <- length(x) / length(data)
return(p)
} #计算向量中非零的比率
calGenePro <- function(data){
rowname <- rownames(data) #数据的行名
ret <- apply(data, 1, cal0)
names(ret) <- rowname
return(ret)
}
MDSC_Gene_Pro <- calGenePro(exp[MDSC_c]) #MDSC类型细胞中表达相应基因的比例
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)
接下来是计算基因在各个分群中表达的比例，因为PDCD1在3个分群中表达比例都很少，小于10%（即基因在某个分群中表达比例小于10%，那么此基因应当在该分群中被剔除），而后面我们要计算CD274-PDCD1在我们数据中的interaction强弱并画图，就没筛选了基因（如要运行，将F改为T），geneMean最终形式为：
# 计算各个配受体对的平均表达量
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) # 临时存储配受体表达量
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))) { # 计算平均表达量
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)
这部分代码即计算候补的配受体对的表达量，如果配受体任何一方表达量为0，那么最终的表达量为0，否则为两者的均值。最后我们会得出means数据框，其内容及形式为：
# 计算各个配受体的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类型细胞
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) #计算MDSC类型细胞的基因表达水平
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))) { # 计算平均表达量
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 # 计算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)
这是细胞通讯分析的关键部分。每次将细胞的标签打乱，随后得到了means_Random，即打乱后计算得到的配受体对的表达量
P值的定义：当原假设为真时，比所得到的样本观察结果更极端的结果出现的概率。将means_Random和means（原结果）进行比较，如上所示，得到一个判定矩阵TEMP，里面是1（大于）和0（不大于）。重复1000次，相加得到的p_values就是比样本观察结果更极端的结果出现的次数，除以1000得到概率，符合P值的定义
而第二个红框的代码，是将means中的0值的P值替换为1，因为0代表不表达，所以没有必要考虑这些P值，故都设为1。p_values内容及形式为：
# 制作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--默认显著性水平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填写数据
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在该数据中的interaction强弱（NK，Tumor，MDSC），如下所示（可以看出interaction并不强）：
我们仿照了CellPhoneDB完成了细胞通讯分析。值得一提的是，关于P值需不需要矫正以及筛选基因规则应该具体问题具体分析，不可一味与别人一致，因为数据是不同的。另外，除了和弦图外，还可用多种形式表达细胞间的关系，有兴趣的小伙伴可以都去试一下
