cellchat:细胞通讯与其可视化利器part1：单个数据集细胞通讯

cellchat教程第一弹，源代码来自于http://www.cellchat.org/
不得不说，这个包的可视化做的真的很漂亮。狠狠喜欢了！

第一部分：数据准备创建对象

CellChat需要两个用户输入:一个是细胞的基因表达数据，另一个是用户分配的细胞标签(即基于标签的模式)或单细胞数据的低维表示(即无标签模式)。对于后者，CellChat根据低维空间或伪时间轨迹空间中的细胞-细胞距离构建共享邻居图，自动对细胞进行分组。

library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
# Here we load a scRNA-seq data matrix and its associated cell meta data
load(url("https://ndownloader.figshare.com/files/25950872")) # This is a combined data from two biological conditions: normal and diseases
data.input = data_humanSkin$data # normalized data matrix

基因表达矩阵信息

meta = data_humanSkin$meta # a dataframe with rownames containing cell mata data
cell.use = rownames(meta)[meta$condition == "LS"] # extract the cell names from disease data# Prepare input data for CelChat analysis
data.input = data.input[, cell.use]
meta = meta[cell.use, ]
# meta = data.frame(labels = meta$labels[cell.use], row.names = colnames(data.input)) # manually create a dataframe consisting of the cell labels
unique(meta$labels) # check the cell labels
#>  [1] Inflam. FIB  FBN1+ FIB    APOE+ FIB    COL11A1+ FIB cDC2
#>  [6] LC           Inflam. DC   cDC1         CD40LG+ TC   Inflam. TC
#> [11] TC           NKT
#> 12 Levels: APOE+ FIB FBN1+ FIB COL11A1+ FIB Inflam. FIB cDC1 cDC2 ... NKT

meta 信息数据框

用户可以从数据矩阵、Seurat或singlecellexperexperiment对象中创建一个新的CellChat对象。如果输入是Seurat或singlecellexperexperiment对象，对象中的元数据将被默认使用，USER必须提供group。通过定义单元组。例如,小组。Seurat对象中的默认单元格标识by = ident。注意:如果用户加载以前计算的CellChat对象(版本<0.5.0)，请通过updateCellChat更新对象。
这个数据有12种细胞类型，还挺多的。

cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")
#> Create a CellChat object from a data matrix
#> Set cell identities for the new CellChat object
#> The cell groups used for CellChat analysis are  APOE+ FIB FBN1+ FIB COL11A1+ FIB Inflam. FIB cDC1 cDC2 LC Inflam. DC TC Inflam. TC CD40LG+ TC NKT

如果在创建CellChat对象时没有添加cell mata信息，USERS还可以稍后使用addMeta添加它，并使用setIdent设置默认cell身份。

cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")
#> Create a CellChat object from a data matrix
#> Set cell identities for the new CellChat object
#> The cell groups used for CellChat analysis are  APOE+ FIB FBN1+ FIB COL11A1+ FIB Inflam. FIB cDC1 cDC2 LC Inflam. DC TC Inflam. TC CD40LG+ TC NKT
cellchat <- addMeta(cellchat, meta = meta)
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在小鼠体内包含2,021个已验证的分子相互作用，包括60%的分泌自分泌/旁分泌信号相互作用，21%的细胞外基质(ECM)-受体相互作用和19%的细胞-细胞接触相互作用。 CellChatDB在人体中包含1939个已验证的分子相互作用，包括61.8%的旁分泌/自分泌信号相互作用，21.7%的细胞外基质(ECM)-受体相互作用和16.5%的细胞-细胞接触相互作用。

用户可以通过添加自己策划的配体-受体对来更新CellChatDB。

CellChatDB <- CellChatDB.human # use CellChatDB.mouse if running on mouse data
pdf("./databasehumancategory.pdf",width = 10,height = 4)
showDatabaseCategory(CellChatDB)
dev.off()# Show the structure of the database
dplyr::glimpse(CellChatDB$interaction)
#> Rows: 1,939
#> Columns: 11
#> $ interaction_name   <chr> "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2", "TGF…
#> $ pathway_name       <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "T…
#> $ ligand             <chr> "TGFB1", "TGFB2", "TGFB3", "TGFB1", "TGFB1", "TGFB…
#> $ receptor           <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGF…
#> $ agonist            <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist", "T…
#> $ antagonist         <chr> "TGFb antagonist", "TGFb antagonist", "TGFb antago…
#> $ co_A_receptor      <chr> "", "", "", "", "", "", "", "", "", "", "", "", ""…
#> $ co_I_receptor      <chr> "TGFb inhibition receptor", "TGFb inhibition recep…
#> $ evidence           <chr> "KEGG: hsa04350", "KEGG: hsa04350", "KEGG: hsa0435…
#> $ annotation         <chr> "Secreted Signaling", "Secreted Signaling", "Secre…
#> $ interaction_name_2 <chr> "TGFB1 - (TGFBR1+TGFBR2)", "TGFB2 - (TGFBR1+TGFBR2…# use a subset of CellChatDB for cell-cell communication analysis
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling") # use Secreted Signaling
# use all CellChatDB for cell-cell communication analysis
# CellChatDB.use <- CellChatDB # simply use the default CellChatDB# set the used database in the object
cellchat@DB <- CellChatDB.use

