10X单细胞（10X空间转录组）细胞通讯分析之InterCell

作者，追风少年i

hello，大家好，今天这个日子相信大家都不会陌生吧，不知道有没有当年复读的童鞋，日子很快，恐怕高考之后，再也没有经历过高中那种炼狱的日子。

今天我们分享一个好的内容，关于细胞通讯，其实关于细胞通讯，最终的目的还是要知道通讯的生物学作用，文章在InterCellDB: A User-Defined Database for Inferring Intercellular Networks,2022年发表在advanced science，单位是浙江大学，看看这个国产分析方法是否能媲美CellChat。

单细胞组学技术为分析细胞间网络提供了足够的分辨率。已经开发了几种分析工具来揭示细胞间相互作用，例如 CellPhoneDB、SingleCellSignalR、CellChat、iTALK 和 NicheNet。但是，该领域仍有未满足的需求。首先，所有当前的分析工具都提供有限类型的交互。实际上，细胞间网络要复杂得多，包括配体-受体、受体-受体、细胞外基质-受体和囊泡-细胞相互作用。此外，所有这些已发表的工具在提供特定功能的细胞间相互作用方面的能力有限。例如，在最近的研究中，发现没有可用的分析方法能够预测调节性 T 细胞（Treg 细胞）和小胶质细胞之间的配体-受体相互作用如何促进缺血性中风后的少突胶质细胞生成和白质修复（说白了就是通讯具体的生物学功能的挖掘）。最后，大多数当前方法都有固定的内置搜索标准，这限制了分析的维度。需要一个新平台来以用户定义的方式分析通讯网络。因此，开发了一个名为 InterCellDB 的新分析平台来解决这些差距，并使用单细胞 RNA 测序 (scRNA-seq) 数据集结合感兴趣的特定生物学功能，定义的细胞间通讯综合分析。

Overview of InterCellDB database and workflow.

首先看看软件特点

• 对配受体基因进行了位置和功能分类
• 差异基因纳入通讯分析，All the DEGs were annotated by matching to the gene database, including their cellular localization, functional classification, and their attribution in GO terms。当然，基因可以自我设定（比如某个GO通路的基因）。
• 置换检验
• 可视化做的不错

参考示例

install.packages("devtools")
install.packages("BiocManager")
BiocManager::install(pkgs = "GO.db")
devtools::install_github("ZJUDBlab/InterCellDB")
library(InterCellDB)  # 0.9.2
library(Seurat)       # 3.2.0
library(dplyr)        # 1.0.5
library(future)       # 1.21.0

1、Data preparation and create InterCell object

InterCellDB requires 2 data as input:

• normalized count matrix
• differentially expressed genes (DEGs) with their belonging clusters

seurat.obj <- readRDS(url("https://figshare.com/ndownloader/files/31497371"))
tmp.counts <- seurat.obj[["RNA"]]@data
colnames(tmp.counts) <- as.character(Idents(seurat.obj))
# show data structure, row names are genes, column names are clusters
tmp.counts[1:5, 1:5]
# 5 x 5 sparse Matrix of class "dgCMatrix"
#            cDC2        NK        NK       cDC2   Myeloid
# Gnai3 .         1.4526836 0.7893417 0.02189368 .        
# Cdc45 0.1342404 0.1345867 .         1.01913912 0.1698247
# Scml2 .         .         .         .          .        
# Apoh  .         .         .         .          .        
# Narf  .         .         .         0.02189368 .

make sure the data normalization and cell clustering are done.Take Seurat for example, the NormalizeData or SCTransform should be processed for getting normalized data, and FindClusters should be processed for getting cell clustering result

差异基因示例

markers.used <- readRDS(url("https://figshare.com/ndownloader/files/31497401"))
#                p_val    LogFC pct.1 pct.2          PVal cluster    gene
# Plbd1   2.436468e-234 1.751750 0.990 0.293 4.703844e-230    cDC2   Plbd1
# Cd209a  1.126670e-219 2.934390 0.642 0.099 2.175148e-215    cDC2  Cd209a
# Ms4a6c  2.003678e-217 1.306534 0.916 0.244 3.868301e-213    cDC2  Ms4a6c
# Alox5ap 2.043584e-213 1.118439 0.922 0.214 3.945343e-209    cDC2 Alox5ap
# Fgr     3.582238e-205 1.177831 0.775 0.154 6.915869e-201    cDC2     Fgr

