Multiple Myeloma (MM) CellChat Analysis

pre-CART: DR vs PD

## 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
## 
## 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
## ── Attaching core tidyverse packages ─────────────────── tidyverse 2.0.0.9000 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ✔ readr     2.1.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ lubridate::%--%()        masks igraph::%--%()
## ✖ BiocGenerics::Position() masks ggplot2::Position(), base::Position()
## ✖ tibble::as_data_frame()  masks igraph::as_data_frame(), dplyr::as_data_frame()
## ✖ Biobase::combine()       masks BiocGenerics::combine(), dplyr::combine()
## ✖ purrr::compose()         masks igraph::compose()
## ✖ tidyr::crossing()        masks igraph::crossing()
## ✖ dplyr::filter()          masks stats::filter()
## ✖ dplyr::lag()             masks stats::lag()
## ✖ purrr::simplify()        masks igraph::simplify()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## [1] "/research/bsi/projects/PI/tertiary/Lin_Yi_yxl07/s215605.10xsc_human_lym_extend/mm_map2ref_bmref/output_mm_all/"
## [1] "file exists!"
## [1] "/research/bsi/projects/PI/tertiary/Lin_Yi_yxl07/s215605.10xsc_human_lym_extend/mm_map2ref_bmref/output_mm_all/cellchat_comparison/preDR-vs-PD/"
## [1] "/research/bsi/projects/PI/tertiary/Lin_Yi_yxl07/s215605.10xsc_human_lym_extend/mm_map2ref_bmref/output_mm_all/cellchat/pre_PD/cellchat_pre-PDr4.2.rds"
## [1] 2
## Merge the following slots: 'data.signaling','images','net', 'netP','meta', 'idents', 'var.features' , 'DB', and 'LR'.
## An object of class CellChat created from a merged object with multiple datasets 
##  933 signaling genes.
##  42850 cells. 
## CellChat analysis of single cell RNA-seq data!

Part I: Predict general principles of cell-cell communication

CellChat starts with the big picture to predict general principles of cell-cell communication. When comparing cell-cell communication among multiple biological conditions, it can answer the following biological questions:

-Whether the cell-cell communication is enhanced or not

-The interaction between which cell types is significantly changed

-How the major sources and targets change from one condition to another

Compare the number of interactions and interaction strength among different cell populations To identify the interaction between which cell populations showing significant changes, CellChat compares the number of interactions and interaction strength among different cell populations.

Differential number of interactions or interaction strength among different cell populations The differential number of interactions or interaction strength in the cell-cell communication network between two datasets can be visualized using circle plot, where red (or blue) colored edges represent increased (or decreased) signaling in the second dataset compared to the first one.

Show differential number of interactions or interaction strength in a greater details using a heatmap. The top colored bar plot represents the sum of column of values displayed in the heatmap (incoming signaling). The right colored bar plot represents the sum of row of values (outgoing signaling). In the colorbar, red (or blue) represents increased (or decreased) signaling in the second dataset compared to the first one.

## Do heatmap based on a merged object 
## 
## Do heatmap based on a merged object

To better control the node size and edge weights of the inferred networks across different datasets, we compute the maximum number of cells per cell group and the maximum number of interactions (or interaction weights) across all datasets.

Compare the major sources and targets in 2D space Comparing the outgoing and incoming interaction strength in 2D space allows ready identification of the cell populations with significant changes in sending or receiving signals between different datasets.

Part II: Identify the conserved and context-specific signaling pathways

CellChat then can identify signaling networks with larger (or less) difference, signaling groups, and the conserved and context-specific signaling pathways based on their cell-cell communication networks among multiple biological conditions.

Identify signaling networks with larger (or less) difference as well as signaling groups based on their functional/structure similarity CellChat performs joint manifold learning and classification of the inferred communication networks based on their functional and topological similarity. NB: Such analysis is applicable to more than two datasets.

Functional similarity: High degree of functional similarity indicates major senders and receivers are similar, and it can be interpreted as the two signaling pathways or two ligand-receptor pairs exhibit similar and/or redundant roles. NB: Functional similarity analysis is not applicable to multiple datsets with different cell type composition.

Structural similarity: A structural similarity was used to compare their signaling network structure, without considering the similarity of senders and receivers. NB: Structural similarity analysis is applicable to multiple datsets with the same cell type composition or the vastly different cell type composition.

#Here we can run the manifold and classification learning analysis based on the functional similarity because the two datasets have the the same cell type composition.

#Identify signaling groups based on their functional similarity

#Identify signaling groups based on structure similarity

Compute and visualize the pathway distance in the learned joint manifold