为了推断细胞状态特异性通讯，我们在一个细胞群中识别过表达的配体或受体，然后在配体或受体过表达时识别过表达的配体-受体相互作用。cellchat也提供一个功能，以投射基因表达数据到蛋白质-蛋白质相互作用(PPI)网络。具体来说，一种扩散过程被用来平滑基因表达值，基于他们的邻居定义在一个高置信度实验验证的蛋白质-蛋白质网络。这种功能在分析单细胞数据时非常有用，因为这种投影可以减少信号基因的缺失效应，特别是对于可能的配体/受体亚基零表达。人们可能会担心这种扩散过程可能引入人工误差，然而，它只会引入非常弱的通信。用户也可以跳过这一步，设置raw。在computeCommunProb()函数中使用= TRUE。

# set the used database in the object
cellchat@DB <- CellChatDB.use# subset the expression data of signaling genes for saving computation cost
cellchat <- subsetData(cellchat) # This step is necessary even if using the whole database
future::plan("multiprocess", workers = 8) # do parallel
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
# project gene expression data onto PPI network (optional)
cellchat <- projectData(cellchat, PPI.human)

第二部分:细胞间通信网络的推理

CellChat通过为每个interacion赋值一个概率值并执行置换检验来推断生物学上重要的细胞-细胞通信。CellChat利用质量作用定律，将基因表达与信号配体、受体和它们的辅因子之间的相互作用的已知知识相结合，从而模拟细胞与细胞之间的通信概率。

推断出的配体-受体对的数目显然取决于计算每个细胞组平均基因表达的方法。默认情况下，CellChat使用一种统计上稳健的平均方法，称为trimean，它产生的interaction比其他方法少。然而，我们发现CellChat在预测更强的interaction方面表现得很好，这对于缩小interaction范围以进行进一步的实验验证非常有帮助。在computeCommunProb中，我们提供了使用其他方法(如5%和10%截断平均值)来计算平均基因表达的选择。值得注意的是，trimean近似于25%的截断平均值，意味着如果一组中表达细胞的百分比小于25%，则平均基因表达为零。要使用10%的截断平均值，USER可以设置type = "truncatedMean"和trim = 0.1。如果研究的生物过程中非常知名的信号通路没有被预测，USER可以使用不同的修剪值尝试truncatedMean。功能computeAveExpr可以帮助检查感兴趣的信号基因的平均表达，如computeAveExpr(cellchat, features =c(“CXCL12”，“CXCR4”)， type = “truncatedMean”， trim = 0.1)。

在分析未排序的单细胞转录组时，假设丰富的细胞群比稀少的细胞群更倾向于集体发送更强的信号，CellChat还可以在概率计算中考虑各细胞群中细胞比例的影响。用户可以设置population。size= TRUE。

cellchat <- computeCommunProb(cellchat)
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat <- filterCommunication(cellchat, min.cells = 10)
test=subsetCommunication(cellchat)

做完这一步，
cellchat对象里会出现一个slot ‘net’，object@net p r o b 是推断的通信概率矩阵，第一个、第二个、第三个维度分布代表 s o u r c e t a r g e t 和受体配体对。 o b j e c t @ n e t prob是推断的通信概率矩阵，第一个、第二个、第三个维度分布代表source target和受体配体对。 object@net prob是推断的通信概率矩阵，第一个、第二个、第三个维度分布代表sourcetarget和受体配体对。object@netpval是每个interaction的p值。
可以通过subsetCommunication(object)访问细胞通讯，返回一个数据框。

CellChat通过总结与每个信号通路相关的所有配体-受体相互作用的通信概率来计算信号通路水平上的通信概率。NB:推断出的每对配体受体和每条信号通路的细胞间通讯网络分别存储在net和netP中。

cellchat <- computeCommunProbPathway(cellchat)

提取函数的用法

df.net <- subsetCommunication(cellchat)
#returns a data frame consisting of all the inferred cell-cell communications at the level of ligands/receptors. Set slot.name = "netP" to access the the inferred communications at the level of signaling pathwaysdf.net <- subsetCommunication(cellchat, sources.use = c(1,2), targets.use = c(4,5))
#gives the inferred cell-cell communications sending from cell groups 1 and 2 to cell groups 4 and 5.df.net <- subsetCommunication(cellchat, signaling = c("WNT", "TGFb"))
#gives the inferred cell-cell communications mediated by signaling WNT and TGFb.