To note, InterCellDB defines several mandatory column names when read in DEGs, which are columns named 'gene', 'cluster', 'LogFC', 'PVal', and their meanings are listed below:

• Column 'gene', the gene name.
• Column 'cluster', the belonging cluster of every DEG.
• Column 'LogFC', the log fold changes of each DEG.
• Column 'PVal', the statistical significance for the gene being defined as DEG.

Then, the InterCell object is created and we name it here as inter.obj. This variable is used for doing all the following analysis.

inter.obj <- CreateInterCellObject(markers.used, species = "mouse", add.exprs = TRUE, exprs.data = tmp.counts)

2、Customize interactions and proteins from databases

For customizing interactions, InterCellDB provides 4 options: evidence source, credibility score, action mode and action effect. To run this example data, we use all experimentally validated interactions, pathway-curated interactions and predicted interactions, and select those physically associated ones.

inter.obj <- SelectDBSubset(inter.obj,
        use.exp = TRUE,   # to use experimentally validated interactions
        exp.score.range = c(1, 1000),  # use all credibility score for 'use.exp'
        use.know = TRUE,  # to use pathway curated interactions
        know.score.range = c(1, 1000),  # use all credibility score for 'use.know'
        use.pred = TRUE,  # to use predicted interactions
        pred.score.range = c(1, 1000),  # use all credibility score for 'use.pred' 
        sel.action.mode = "binding",  # select physically associated ones
        sel.action.effect = "ALL")    # not limiting action effects

Details about 4 options on customizing interaction database are given in the supporting information of the article. Here, we list some of the most commonly used options.

• evidence sources: experimentally validated (use.exp), pathway curated (use.know), predicted (use.pred).
• credibility score: ranging from 1 to 1000. Highest confidence should >= 900, and high confidence should >= 700, and medium >= 400, while low < 400.
• action mode: giving the relation of interacting proteins. Commonly used one is 'binding', which selects physically associated protein interactions.
• action effect: giving the effect of protein interaction. Commonly used ones are 'positive' and 'negative'. 'positive' means one protein promotes the expression of the other, while 'negative' means the inhibition of the other.

For customizing included proteins, InterCellDB also provides 4 options: protein expression change, protein subcellular location, function and relevant biological process.

# fetch genes of interest. GO.db version is v3.8.2 here, different version may have non-identical result.
genes.receiver <- FetchGeneOI(inter.obj, 
        sel.location = "Plasma Membrane",  # fetch proteins located in plasma membrane
        sel.location.score = c(4:5),       # 4 and 5 are of high confidence
        sel.type = "Receptor",             # fetch those receptors
        sel.go.terms = "GO:0006955"        # fetch [GO:0006955]immune response related proteins
        )
# [1] "Fetch 160 genes of interest."
genes.sender <- FetchGeneOI(inter.obj, 
        sel.location = "Extracellular Region",  # fetch proteins located in extracellular region
        sel.location.score = c(4:5),            # 4 and 5 are of high confidence
        sel.go.terms = "GO:0006955"             # fetch [GO:0006955]immune response related proteins
        )
# [1] "Fetch 281 genes of interest."

To note, as the protein expression change is highly associated with the input data, InterCellDB put it to be done in later steps, see parameter sel.exprs.change in AnalyzeInterInFullView. If users want to explore more selections on the other 3 options, InterCellDB provides several functions to list all items for them.

• protein subcellular location: ListAllGeneLocation
• protein function: ListAllGeneType and ListAllGeneMergeType
• protein involved biological process: ListAllGeneGOTerm

3、Perform full network analysis

Full network analysis summarizes the interactions between every 2 cell clusters, and infers main participants by aggregated power and total count of gene pairs. The power of every gene pair is calculated by the product of expressions of 2 participating genes. The selected gene pairs are those with p-value < 0.05 in cell label permutation test (the p-value cutoff can also be set by users).

