hello，大家好，今天我们来分享一个更加重要的通讯方法，cellcall，规则很简单，就是我们在对单细胞进行通讯分析的时候，通讯到底有没有起到作用？？？之前的我分享的文章10X单细胞通讯分析之scMLnet（配受体与TF，差异基因（靶基因）网络通讯分析）、10X单细胞（10X空间转录组）通讯分析之NicheNet，都提到过这个问题，今天，我们就要把配受体和TF因子联合起来分析，得到真正起到通讯作用的配受体，真正了解细胞之间的交流。文章在CellCall: integrating paired ligand–receptor and transcription factor activities for cell–cell communication,2021年8月2号发表于Nucleic Acids Research，影响因子16分。

Overview of CellCall

CellCall 是一个通过整合细胞内和细胞间信号来推断细胞间通讯网络和内部调节信号的工具。

特点

• 1、 CellCall collects ligand-receptor-transcript factor (L-R-TF) axis datasets based on KEGG pathways（最大的亮点）。

• 2、 According to prior knowledge of L-R-TF interactions, CellCall infers intercellular communication by combining the expression of ligands/receptors and downstream TF activities for certain L-R pairs. （就是要拿到真正起到作用的配受体对）。

• 3、CellCall embeds a pathway activity analysis method to identify the crucial pathways involved in communications between certain cell types. （这一点也很重要，上升到通路的高度）。

• 4、CellCall offers a rich suite of visualization options (Circos plot, Sankey plot, bubble plot, ridge plot, etc.) to intuitively present the analysis results.（可视化，最后的呈现方式）。

The overview figure of CellCall is shown as follows.

图片.png

ligand-receptor-transcript factor (L-R-TF)的分析模式，非常重要，比原来的方法进了一大步。

Main functions of CellCall

CellCall provides a variety of functions including intercellular communication analysis, pathway activity analysis and a rich suite of visualization tools to intuitively present the results of the analysis (including Heatmap, Circos plot, Bubble plot, Sankey plot, TF enrichment plot and Ridge plot).

首先看看Intercellular communication analysis

加载数据（说白了就是矩阵，seurat的分析结果rds完全匹配）

The format of the input file is as follow table:

• 1. The row names: gene symbols.
• 2. The column names: cell IDs. The colnames can't contain punctuation such as commas, periods, dashes, etc. Using underline to connect barcoder and cell type is recommended. Take the input format below as an example, the column name is made up of index and cell type. Users should set names.field=2 and names.delim="_" in the function CreateNichConObject(). After that, cell type information is obtained and stored in the S4 object for later analysis. Because method in this paper depends on the cell type information, obtaining celltype information correctly is important.
• 3. Other place: the expression values (counts or TPM) for a gene in a cell.

图片.png

注意这里的行名，写上了细胞类型。

Create object

 mt <- CreateNichConObject(data=data, min.feature = 3,
                            names.field = 2,
                            names.delim = "_",
                            source = "TPM",
                            scale.factor = 10^6,
                            Org = "Homo sapiens",
                            project = "Microenvironment")

Arguments
data	A dataframe with row of gene and column of sample and the value must be numeric. Meanwhile what matters is that the colnames of dataframe should be in line with the paramter 'names.delim' and 'names.field', the former for pattern to splite every colnames, the latter for setting which index in splited colnames is cell type information. The function can get the 'CELLTYPE' information from the colnames 'BARCODE_CLUSTER_CELLTYPE' with names.delim="_" and names.field='3', and then stored in slot meta.data of CreateNichCon. Cell type annotation from every cell is essential for scoring cell communication. If the colnames of data don't coincide with the paramter 'names.delim' and 'names.field', CreateNichCon object may fail to create.
min.feature	Include cells where enough features equalling min.feature are detected. It's a preprocess which is the same as Seurat and set min.feature=0, if you don't want to filter cell. This parameter depends on the sequencing technology of the input data.
names.delim	Set the pattern to splite column name into vector. If the column name of the input matrix is BARCODE_CLUSTER_CELLTYPE, set names.delim="_" to get CELLTYPE of BARCODE_CLUSTER_CELLTYPE with names.field=3.
names.field	Set the index of column name vector which is splited by parameter names.delim to get cell type information. If the column name of the input matrix is BARCODE_CLUSTER_CELLTYPE, set names.field=3 to get CELLTYPE of BARCODE_CLUSTER_CELLTYPE with names.delim="_".
source	The type of expression dataframe, eg "UMI", "fullLength", "TPM", or "CPM". When the source of input data is "TPM" or "CPM", no transformation on the data. Otherwise, we transform the data to TPM with the parameter source="fullLength" and to CPM with source="UMI".
scale.factor	Sets the scale factor for cell-level normalization, default "10^6", if the parameter is "UMI" or "fullLength". Otherwise this parameter doesn't work.
Org	Set the species source of gene, eg "Homo sapiens", "Mus musculus". This parameter decides the paired ligand-receptor dataset and the transcript length which is needed in "TPM" transformation.
project	Sets the project name for the CreateNichCon object.