第一行和第二列是细胞类型，第三列是信号通路

配体受体水平

计算细胞-细胞聚合通信网络：我们可以通过计算连接数或汇总通信概率来计算胞-胞聚合通信网络。USER还可以通过设置源来计算单元组子集之间的聚合网络。使用souce.use和targets.use。

cellchat <- aggregateNet(cellchat)

Return an updated CellChat object:

‘object@net$count’ is a matrix: rows and columns are sources and targets respectively, and elements are the number of interactions between any two cell groups. USER can convert a matrix to a data frame using the function ‘reshape2::melt()’

count里面行列分别是信号发出、接受的细胞群，值是两个细胞群之间的interaction数目。
‘object@net$weight’ is also a matrix containing the interaction weights between any two cell groups

weights是权重。权重怎么算的？
我们也可以看到聚集的细胞之间的通信网络。例如，使用圆图显示任意两个细胞群之间的相互作用数量或总相互作用强度(权重)。

groupSize <- as.numeric(table(cellchat@idents))
pdf("./circle_plot.pdf",width = 8,height = 4)
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")
dev.off()

仔细看这个图，是有箭头的，颜色和source是一致的，圈圈的大小是每个细胞群细胞的数量。
由于细胞通信网络的复杂性，我们可以对每个细胞群发出的信号进行检测。这里我们还控制参数edge.weight.max这样我们就可以比较不同网络之间的边权值。

第三部分:细胞-细胞通信网络的可视化

在对细胞-细胞通信网络进行推断的基础上，CellChat为进一步的数据探索、分析和可视化提供了多种功能。

它提供了多种可视化细胞-细胞通信网络的方法，包括分层图、圆图、弦图和气泡图。
它提供了一个易于使用的工具来提取和可视化推断网络的高阶信息。例如，它允许对细胞群体的主要信号输入和输出进行预测，以及这些群体和信号如何协调发挥作用。
通过结合社会网络分析、模式识别和多种学习方法，它可以定量地描述和比较推断出的细胞-细胞通信网络。

cellchat@netP$pathways可以看所有pathways
层次图：
和弦图：
圈圈图：
热图：
计算每个配体-受体对整个信号通路的贡献，并将单个配体-受体对介导的细胞-细胞通信可视化

第四部分：细胞通讯系统分析

为了便于解释复杂的细胞间通信网络，CellChat通过从图论、模式识别和流形学习中抽象出来的方法对网络进行定量测量。

它可以利用网络分析中的中心性度量确定给定信令网络中的主要信令源和目标，以及中介和影响者
它可以预测特定细胞类型的关键输入和输出信号，并利用模式识别方法协调不同细胞类型之间的反应（这个预测有点迷）
它可以通过定义相似性度量来分组信号通路，并从功能和拓扑的角度进行manifold learning。
它可以通过多个网络的联合流形学习来描述保守的和context-specific的信号通路。

确定signaling角色(例如，主要的发送者，接收者)以及主要的贡献singnaling

# Compute the network centrality scores
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP") # the slot 'netP' means the inferred intercellular communication network of signaling pathways
# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
pdf("./cxcl.signalingrole.pdf",width = 8,height = 2.5)
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)
dev.off()

在2D空间中可视化主要的发送者(源)和接收者(目标)

# Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways
gg1 <- netAnalysis_signalingRole_scatter(cellchat)
#> Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways
# Signaling role analysis on the cell-cell communication networks of interest
gg2 <- netAnalysis_signalingRole_scatter(cellchat, signaling = c("CXCL", "CCL"))
#> Signaling role analysis on the cell-cell communication network from user's input
pdf("./cxcl.signalingrole2d.pdf",width = 8,height = 2.5)
gg1 + gg2
dev.off()

识别对某些细胞群的传出或传入信号贡献最大的信号

# Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways
ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming")
pdf("./outin.cell.pdf",width = 8,height = 6)
ht1 + ht2
dev.off()

确定全局细胞通讯模式，以探索如何多种细胞类型和信号通路协调在一起

除了探索单个通路的详细通信，一个重要的问题是多个细胞群和信号通路如何协调发挥作用。CellChat使用一种模式识别方法来识别全局通信模式。随着**模式(pattern)**数量的增加，可能会出现冗余模式，这使得解释通信模式变得困难。我们选择了5个默认模式。一般来说，当模式的数量大于2时，它在生物学上是有意义的。此外，我们还提供了一个函数selectK来推断模式的数量，该函数基于NMF R包中实现的两个指标，包括Cophenetic和Silhouette。这两个度量都是基于共识矩阵的层次聚类来度量特定数量模式的稳定性。对于一个模式数量的范围，一个合适的模式数量是在Cophenetic和Silhouette值开始突然下降的一个。