从基因表达解读细胞之间的相互作用和交流

2015/05/02来源：CHY Blog评论845

目前单细胞转录组技术的快速发展，关于单细胞转录组数据的应用也越来越宽泛，其中细胞相互作用的探索尤为重要。

从基因表达解读细胞之间的相互作用和交流

细胞互作存在多种形式

自分泌信号转导是指细胞内通讯，细胞分泌配体，这些配体用于通过同源受体诱导同一细胞上表达的那些分子的细胞应答。
旁分泌细胞间的通讯不需要细胞间的接触，而是取决于信号分子在分泌后从一个细胞扩散到另一个细胞。
近分泌，即依赖于接触的细胞间通讯依赖于间隙连接或膜纳米管等其他结构，使信号分子直接在细胞之间传递，而不会分泌到细胞外。
内分泌细胞间的通讯代表细胞间的通讯，信号分子被分泌并通过诸如血浆的细胞外液传播很长一段距离。

细胞互作分析步骤

步骤1：通过转录组学分析样品或细胞，以测量基因的表达；
步骤2：然后对生成的数据进行预处理以构建基因表达矩阵，其中包含跨不同样品或细胞的每个基因的转录水平；
步骤3：从其他来源生成或获得参与细胞间通讯的相互作用蛋白列表，通常包括分泌蛋白和膜结合蛋白(分别为配体和受体)之间的相互作用；
步骤4：在基因表达矩阵中仅保留与相互作用蛋白相关的基因。
步骤5：它们的表达水平用作输入，使用评分函数[函数f (L，R)]，其中L和R分别是配体和受体的表达值来计算每个配体-受体对的交流得分。可以使用聚合函数[函数g (Cell 1，Cell 2)]，其中Cell 1和Cell 2都是这些细胞或相应样本的所有通讯得分，可以汇总这些通信得分以计算各个样本或细胞之间的总体交互状态；
步骤6：最后，可以通过Circos图和网络可视化来表示交流和汇总分数，以方便对结果进行分析解释。

Plant cell-cell communication

植物由于存在细胞壁，成为了一道阻隔细胞交流的物理屏障；但是目前研究发现植物还是存在两种机制进行细胞交流：

小分子（如肽和植物激素）通过细胞壁和质膜受体的扩散
在胞质分裂的过程中，通过plasmodesmata保持细胞的连接(可以将营养素，激素，调节蛋白和RNA等分子从一个细胞转运到另一个细胞)，包括primary plasmodesmata、secondary plasmodesmata

能够在植物细胞之间移动的蛋白称为：non-cell-autonomous proteins (NCAPs)
CPC is transported from atrichoblasts, where it is expressed, to trichoblasts, and accumulates in their nuclei.

植物中也支持膜受体配体介导的细胞间信号转导；通过转录因子进行细胞间交流、small RNA介导的细胞信号传导、

细胞配受体通识

主要介绍动物相关细胞信号交流知识；细胞通常使用化学信号进行交流。这些化学信号是由发送细胞产生的蛋白质或其他分子，通常由细胞分泌并释放到细胞外空间。在那里，它们可以像漂流瓶一样漂浮到邻近的细胞。一个受体只能识别一个(或几个)特定的配体，一个配体只能与一个(或几个)目标受体结合。配体与受体结合会改变受体的形状或活性，使其能够传递信号或直接在细胞内部产生变化。
受体分为：胞内受体、细胞表面受体
细胞表面受体：配体门控离子通道，G蛋白偶联受体，受体酪氨酸激酶

CellChat

输入文件：1.细胞基因表达数据；2.细胞标签
输入是均一化的数据(Seurat@assay$RNA@data)；如果用户提供counts数据，可以用normalizeData函数来均一化。对于细胞的信息，需要一个带有rownames的数据格式作为CellChat的输入。

library(CellChat)
library(ggplot2)
library(ggalluvial)
library(svglite)
library(Seurat)
library(SeuratData)
options(stringsAsFactors = FALSE)

# 根据Seurat对象准备CellChat输入数据
data.input  <- pbmc3k.final@assays$RNA@data
identity = data.frame(group =pbmc3k.final$seurat_annotations   , row.names = names(pbmc3k.final$seurat_annotations)) # create a dataframe consisting of the cell labels
unique(identity$group) # check the cell labels

# 构建cellchat对象
cellchat <- createCellChat(data = data.input)
cellchat

# 可视化数据结构
library(mindr)
(out <- capture.output(str(cellchat)))
out2 <- paste(out, collapse="n")
mm(gsub("\.\.@","# ",gsub("\.\. ","#",out2)),type ="text")

# 添加metadata到cellchat对象
cellchat <- addMeta(cellchat, meta = identity, meta.name = "labels")
cellchat <- setIdent(cellchat, ident.use = "labels") # set "labels" as default cell identity
levels(cellchat@idents) # show factor levels of the cell labels
groupSize <- as.numeric(table(cellchat@idents)) # number of cells in each cell group