Infer the cell-cell communication score

The communication score of an L-R interaction between cell types is evaluated by integrating the L2- norm of the L-R interaction and the activity score of the downstream TFs.

mt <- TransCommuProfile(object = mt,
                        pValueCor = 0.05,
                        CorValue = 0.1,
                        topTargetCor=1,
                        p.adjust = 0.05,
                        use.type="median",
                        probs = 0.9,
                        method="weighted",
                        IS_core = TRUE,
                        Org = 'Homo sapiens')

参数分析

Pathway activity analysis

CellCall embeds a pathway activity analysis method to help explore the main pathways involved in communication between certain cells.

n <- mt@data$expr_l_r_log2_scale

pathway.hyper.list <- lapply(colnames(n), function(i){
    print(i)
    tmp <- getHyperPathway(data = n, object = mt, cella_cellb = i, Org="Homo sapiens")
    return(tmp)
})

For pathway activity analysis, Bubble plot is adopted to present the analysis results. Function of getForBubble is used to merge the data and plotBubble is used to draw the bubble plot.

myPub.df <- getForBubble(pathway.hyper.list, cella_cellb=colnames(n))
p <- plotBubble(myPub.df)

图片.png

最后的可视化Visualization

Circle plot

  cell_color <- data.frame(color=c('#e31a1c','#1f78b4',
                                   '#e78ac3','#ff7f00'), stringsAsFactors = FALSE)
  rownames(cell_color) <- c("SSC", "SPGing", "SPGed", "ST")

Plotting circle with CellCall object:

ViewInterCircos(object = mt, font = 2, cellColor = cell_color, 
                lrColor = c("#F16B6F", "#84B1ED"),
                arr.type = "big.arrow",arr.length = 0.04,
                trackhight1 = 0.05, slot="expr_l_r_log2_scale",
                linkcolor.from.sender = TRUE,
                linkcolor = NULL, gap.degree = 2,
                order.vector=c('ST', "SSC", "SPGing", "SPGed"),
                trackhight2 = 0.032, track.margin2 = c(0.01,0.12), DIY = FALSE)

图片.png

Plotting circle with DIY dataframe of mt@data$expr_l_r_log2_scale:

  ViewInterCircos(object = mt@data$expr_l_r_log2_scale, font = 2, 
                  cellColor = cell_color,
                  lrColor = c("#F16B6F", "#84B1ED"),
                  arr.type = "big.arrow",arr.length = 0.04,
                  trackhight1 = 0.05, slot="expr_l_r_log2_scale",
                  linkcolor.from.sender = TRUE,
                  linkcolor = NULL, gap.degree = 2,
                  order.vector=c('ST', "SSC", "SPGing", "SPGed"),
                  trackhight2 = 0.032, track.margin2 = c(0.01,0.12), DIY = T)

图片.png

图片.png

Pheatmap plot

viewPheatmap(object = mt, slot="expr_l_r_log2_scale", show_rownames = T,
             show_colnames = T,treeheight_row=0, treeheight_col=10,
             cluster_rows = T,cluster_cols = F,fontsize = 12,angle_col = "45",  
             main="score")

图片.png

图片.png

Sankey plot

  mt <- LR2TF(object = mt, sender_cell="ST", recevier_cell="SSC",
              slot="expr_l_r_log2_scale", org="Homo sapiens")
  head(mt@reductions$sankey)
 if(!require(networkD3)){
        BiocManager::install("networkD3")
  }
  
  sank <- LRT.Dimplot(mt, fontSize = 8, nodeWidth = 30, height = NULL, width = 1200,         
                      sinksRight=FALSE, DIY.color = FALSE)
  networkD3::saveNetwork(sank, "~/ST-SSC_full.html")