#We can identify the signaling networks with larger (or less) difference based on their Euclidean distance in the shared two-dimensions space. Larger distance implies larger difference of the communication networks between two datasets in terms of either functional or structure similarity. NB: We only compute the distance of overlapped signaling pathways between two datasets. Those signaling pathways that are only identified in one dataset are not considered here. If there are more than three datasets, one can do pairwise comparisons by defining comparison in the function rankSimilarity.

Identify and visualize the conserved and context-specific signaling pathways

By comparing the information flow/interaction strengh of each signaling pathway, we can identify signaling pathways, (i) turn off, (ii) decrease, (iii) turn on or (iv) increase, by change their information flow at one condition as compared to another condition.

Compare the overall information flow of each signaling pathway We can identify the conserved and context-specific signaling pathways by simply comparing the information flow for each signaling pathway, which is defined by the sum of communication probability among all pairs of cell groups in the inferred network (i.e., the total weights in the network).

This bar graph can be plotted in a stacked mode or not. Significant signaling pathways were ranked based on differences in the overall information flow within the inferred networks between the two conditions.

Compare outgoing (or incoming) signaling associated with each cell population

The above analysis summarize the information from the outgoing and incoming signaling together. We can also compare the outgoing (or incoming) signaling pattern between two datasets, allowing to identify signaling pathways/ligand-receptors that exhibit different signaling patterns.

We can combine all the identified signaling pathways from different datasets and thus compare them side by side, including outgoing signaling, incoming signaling and overall signaling by aggregating outgoing and incoming signaling together. NB: rankNet also shows the comparison of overall signaling, but it does not show the signaling strength in specific cell populations.

## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.14.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))
## ========================================

Part III: Identify the upgulated and down-regulated signaling ligand-receptor pairs

Identify dysfunctional signaling by comparing the communication probabities We can compare the communication probabilities mediated by ligand-receptor pairs from some cell groups to other cell groups. This can be done by setting comparison in the function netVisual_bubble.

## Comparing communications on a merged object

Moreover, identify the up-regulated (increased) and down-regulated (decreased) signaling ligand-receptor pairs in one dataset compared to the other dataset. This can be done by specifying max.dataset and min.dataset in the function netVisual_bubble. The increased signaling means these signaling have higher communication probability (strength) in one dataset compared to the other dataset.

## Comparing communications on a merged object 
## 
## Comparing communications on a merged object

NB: The ligand-receptor pairs shown in the bubble plot can be accessed via signaling.LSIncreased = gg1$data.

Identify dysfunctional signaling by using differential expression analysis

Alternative, we can identify the upgulated and down-regulated signaling ligand-receptor pairs based on the differential gene expression analysis. Specifically, we perform differential expression analysis between two biological conditions for each cell group, and then obtain the up-regulated and down-regulated signaling based on the fold change of ligands in the sender cells and receptors in the receiver cells. Such analysis can be done as follows.

## Use the joint cell labels from the merged CellChat object

Since the signaling genes in the net.up and net.down might be complex with multi-subunits, further deconvolution is done to obtain the individual signaling genes. We then visualize the upgulated and down-regulated signaling ligand-receptor pairs using bubble plot or chord diagram.

## Comparing communications on a merged object 
## 
## Comparing communications on a merged object

Visualize the upgulated and down-regulated signaling ligand-receptor pairs using Chord diagram

#Visualize the enriched ligands, signaling,or ligand-receptor pairs in one condition compared to another condition using wordcloud

Part IV: Visually compare cell-cell communication using Hierarchy plot, Circle plot or Chord diagram

Visualize the cell-cell communication network using Hierarchy plot, Circle plot or Chord diagram.

Edge color/weight, node color/size/shape: In all visualization plots, edge colors are consistent with the sources as sender, and edge weights are proportional to the interaction strength. Thicker edge line indicates a stronger signal. In the Hierarchy plot and Circle plot, circle sizes are proportional to the number of cells in each cell group. In the hierarchy plot, solid and open circles represent source and target, respectively. In the Chord diagram, the inner thinner bar colors represent the targets that receive signal from the corresponding outer bar. The inner bar size is proportional to the signal strength received by the targets. Such inner bar is helpful for interpreting the complex chord diagram. Note that there exist some inner bars without any chord for some cell groups, please just igore it because this is an issue that has not been addressed by circlize package.

## Comparing communications on a merged object 
## 
## Comparing communications on a merged object

separate plot

## Comparing communications on a merged object 
## 
## Comparing communications on a merged object 
## 
## Comparing communications on a merged object 
## 
## Comparing communications on a merged object

Using chord diagram, CellChat provides two functions netVisual_chord_cell and netVisual_chord_gene for visualizing cell-cell communication with different purposes and different levels. netVisual_chord_cell is used for visualizing the cell-cell communication between different cell groups (where each sector in the chord diagram is a cell group), and netVisual_chord_gene is used for visualizing the cell-cell communication mediated by mutiple ligand-receptors or signaling pathways (where each sector in the chord diagram is a ligand, receptor or signaling pathway.)

