library(CellChat)
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: igraph
## Warning: package 'igraph' was built under R version 4.3.2
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:dplyr':
##
## as_data_frame, groups, union
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.2
library(dplyr)
library(Seurat)
## Warning: package 'Seurat' was built under R version 4.3.2
## Loading required package: SeuratObject
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.3.2
## 'SeuratObject' was built with package 'Matrix' 1.6.3 but the current
## version is 1.6.5; it is recomended that you reinstall 'SeuratObject' as
## the ABI for 'Matrix' may have changed
##
## Attaching package: 'SeuratObject'
## The following object is masked from 'package:BiocGenerics':
##
## intersect
## The following object is masked from 'package:base':
##
## intersect
##
## Attaching package: 'Seurat'
## The following object is masked from 'package:igraph':
##
## components
library(patchwork)
library(ComplexHeatmap)
## Warning: package 'ComplexHeatmap' was built under R version 4.3.1
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.18.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(RColorBrewer)
library(ggplot2)
library(compiler)
library(presto)
## Loading required package: Rcpp
## Loading required package: data.table
## Warning: package 'data.table' was built under R version 4.3.2
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
#Import seurat object for cases
Pap1 <- readRDS("~/Desktop/Mureen.dec6/Pap1.rds")
Pap1
## An object of class Seurat
## 54918 features across 27446 samples within 3 assays
## Active assay: integrated (2000 features, 2000 variable features)
## 2 layers present: data, scale.data
## 2 other assays present: RNA, SCT
## 3 dimensional reductions calculated: pca, umap, harmony
Control <- readRDS("~/Desktop/Mureen.dec6/Control.rds")
Control
## An object of class Seurat
## 54918 features across 27392 samples within 3 assays
## Active assay: integrated (2000 features, 2000 variable features)
## 2 layers present: data, scale.data
## 2 other assays present: RNA, SCT
## 3 dimensional reductions calculated: pca, umap, harmony
#CellChat necessitates two types of user input: firstly, the gene expression data of the cell, and secondly, either a user-defined cell label (referred to as label-based mode) or a reduced-dimensional representation of the single-cell data (known as label-free mode). In the latter scenario, CellChat automatically clusters cells by generating a shared neighbor graph based on their proximity in reduced-dimensional space or pseudo-temporal trajectory space.
#Data Download
#For a gene expression data matrix, genes should be organized in rows with associated row names, while cells should be arranged in columns with corresponding column names. Input for CellChat analysis requires normalized data, such as library size normalization followed by a logarithmic transformation with a pseudo-count of 1. In the event of count data provided by the user, we offer a normalizeData function to adjust for library size before conducting the logarithmic transformation. For information regarding cell grouping, a data frame with row names is necessary as input for CellChat.
#The user can create a new CellChat object from a data matrix, Seurat or SingleCellExperiment object. If the input is a Seurat or SingleCellExperiment object, the metadata in the object is used by default, and the user must provide group.by to define the cell group. For example, for the default cell identity in a Seurat object, group.by = "ident".
#Note: If the user loads a previously calculated CellChat object (version < 0.5.0), please update the object via updateCellChat.
meta <- Pap1@meta.data
labels <- as.factor(Pap1@meta.data$labels)
cellchat <- createCellChat(Pap1$RNA@data, meta = Pap1@meta.data, group.by = "labels")
## [1] "Create a CellChat object from a data matrix"
## Warning in createCellChat(Pap1$RNA@data, meta = Pap1@meta.data, group.by =
## "labels"): The 'meta' data does not have a column named `samples`. We now add
## this column and all cells are assumed to belong to `sample1`!
## Set cell identities for the new CellChat object
## The cell groups used for CellChat analysis are Astrocyte, Endothelial Cell, GABA-ergic Neuron, Glutamatergic Neuron, Microglia, Oligodendrocyte, Oligodendrocyte Precursor Cell, Unclassified
Pap1.cellchat <- cellchat
Pap1.cellchat
## An object of class CellChat created from a single dataset
## 28608 genes.
## 27446 cells.
## CellChat analysis of single cell RNA-seq data!
#If the cell meta information is not added when creating the CellChat object, the user can also add it later using addMeta and set the default cell identity using setIdent.
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
## [1] "Astrocyte" "Endothelial Cell"
## [3] "GABA-ergic Neuron" "Glutamatergic Neuron"
## [5] "Microglia" "Oligodendrocyte"
## [7] "Oligodendrocyte Precursor Cell" "Unclassified"
groupSize <- as.numeric(table(cellchat@idents)) # number of cells in each cell group
groupSize
## [1] 734 17 4550 16745 836 3581 707 276
#Before users can use CellChat to infer cell-to-cell communication, they need to set up a ligand-receptor interaction database and identify overexpressed ligands or receptors.
#CellChatDB is a hand-curated database of human and mouse literature supporting ligand-receptor interactions. CellChatDB v2 contains ~3,300 validated molecular interactions, including ~40% secretory autocrine/paracrine signaling interactions, ~17% extracellular matrix (ECM)-receptor interactions, ~13% cell- Cell contact interactions and approximately 30% of non-protein signaling. Compared with CellChatDB v1, CellChatDB v2 adds more than 1000 protein and non-protein interactions, such as metabolism and synaptic signaling. It should be noted that for gene-related molecules that are not directly measured in single-cell RNA sequencing, CellChat v2 estimates the expression of ligands and receptors, using key mediators or enzymes of these molecules, in order to potentially pass through non-protein mediated communication.
#Users can update CellChatDB by adding their own curated ligand-receptor pairs.
CellChatDB <- CellChatDB.mouse # use CellChatDB.human if running on mouse data
showDatabaseCategory(CellChatDB)
# Show the structure of the database
dplyr::glimpse(CellChatDB$interaction)
## Rows: 3,379
## Columns: 28
## $ interaction_name <chr> "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2",…
## $ pathway_name <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb…
## $ ligand <chr> "Tgfb1", "Tgfb2", "Tgfb3", "Tgfb1", "Tgfb1", …
## $ receptor <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1…
## $ agonist <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist…
## $ antagonist <chr> "TGFb antagonist", "TGFb antagonist", "TGFb a…
## $ co_A_receptor <chr> "", "", "", "", "", "", "", "", "", "", "", "…
## $ co_I_receptor <chr> "TGFb inhibition receptor", "TGFb inhibition …
## $ evidence <chr> "KEGG: mmu04350", "KEGG: mmu04350", "KEGG: mm…
## $ annotation <chr> "Secreted Signaling", "Secreted Signaling", "…
## $ interaction_name_2 <chr> "Tgfb1 - (Tgfbr1+Tgfbr2)", "Tgfb2 - (Tgfbr1+…
## $ is_neurotransmitter <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL…
## $ ligand.symbol <chr> "Tgfb1", "Tgfb2", "Tgfb3", "Tgfb1", "Tgfb1", …
## $ ligand.family <chr> "TGF-beta", "TGF-beta", "TGF-beta", "TGF-beta…
## $ ligand.location <chr> "Extracellular matrix, Secreted, Extracellula…
## $ ligand.keyword <chr> "Disease variant, Signal, Reference proteome,…
## $ ligand.secreted_type <chr> "growth factor", "growth factor", "cytokine;g…
## $ ligand.transmembrane <lgl> FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALS…
## $ receptor.symbol <chr> "Tgfbr1, Tgfbr2", "Tgfbr1, Tgfbr2", "Tgfbr1, …
## $ receptor.family <chr> "Protein kinase superfamily, TKL Ser/Thr prot…
## $ receptor.location <chr> "Cell membrane, Secreted, Membrane raft, Cell…
## $ receptor.keyword <chr> "Membrane, Secreted, Disulfide bond, Kinase, …
## $ receptor.surfaceome_main <chr> "Receptors", "Receptors", "Receptors", "Recep…
## $ receptor.surfaceome_sub <chr> "Act.TGFB;Kinase", "Act.TGFB;Kinase", "Act.TG…
## $ receptor.adhesome <chr> "", "", "", "", "", "", "", "", "", "", "", "…
## $ receptor.secreted_type <chr> "", "", "", "", "", "", "", "", "", "", "", "…
## $ receptor.transmembrane <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRU…
## $ version <chr> "CellChatDB v1", "CellChatDB v1", "CellChatDB…
# use a subset of CellChatDB for cell-cell communication analysis
#CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling", key = "annotation") # 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
## An object of class CellChat created from a single dataset
## 28608 genes.
## 27446 cells.
## CellChat analysis of single cell RNA-seq data!
# subset the expression data of signaling genes for saving computation cost
cellchat <- subsetData(cellchat) # This step is necessary even if using the whole database
#futurhttps://rpubs.com/marswh/boxplote::plan("multisession", workers = 4) # do parallel
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
## The number of highly variable ligand-receptor pairs used for signaling inference is 2520
#execution.time = Sys.time() - ptm
#print(as.numeric(execution.time, units = "secs"))
#project gene expression data onto PPI (Optional: when running it, USER should set `raw.use = FALSE` in the function `computeCommunProb()` in order to use the projected data)
#To infer communication for a specific cell state, CellChat identifies overexpressed ligands or receptors in a population of cells and then identifies overexpressed ligand-receptor interactions if the ligand or receptor is overexpressed.
#Author of CellChat also provide a function to project gene expression data onto a protein-protein interaction (PPI) network. Specifically, a diffusion process is used to smooth the expression values of genes based on neighbors defined in a high-confidence experimentally validated protein-protein network. This feature is useful when analyzing single-cell data at shallow sequencing depth, as projection reduces dropout effects on signaling genes, especially for ligand/receptor subunits that may have zero expression. There may be concerns about possible human error introduced by this diffusion process, however, it will only introduce very weak communication. By default, CellChat uses raw data (i.e. object@data.signaling) instead of projected data. To use projected data, users should run the projectData function before running computeCommunProb and then set raw.use = FALSE when running computeCommunProb.
# cellchat <- projectData(cellchat, PPI.human)
#CellChat infers biologically important intercellular communication by assigning a probability value to each interaction and performing permutation tests. CellChat uses the law of mass action to model the probability of cell-to-cell communication by integrating known interactions between gene expression and signaling ligands, receptors, and their cofactors.
###
#The number of inferred ligand-receptor pairs obviously depends on the method used to calculate average gene expression per cell group. By default, CellChat uses a statistically robust averaging method called "trimean", which produces fewer interactions than other methods. However, we found that CellChat performed well in predicting stronger interactions, which was very helpful in further narrowing down the experimentally validated interactions. In computeCommunProb, we provide an option to use other methods to calculate average gene expression, such as 5% and 10% truncated averages. It is worth noting that "trimean" approximates a 25% truncated mean, meaning that if the percentage of expressing cells in a group is less than 25%, the average gene expression is zero. To use a 10% truncated mean, the user can set type = “truncatedMean” and trim = 0.1. To determine the appropriate trim value, CellChat provides a function computeAveExpr, which can help examine the average expression of a signaling gene of interest, e.g.
#cellchat <- computeCommunProb(cellchat, type = "triMean") #(run this when you are compiling this pipeline, I done this and save as rds file so for this tutorial I load this as rds file)
cellchat <- readRDS("/Users/usri/Desktop/NeuroChat.feb22/March6/Pap1.cellchat.rds")
#saveRDS(cellchat, file = "/Users/usri/Desktop/NeuroChat.feb22/March6/Pap1.cellchat.rds")
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat
## An object of class CellChat created from a single dataset
## 28608 genes.
## 27446 cells.
## CellChat analysis of single cell RNA-seq data!
#We provide a subsetCommunication function that provides easy access to inferred intercellular communications of interest. For example
df.net <- subsetCommunication(cellchat) # it returns a data frame consisting of all inferred cell-to-cell communications at the ligand/receptor level. Set slot.name = "netP" to access inferred communication at the signal path level.
saveRDS(df.net, file = "/Users/usri/Desktop/NeuroChat.feb22/March6/Pap1.df.net.1.rds")
df.net <- subsetCommunication(cellchat, sources.use = c(1,2), targets.use = c(4,5)) # it returns the inferred intercellular communication sent from cell groups 1 and 2 to cell groups 4 and 5
saveRDS(df.net, file = "/Users/usri/Desktop/NeuroChat.feb22/March6/Pap1.df.net.2.rds")
#CellChat calculates communication probabilities at the signaling pathway level by summing the communication probabilities for all ligand-receptor interactions associated with each signaling pathway.
#NOTE: The inferred intercellular communication network for each ligand-receptor pair and each signaling pathway is stored in ‘net’ and ‘netP’ respectively.
cellchat <- computeCommunProbPathway(cellchat)
#CellChat can calculate aggregated inter-cell communication networks by counting the number of links or summarizing communication probabilities. Users can also calculate aggregate networks between a subset of cell groups by setting sources.use and targets.use.
cellchat <- aggregateNet(cellchat)
#execution.time = Sys.time() - ptm
#print(as.numeric(execution.time, units = "secs"))
#CellChat can also visualize aggregated intercellular communication networks. For example, use a circle plot to show the number of interactions or the total interaction strength (weight) between any two groups of cells.
ptm = Sys.time()
groupSize <- as.numeric(table(cellchat@idents))
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")
#Due to the complexity of intercellular communication networks, we can examine the signals sent by each group of cells. Here we also control the parameter edge.weight.max so that we can compare edge weights between different networks.
mat <- cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat)) {
mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
mat2[i, ] <- mat[i, ]
netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}
#After inferring the intercellular communication network, CellChat provides various functions for further data exploration, analysis, and visualization.
######
#It provides multiple ways to visualize intercellular communication networks, including layer diagrams, circle diagrams, chord diagrams, and bubble diagrams.
######
#It provides an easy-to-use tool for extracting and visually inferring high-order information of networks. For example, it allows easy prediction of the main signaling inputs and outputs of cell populations, and how these populations and signals work together to achieve function.
#####
##It enables the quantitative description and comparison of inferred intercellular communication networks using an integrated approach that combines social network analysis, pattern recognition, and manifold learning methods.
##
##The user should define vertex.receiver, which is a numeric vector giving the index of the group of cells in the left part of the hierarchical diagram that is targeted. This hierarchical diagram consists of two parts: the left part shows the autocrine and paracrine signals to certain groups of cells of interest (i.e. the defined vertex.receiver), and the right part shows the autocrine and paracrine signals to the remaining cells in the dataset Cell group. Hierarchical plots therefore provide an informative and intuitive way to visualize autocrine and paracrine communication between groups of cells of interest. For example, when studying cellular communication between fibroblasts and immune cells, the user can define all fibroblast groups as vertex.receiver.
#######
##CellChat provides two functions, netVisual_chord_cell and netVisual_chord_gene, to visualize inter-cell communication for different purposes and at different levels. netVisual_chord_cell is used to visualize intercellular communication between different cell groups (each sector in the chord diagram is a cell group), and netVisual_chord_gene is used to visualize intercellular communication mediated by multiple ligand-receptors or signaling pathways ( Each sector in the chord diagram is a ligand, receptor, or signaling pathway).
#########
##In all visualizations, the edge color is consistent with the source as the sender, and the edge weight is proportional to the interaction strength. Thicker edges indicate stronger signals. In hierarchical and circular plots, circle size is proportional to the number of cells in each cell group. In the hierarchy diagram, solid and open circles represent sources and destinations respectively. In a chord diagram, the inner thinner bar color represents the target receiving the signal from the corresponding outer bar. The size of the inner bar is proportional to the signal strength received by the target. This internal bar helps explain complex chord diagrams. Please note that for some cell groups there are some internal bars without any chords, please ignore it as this is an issue not yet solved by the circlize package.
#######
#You can use netVisual_aggregate to visualize the inferred communication network of a signaling pathway, and netVisual_individual to visualize the inferred communication network of individual ligand-receptor pairs associated with that signaling pathway.
#Here we take the input of a signal path as an example. All signaling pathways showing significant communication can be accessed via cellchat@netP$pathways.
##########
cellchat@netP$pathways
## [1] "NRXN" "Glutamate" "ADGRL" "NCAM" "NRG" "PTPR"
## [7] "CADM" "NEGR" "PTN" "CNTN" "NGL" "PTPRM"
## [13] "UNC5" "SLIT" "Netrin" "GABA-A" "EPHA" "SEMA3"
## [19] "SEMA6" "CDH" "FGF" "GABA-B" "LAMININ" "MAG"
## [25] "PROS" "APP" "FLRT" "CypA" "BMP" "CLDN"
## [31] "GAP" "VEGF" "CDH5" "JAM"
pathways.show <- c("NRXN")
# Hierarchy plot
# Here we define `vertex.receive` so that the left portion of the hierarchy plot shows signaling to fibroblast and the right portion shows signaling to immune cells
vertex.receiver = seq(1,4) # a numeric vector.
netVisual_aggregate(cellchat, signaling = pathways.show, vertex.receiver = vertex.receiver)
# Circle plot for NRXN pathway
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle")
# Chord diagram for NRXN pathway
#CellChat has a separate function netVisual_chord_cell to flexibly visualize the signal network by adjusting different parameters in the circlize package. For example, we can define a named character vector group to create multiple sets of chord diagrams, for example, grouping clusters of cells into different cell types.
par(mfrow=c(1,3))
#netVisual_aggregate(cellchat, signaling = pathways.show, layout = "chord")
# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
## Do heatmap based on a single object
netAnalysis_contribution(cellchat, signaling = pathways.show)
#We can also visualize intercellular communication mediated by individual ligand-receptor pairs. We provide a function extractEnrichedLR to extract all significant interactions (L-R pairs) and associated signaling genes for a given signaling pathway.
pairLR.NRXN <- extractEnrichedLR(cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.NRXN[1,] # show one ligand-receptor pair
# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, vertex.receiver = vertex.receiver)
## [[1]]
# Circle plot
netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
## [[1]]
# Access all the signaling pathways showing significant communications
pathways.show.all <- cellchat@netP$pathways
# check the order of cell identity to set suitable vertex.receiver
levels(cellchat@idents)
## [1] "Astrocyte" "Endothelial Cell"
## [3] "GABA-ergic Neuron" "Glutamatergic Neuron"
## [5] "Microglia" "Oligodendrocyte"
## [7] "Oligodendrocyte Precursor Cell" "Unclassified"
setwd("/Users/usri/Desktop/NeuroChat.feb22/Neurochat.March2/March.4")
vertex.receiver = seq(1,4)
for (i in 1:length(pathways.show.all)) {
# Visualize communication network associated with both signaling pathway and individual L-R pairs
netVisual(cellchat, signaling = pathways.show.all[i], vertex.receiver = vertex.receiver, layout = "hierarchy")
# Compute and visualize the contribution of each ligand-receptor pair to the overall signaling pathway
gg <- netAnalysis_contribution(cellchat, signaling = pathways.show.all[i])
ggsave(filename=paste0(pathways.show.all[i], "_L-R_contribution.pdf"), plot=gg, width = 3, height = 2, units = 'in', dpi = 300)
}
#We can also use netVisual_bubble to display all significant interactions (L-R pairs) from some cell groups to other cell groups.
# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11), remove.isolate = FALSE)
## Comparing communications on a single object
# show all the significant interactions (L-R pairs) associated with certain signaling pathways
pathways.show.all <- cellchat@netP$pathways
pathways.show.all
## [1] "NRXN" "Glutamate" "ADGRL" "NCAM" "NRG" "PTPR"
## [7] "CADM" "NEGR" "PTN" "CNTN" "NGL" "PTPRM"
## [13] "UNC5" "SLIT" "Netrin" "GABA-A" "EPHA" "SEMA3"
## [19] "SEMA6" "CDH" "FGF" "GABA-B" "LAMININ" "MAG"
## [25] "PROS" "APP" "FLRT" "CypA" "BMP" "CLDN"
## [31] "GAP" "VEGF" "CDH5" "JAM"
netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11), signaling = pathways.show.all, remove.isolate = FALSE)
## Comparing communications on a single object
# show all the significant interactions (L-R pairs) based on user's input (defined by `pairLR.use`)
pairLR.use <- extractEnrichedLR(cellchat, signaling = pathways.show.all)
netVisual_bubble(cellchat, sources.use = c(1,7), targets.use = c(1:6), pairLR.use = pairLR.use, remove.isolate = TRUE)
## Comparing communications on a single object
# show all the significant interactions (L-R pairs) associated with certain signaling pathways
netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11), signaling = c("NRXN" ,"Glutamate", "ADGRL", "NRG"), remove.isolate = FALSE)
## Comparing communications on a single object
#Similar to bubble charts, CellChat provides a function netVisual_chord_gene for drawing chord charts.
#Show all interactions (L-R pairs or signaling pathways) from some groups of cells to other groups of cells. Two special cases: one showing all interactions sent from a cell group and the other showing all interactions received by a cell group.
#Demonstrate interactions entered by the user or certain signaling pathways defined by the user.
# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(cellchat, sources.use = 4, targets.use = c(5:11), lab.cex = 0.5,legend.pos.y = 30)
# show all the interactions received by Inflam.DC
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = 8, legend.pos.x = 15)
# show all the significant interactions (L-R pairs) associated with certain signaling pathways
#netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("NRXN", "Glutamate", "ADGRL", "NCAM", "NRG", "PTPR", "CADM"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("NRXN"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("Glutamate"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("ADGRL"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("NCAM"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("NRG"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("PTPR"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("CADM"),legend.pos.x = 8)
#We can use a Seurat wrapper function plotGeneExpression to plot the gene expression distribution of signaling genes associated with L-R pairs or signaling pathways.
plotGeneExpression(cellchat, signaling = "PTPR")
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
plotGeneExpression(cellchat, signaling = "NCAM")
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
plotGeneExpression(cellchat, signaling = "NCAM", enriched.only = FALSE)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
#CellChat allows easy identification of the main senders, receivers, mediators, and influencers in intercellular communication networks by calculating several network centrality measures for each cell group. Specifically, we used metrics in weighted directed networks, including out-degree, in-degree, flow centrality, and information centrality, to identify the main senders, receivers, mediators, and influencers of inter-cell communication, respectively. In a weighted directed network with calculated communication probabilities as weights, out-degree (calculated as the sum of communication probabilities sent by a cell group) and in-degree (calculated as the sum of communication probabilities received by a cell group), can be respectively Used to identify major cellular senders and receivers of signaling networks.
# 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
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)
#you can also intuitive way to visualize the main senders (sources) and receivers (targets) in two-dimensional space using scatter plots.
# 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 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("NCAM", "PTPR", "NCAM", "NRXN"))
## Signaling role analysis on the cell-cell communication network from user's input
#> Signaling role analysis on the cell-cell communication network from user's input
gg1 + gg2
#You can also answer the question of which signals contribute most to the emitted or received signals of certain groups of cells.
# 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")
ht1 + ht2
# Signaling role analysis on the cell-cell communication networks of interest
ht <- netAnalysis_signalingRole_heatmap(cellchat, signaling = c("NCAM", "PTPR", "NCAM", "NRXN"))
ht
#In addition to exploring the detailed communication of individual pathways, an important question is how multiple cell groups and signaling pathways work in coordination. CellChat uses pattern recognition methods to identify global communication patterns.
#As the number of patterns increases, redundant patterns may appear, making it difficult to interpret communication patterns. We select five modes by default. Typically, patterns greater than 2 are biologically meaningful. Additionally, we provide a selectK function to infer the number of patterns based on two metrics implemented in the NMF R package, including Cophenetic and Silhouette. These two indicators measure the stability of a specific number of patterns based on hierarchical clustering of the consensus matrix. For a range of mode numbers, the appropriate mode number is where the Cophenetic and Silhouette values start to suddenly drop off.
#Outgoing patterns reveal how sending cells (i.e., cells that serve as signal sources) coordinate with each other and how they coordinate with certain signaling pathways to drive communication.
##To visualize the association of potential patterns with groups of cells and ligand-receptor pairs or signaling pathways, we used alluvial plots. We first normalize each row of W and each column of H to [0,1], then set elements in W and H to zero if they are less than 0.5. Such thresholding allows revealing the most enriched cell groups and signaling pathways with each inferred pattern, i.e., each cell group or signaling pathway is associated with only one inferred pattern. These thresholded matrices W and H are used as inputs to create river plots.
library(NMF)
## Loading required package: registry
## Loading required package: rngtools
## Loading required package: cluster
## NMF - BioConductor layer [OK] | Shared memory capabilities [NO: bigmemory] | Cores 2/2
## To enable shared memory capabilities, try: install.extras('
## NMF
## ')
##
## Attaching package: 'NMF'
## The following objects are masked from 'package:igraph':
##
## algorithm, compare
library(ggalluvial)
selectK(cellchat, pattern = "outgoing")
# Both Cophenetic and Silhouette values begin to drop suddenly when the number of outgoing patterns is 3.
nPatterns = 3
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = nPatterns)
# river plot
netAnalysis_river(cellchat, pattern = "outgoing")
## Please make sure you have load `library(ggalluvial)` when running this function
# dot plot
netAnalysis_dot(cellchat, pattern = "outgoing")
#Shows how target cells (i.e. cells that are signal receivers) coordinate with each other and how they coordinate with certain signaling pathways to respond to received signals.
selectK(cellchat, pattern = "incoming")
nPatterns = 3
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = nPatterns)
# river plot
netAnalysis_river(cellchat, pattern = "incoming")
## Please make sure you have load `library(ggalluvial)` when running this function
# dot plot
netAnalysis_dot(cellchat, pattern = "incoming")
#Additionally, CellChat is able to quantify the similarities between all significant signaling pathways and then group them based on their cellular communication network similarities. Grouping can be based on functional similarity or structural similarity.
##
#Functional similarity: A high degree of functional similarity indicates that the primary sender and receiver are similar, which can be interpreted as two signaling pathways or two ligand-receptor pairs exhibiting similar and/or redundant roles. Functional similarity analysis requires identical cell population composition between the two data sets.
##
#Structural similarity: Structural similarity is used to compare their signaling network structures, regardless of the similarity of the sender and receiver.
#Identification of signal groups based on functional similarity.
cellchat <- computeNetSimilarity(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
## Manifold learning of the signaling networks for a single dataset
#> Manifold learning of the signaling networks for a single dataset
cellchat <- netClustering(cellchat, type = "functional")
## Classification learning of the signaling networks for a single dataset
#> Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(cellchat, type = "functional", label.size = 3.5)
netVisual_embeddingZoomIn(cellchat, type = "functional", nCol = 2)
##Identification of signal groups based on structural similarities.
cellchat <- computeNetSimilarity(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
## Manifold learning of the signaling networks for a single dataset
#> Manifold learning of the signaling networks for a single dataset
cellchat <- netClustering(cellchat, type = "structural")
## Classification learning of the signaling networks for a single dataset
#> Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(cellchat, type = "structural", label.size = 3.5)
netVisual_embeddingZoomIn(cellchat, type = "structural", nCol = 2)
#saveRDS(cellchat, file = "cellchat_snRNAseq.Pap1.rds")