图片.png

The first pillar is ligand，the second pillar is receptor and the last pillar is TF. And the color of left and right flow are consistent with ligand and TF respectively.

图片.png

The color depends on ligand and receptor (by function sankey_graph with isGrandSon = FALSE)

library(magrittr)
library(dplyr)
tmp <- mt@reductions$sankey
tmp1 <- dplyr::filter(tmp, weight1>0) ## filter triple relation with weight1 (LR score)
tmp.df <- trans2tripleScore(tmp1)  ## transform weight1 and weight2 to one value (weight)
head(tmp.df)

## set the color of node in sankey graph
mycol.vector = c('#5d62b5','#29c3be','#f2726f','#62b58f','#bc95df', '#67cdf2', '#ffc533', '#5d62b5', '#29c3be')  
elments.num <-  tmp.df %>% unlist %>% unique %>% length()
mycol.vector.list <- rep(mycol.vector, times=ceiling(elments.num/length(mycol.vector)))
sankey_graph(df = tmp.df, axes=1:3, mycol = mycol.vector.list[1:elments.num], nudge_x = NULL, font.size = 4, boder.col="white", isGrandSon = F)

The first pillar is ligand，the second pillar is receptor and the last pillar is TF. And the color of left and right flow are consistent with ligand and receptor respectively.

图片.png

The color depends on ligand (by function sankey_graph with isGrandSon = TRUE)

library(magrittr)
library(dplyr)
tmp <- mt@reductions$sankey
tmp1 <- dplyr::filter(tmp, weight1>0)  ## filter triple relation with weight1 (LR score)
tmp.df <- trans2tripleScore(tmp1)  ## transform weight1 and weight2 to one value (weight)

## set the color of node in sankey graph
mycol.vector = c('#9e0142','#d53e4f','#f46d43','#fdae61','#fee08b','#e6f598','#abdda4','#66c2a5','#3288bd','#5e4fa2')
elments.num <-  length(unique(tmp.df$Ligand))
mycol.vector.list <- rep(mycol.vector, times=ceiling(elments.num/length(mycol.vector)))

sankey_graph(df = tmp.df, axes=1:3, mycol = mycol.vector.list[1:elments.num], 
             isGrandSon = TRUE, nudge_x = nudge_x, font.size = 2, boder.col="white",            
             set_alpha = 0.8)

The first pillar is ligand，the second pillar is receptor and the last pillar is TF. And the color of left and right flow are consistent with ligand.

图片.png

图片.png

TF enrichment plot

TF enrichment plot is adopted to present the TF activities in receiver cells.

Obtain the gene sets (TGs of TF):

For some biologists, they pay more attention to the TGs of the TF. This tool provide options to present or export the details of the TG list, stored in the NichConObject@datacell_type@geneSets.

mt@data$gsea.list$SSC@geneSets

The figure presents a part of result stored in the NichConObject@datacell_type@geneSets.

图片.png

Show all TFs in the SSC:

  ssc.tf <- names(mt@data$gsea.list$SSC@geneSets)
  ssc.tf

图片.png

Draw the TF enrichment plot:

getGSEAplot(gsea.list=mt@data$gsea.list, geneSetID=c("CREBBP", "ESR1", "FOXO3"), 
            myCelltype="SSC", fc.list=mt@data$fc.list,  
            selectedGeneID = mt@data$gsea.list$SSC@geneSets$CREBBP[1:10],
            mycol = NULL)

图片.png

图片.png

Ridge plot

Ridge plot is adopted to present FC distribution of TGs of activated TFs.

## gsea object
egmt <- mt@data$gsea.list$SSC

## filter TF
egmt.df <- data.frame(egmt)
head(egmt.df[,1:6])
flag.index <- which(egmt.df$p.adjust < 0.05)

ridgeplot.DIY(x=egmt, fill="p.adjust", showCategory=flag.index, core_enrichment = T,
                orderBy = "NES", decreasing = FALSE)

图片.png

图片.png

生活很好，有你更好