NB: Please ignore the note when generating the plot such as “Note: The first link end is drawn out of sector ‘MIF’.”. If the gene names are overlapped, you can adjust the argument small.gap by decreasing the value.

Save the merged CellChat object

## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux 8.6 (Ootpa)
## 
## Matrix products: default
## BLAS:   /usr/lib64/libblas.so.3.8.0
## LAPACK: /usr/lib64/liblapack.so.3.8.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C         LC_TIME=C           
##  [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=C       
##  [7] LC_PAPER=C           LC_NAME=C            LC_ADDRESS=C        
## [10] LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C 
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ComplexHeatmap_2.14.0 lubridate_1.9.2       forcats_1.0.0        
##  [4] stringr_1.5.0         purrr_1.0.1           readr_2.1.4          
##  [7] tidyr_1.3.0           tibble_3.2.1          tidyverse_2.0.0.9000 
## [10] patchwork_1.1.2       CellChat_1.6.1        synchronicity_1.3.5  
## [13] bigmemory_4.6.1       Biobase_2.58.0        BiocGenerics_0.44.0  
## [16] ggplot2_3.4.2         igraph_1.5.0          dplyr_1.1.2          
## 
## loaded via a namespace (and not attached):
##   [1] bigmemory.sri_0.1.6  colorspace_2.1-0     ggsignif_0.6.4      
##   [4] rjson_0.2.21         circlize_0.4.15      GlobalOptions_0.1.2 
##   [7] BiocNeighbors_1.16.0 clue_0.3-64          rstudioapi_0.14     
##  [10] farver_2.1.1         ggpubr_0.6.0         listenv_0.9.0       
##  [13] ggrepel_0.9.3        RSpectra_0.16-1      fansi_1.0.4         
##  [16] codetools_0.2-18     doParallel_1.0.17    cachem_1.0.8        
##  [19] knitr_1.43           jsonlite_1.8.7       Cairo_1.6-0         
##  [22] broom_1.0.5          gridBase_0.4-7       cluster_2.1.4       
##  [25] png_0.1-8            BiocManager_1.30.21  compiler_4.2.2      
##  [28] backports_1.4.1      Matrix_1.5-4.1       fastmap_1.1.1       
##  [31] cli_3.6.1            htmltools_0.5.5      tools_4.2.2         
##  [34] coda_0.19-4          gtable_0.3.3         glue_1.6.2          
##  [37] reshape2_1.4.4       Rcpp_1.0.10          carData_3.0-5       
##  [40] statnet.common_4.9.0 jquerylib_0.1.4      NMF_0.26            
##  [43] vctrs_0.6.3          svglite_2.1.1        iterators_1.0.14    
##  [46] ggalluvial_0.12.5    xfun_0.39            globals_0.16.2      
##  [49] network_1.18.1       timechange_0.2.0     lifecycle_1.0.3     
##  [52] irlba_2.3.5.1        rngtools_1.5.2       rstatix_0.7.2       
##  [55] future_1.32.0        scales_1.2.1         hms_1.1.3           
##  [58] parallel_4.2.2       RColorBrewer_1.1-3   yaml_2.3.7          
##  [61] reticulate_1.30      pbapply_1.7-2        ggnetwork_0.5.12    
##  [64] sass_0.4.6           stringi_1.7.12       highr_0.10          
##  [67] S4Vectors_0.36.2     foreach_1.5.2        BiocParallel_1.32.6 
##  [70] shape_1.4.6          rlang_1.1.1          pkgconfig_2.0.3     
##  [73] systemfonts_1.0.4    matrixStats_1.0.0    evaluate_0.21       
##  [76] lattice_0.20-45      labeling_0.4.2       cowplot_1.1.1       
##  [79] tidyselect_1.2.0     parallelly_1.36.0    plyr_1.8.8          
##  [82] magrittr_2.0.3       R6_2.5.1             magick_2.7.4        
##  [85] IRanges_2.32.0       generics_0.1.3       sna_2.7-1           
##  [88] pillar_1.9.0         withr_2.5.0          abind_1.4-5         
##  [91] future.apply_1.11.0  crayon_1.5.2         car_3.1-2           
##  [94] uuid_1.1-0           utf8_1.2.3           tzdb_0.4.0          
##  [97] rmarkdown_2.22       GetoptLong_1.0.5     FNN_1.1.3.2         
## [100] digest_0.6.32        stats4_4.2.2         munsell_0.5.0       
## [103] registry_0.5-1       bslib_0.5.0