Full network analysis goes with 2 steps:

• generate cell label permutation list
• network analysis and filter statistically significant gene pairs

plan("multiprocess", workers = 2)  # package future provides the parallel interface, here create 2 parallel processes  
# using options(future.globals.maxSize=1024^2 * 1000) to get more space to run parallel process
tmp.permlist <- Tool.GenPermutation(inter.obj, tmp.counts, perm.times = 100)
plan("sequential")  # close the parallel processes

推断网络

plan("multiprocess", workers = 4)
inter.obj <- AnalyzeInterInFullView(inter.obj, 
        sel.some.genes.X = genes.receiver,  # set the genes for receiver cells
        sel.some.genes.Y = genes.sender,    # set the genes for sender cells
        sel.exprs.change = "Xup.Yup",       # select genes up-regulated for both sender cells and receiver cells
        run.permutation = TRUE,             # run statistical test with permutation list
        perm.expression = tmp.permlist,     # given the permutation list, generated by `Tool.GenPermutation`
        perm.pval.sides = "one-side-large", # use one-tailed test
        perm.pval.cutoff = 0.05)            # p-value cutoff 0.05
plan("sequential")

The results of full network analysis can be further collected by GetResultFullView for visualization or table output. Users can define their clusters of interest. Here, we select the interactions among clusters in tumor from example data, which exclude interactions involving 'LN T cell', 'Endo lymphatic', 'Endo LN' and 'Fibroblast LN' clusters.

torm.LN.clusters <- c("LN T cell", "Endo lymphatic", "Endo LN", "Fibroblast LN")
all.clusters <- ListAllClusters(inter.obj)                # list all clusters
used.clusters <- setdiff(all.clusters, torm.LN.clusters)  # select those clusters in tumor
used.clusters <- used.clusters[order(used.clusters)]
# show the result of full network analysis
fullview.result <- GetResultFullView(inter.obj, 
        show.clusters.in.x = used.clusters,
        show.clusters.in.y = used.clusters,
        plot.axis.x.name = "signal receiving cells",
        plot.axis.y.name = "signal sending cells",
        nodes.size.range = c(1,8))
Tool.ShowGraph(fullview.result)  # visualization
Tool.WriteTables(fullview.result, dir.path = "./")  # write result into csv file

图片.png

4、Perform intercellular analysis

Intercellular analysis focuses on one cell-cell interaction, and explores the results in 4 aspects:

• summarize action mode and action effect
• subset and rank gene pairs
• evaluate the specificity of gene pairs
• summarize gene pairs in spatial pattern

4.1 Summarize action mode and action effect

inter.obj <- FetchInterOI(inter.obj, cluster.x = "Myeloid", cluster.y = "CAF 1")
inter.obj <- AnalyzeInterInAction(inter.obj)
used.action.mode <- c("activation", "inhibition", "catalysis", "reaction", "expression", "ptmod")  # action mode: not showing mode 'binding' and 'other'
used.color.mode <- c("#FB8072", "#80B1D3", "#8DD3C7", "#FFFFB3", "#BEBADA", "#FDB462")  # customized color, one-to-one corresponding to those in `used.action.mode`
action.mode.result <- GetResultPieActionMode(inter.obj, 
        limits.exprs.change = c("Xup.Yup"),
        limits.action.mode = used.action.mode,
        color.action.mode = used.color.mode)
Tool.ShowGraph(action.mode.result)

图片.png

action.effect.result <- GetResultPieActionEffect(inter.obj, limits.exprs.change = c("Xup.Yup"))
Tool.ShowGraph(action.effect.result)

图片.png

4.2 Subset and rank gene pairs

Circumstances exist when too many gene pairs are returned or the initial selection of gene subsets are not well satisfied to users' needs. We provide additional function SelectInterSubset for picking up subsets of gene pairs in intercellular scale. All subset selection given in full network analysis are applicable in this function. In addition, action mode and action effect can be customized.

Here, we filter activated gene pairs, i.e. selecting those with either 'activation' action mode or 'positive' action effect.

inter.obj <- SelectInterSubset(inter.obj, 
        sel.action.mode = "activation",
        sel.action.effect = "positive",
        sel.action.merge.option = "union"
        )