# 导入配受体库
# 查看数据库结构
CellChatDB <- CellChatDB.human 
(out3 <- capture.output(str(CellChatDB)))
out4 <- paste(out3, collapse="n")
mm(gsub("\$","# ",gsub("\.\. ","#",out4)),type ="text")

# 配受体库提取
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling") # use Secreted Signaling for cell-cell communication analysis
cellchat@DB <- CellChatDB.use # set the used database in the object

# 表达数据预处理
# 在一个细胞组中识别过表达的配体或受体
cellchat <- subsetData(cellchat) # subset the expression data of signaling genes for saving computation cost
future::plan("multiprocess", workers = 4) 
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
cellchat <- projectData(cellchat, PPI.human) 

# 为每个相互作用分配一个概率值并进行置换检验来推断生物意义上的细胞-细胞通信
# cellchat <- computeCommunProb(cellchat)  注意这个函数如果你可以用就用
mycomputeCommunProb <-edit(computeCommunProb)  
environment(mycomputeCommunProb) <- environment(computeCommunProb)
cellchat <- mycomputeCommunProb(cellchat)  

# 聚合通信网络
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

# 层次图或圈图可视化
# 实体圆和空心圆分别表示源和目标。圆的大小与每个细胞组的细胞数成比例。边缘颜色与信源一致。线越粗，信号越强。
cellchat@netP$pathways
levels(cellchat@idents) 
vertex.receiver = seq(1,4)
pathways.show <- "MIF"
mynetVisual_aggregate(cellchat, signaling = pathways.show,  vertex.receiver = vertex.receiver, vertex.size = groupSize)
mynetVisual_aggregate(cellchat, signaling = c("MIF"), layout = "circle", vertex.size = groupSize,pt.title=20,vertex.label.cex = 1.7)
netAnalysis_contribution(cellchat, signaling = pathways.show)

netVisual_signalingRole(cellchat, signaling = pathways.show, width = 12, height = 2.5, font.size = 10)

# 识别分泌细胞外向交流模式
nPatterns = 5 
# 同样在这里遇到了bug，难道说是我没有安装好吗，de了它。
# cellchat <- myidentifyCommunicationPatterns(cellchat, pattern = "outgoing", k = nPatterns)  
myidentifyCommunicationPatterns <- edit(identifyCommunicationPatterns)
environment(myidentifyCommunicationPatterns) <- environment(identifyCommunicationPatterns)
cellchat <- myidentifyCommunicationPatterns(cellchat, pattern = "outgoing", k = nPatterns)
# Visualize the communication pattern using river plot
netAnalysis_river(cellchat, pattern = "outgoing")
# Visualize the communication pattern using dot plot
netAnalysis_dot(cellchat, pattern = "outgoing")

# 识别目标细胞的传入(incoming)通信模式
netAnalysis_river(cellchat, pattern = "incoming")
netAnalysis_dot(cellchat, pattern = "incoming")

celltalker

cellphonedb

library(Seurat)
library(SeuratData)
pbmc3k
head(pbmc3k@meta.data)

# 根据Seurat对象准备cellphonedb文件
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[,'seurat_annotations', drop=F])  
meta_data <- as.matrix(meta_data)
meta_dat[is.na(meta_data)] = "Unkown" #  细胞类型中不能有NA
write.table(meta_data, 'cellphonedb_meta.txt', sep='t', quote=F, row.names=F)

# linux运行
cellphonedb method statistical_analysis  cellphonedb_meta.txt  cellphonedb_count.txt --counts-data=gene_name
cellphonedb plot dot_plot 
cellphonedb plot heatmap_plot cellphonedb_meta.txt   

# R中绘图
pbmc='D:\SingleCell\out\' ##  outs 文件放在这里了。

library(psych)
library(qgraph)
library(igraph)

netf<- "count_network.txt"
mynet <- read.delim(paste0(pbmc,"count_network.txt"), check.names = FALSE)
head(mynet)
net<- graph_from_data_frame(mynet)
plot(net)

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)
coords <- layout_in_circle(net, order =
                             order(membership(karate_groups)))  # 设置网络布局

E(net)$width  <- E(net)$count/10  # 边点权重（粗细）
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) 

net2 <- net  # 复制一份备用

for (i in 1: length(unique(mynet$SOURCE)) ){
  E(net)[map(unique(mynet$SOURCE),function(x) {
    get.edge.ids(net,vp = c(unique(mynet$SOURCE)[i],x))
  })%>% unlist()]$color <- allcolour[i]
}  # 这波操作谁有更好的解决方案？ 

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) 

plot(net, edge.arrow.size=.1, 
     edge.curved=0.2, # 只是调了这个参数
     vertex.color=allcolour,
     vertex.frame.color="#555555",
     vertex.label.color="black",
     layout = coords,
     vertex.label.cex=.7) 

# 贝壳图
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)[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) 
}
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, # 绘制边的权重
       vertex.color=allcolour,
       vertex.frame.color="#555555",
       vertex.label.color="black",
       layout = coords,
       vertex.label.cex=1
  ) 

}

# 点图
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))