result.inter.pairs <- GetResultTgCrosstalk(inter.obj, 
        func.to.use = list(
            Power = function(x) { log1p(x) },        # further log-transform the power
            Confidence = function(x) { 1/(1+x) }     # invert the confidence to show
        ),  
        axis.order.xy = c("Power", "Power"),       # order genes by power
        axis.order.xy.decreasing = c(TRUE, TRUE),  # select if ordering is decreasing
        nodes.size.range = c(1, 8))
Tool.ShowGraph(result.inter.pairs)

图片.png

We rank gene pairs by their power, and select those top ranked genes

inter.obj <- SelectInterSubset(inter.obj, 
        sel.some.genes.X = c("Itgb2","Itgam","C5ar1","Ccr2","C3ar1","Ccr1","Ccr5"),
        sel.some.genes.Y = c("C3","Ccl11","Cxcl12","Cxcl1","Ccl2","Ccl7","C4b")
        )

4.3 Evaluate the specificity of gene pairs

Besides the power rank, we also evaluate the specificity of the gene pairs

tmp.target.cluster.groups <- ListClusterGroups(inter.obj, 
        use.former = TRUE,
        cluster.former = c("Myeloid")  # list all crosstalk with Myeloid as receiver cell, this order is aligned with the X and Y axis in fullview result. In this example, clusters listed in X axis are the receiver cells.
        )
tmp.target.cluster.groups <- setdiff(tmp.target.cluster.groups, c("Myeloid~Myeloid", "Myeloid~Endo LN", "Myeloid~Endo lymphatic", "Myeloid~Fibroblast LN", "Myeloid~LN T cell"))

The analysis is processed in AnalyzeInterSpecificity and users can visualize the result by GetResultTgSpecificity.

inter.obj <- AnalyzeInterSpecificity(inter.obj, 
        to.cmp.cluster.groups = tmp.target.cluster.groups
        )

result.specificity <- GetResultTgSpecificity(inter.obj,
        sel.uq.cnt.options = seq_along(tmp.target.cluster.groups),
        plot.uq.cnt.merged = TRUE, 
        plot.name.X.to.Y = FALSE,
        func.to.use = list(Power = function(x) { log1p(x) },
                       Confidence = function(x) { 1/(1+x) }),
        dot.size.range = c(2,8),
        dot.colour.seq = c("#00809D", "#EEEEEE", "#C30000"),
        dot.colour.value.seq = c(0, 0.4, 1)
        )
Tool.ShowGraph(result.specificity)

图片.png

4.4 Summarize gene pairs in spatial pattern

蛋白质-蛋白质相互作用受到它们空间接近性的限制。 InterCellDB 为基因提供了相应的基因产物亚细胞位置，还提供了可视化方法以空间模式显示基因对。 可以将作用方式、作用效果、蛋白质亚细胞定位、基因表达变化整合到现实情况中。
Here, we choose only the plasma membrane protein and extracellular region protein to show the paracrine protein-protein interactions.

tmp.hide.locations <- setdiff(ListAllGeneLocation(inter.obj), c("Plasma Membrane", "Extracellular Region"))  # select the subcellular locations allowed to show
set.seed(101L)  # set seed to make reproducible result, or gene will be re-arranged every time re-running GetResultTgCellPlot
result.sptialpattern <- GetResultTgCellPlot(inter.obj, 
        area.extend.times = 20,  # control the size of plotting. Increase the value by 10 when warnings ask to
        hide.other.area = TRUE,
        hide.locations.X = tmp.hide.locations,
        hide.locations.Y = tmp.hide.locations,
    expand.gap.radius.list = list(ECM = 8, CTP = 2, NC = 2, OTHER = 2),
        link.size = 0.3,
        link.alpha = 0.8,
        legend.manual.left.spacing = grid::unit(0.1, "cm")
        )
Tool.ShowGraph(result.sptialpattern)

图片.png

生活很好，有你更好，百度文库出现了大量抄袭我的文章，对此我深表无奈，我写的文章，别人挂上去赚钱，抄袭可耻，挂到百度文库的人更可耻