CellChat_NetworkCentrality.Rmd

---
title: "CellChat_NetworkCentrality"
author: "Daniel Lee"
date: "2023-10-13"
output: pdf_document
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = TRUE,  
                      warning = FALSE, message = FALSE, 
                      fig.align = "center",
                      R.options = list(max.print=100))
```

```{r}
# running renv
# install.packages("renv")
library(renv)
# renv::init()
# renv::snapshot()
# renv::status()
# renv::clean()
# renv::restore()
# renv::history()
```

```{r}
library(CellChat)
library(Seurat)
```

```{r}
# Load your data
setwd("/Volumes/biologos/Bioinformatics/Using-CellChat-for-Ligand-Receptor-Network-Analysis/Data")
myeloid <- readRDS("myeloid_srt_int.rds")
lymphoid <- readRDS("lymphoid_srt_int.rds")

# Assign cell types to myeloid clusters
myeloid$cell_type <- factor(
  recode(myeloid$seurat_clusters,
         `0` = "IL1B+ APOE+ Macrophages",
         `1` = "IL1B+ APOE- Macrophages",
         `2` = "Monocytes",
         `3` = "KIT+ Mast cells",
         `4` = "Dendritic cells",
         `5` = "IL1B- APOE+ Macrophages",
         `6` = "LYVE1+ Macrophages",
         `7` = "IFN-associated Myeloid cells"
  )
)

# Assign cell types to lymphoid clusters
lymphoid$cell_type <- factor(
  recode(lymphoid$seurat_clusters,
         `0` = "IL7R+ CD4+ Effector T cells",
         `1` = "GZMK+ CD8+ T cells",
         `2` = "GZMH+ CD8+ T cells",
         `3` = "Natural Killer cells",
         `4` = "CXCL8+ T cells",
         `5` = "CD79A+ B cells type 1",
         `6` = "CD79A+ B cells type 2",
         `7` = "TIGIT+ CD4+ Treg cells",
         `8` = "Plasma cells"
  )
)

# Designate the specific clusters
# myeloid_clusters <- c("IL1B+ APOE+ Macrophages",
#                       "IL1B+ APOE- Macrophages",
#                       "Dendritic cells",
#                       "IL1B- APOE+ Macrophages",
#                       "LYVE1+ Macrophages")
# 
# lymphoid_clusters <- c("IL7R+ CD4+ Effector T cells",
#                        "GZMK+ CD8+ T cells",
#                        "GZMH+ CD8+ T cells",
#                        "CXCL8+ T cells",
#                        "TIGIT+ CD4+ Treg cells")

# Subset the Seurat objects based on the cell type
# myeloid_subset <- subset(myeloid, subset = cell_type %in% myeloid_clusters)
#
# lymphoid_subset <- subset(lymphoid, subset = cell_type %in% lymphoid_clusters)

# Check the subset step
# print(unique(myeloid_subset$cell_type))
# print(unique(lymphoid_subset$cell_type))

# Offset the seurat_clusters in the lymphoid dataset
lymphoid$seurat_clusters <- as.numeric(lymphoid$seurat_clusters) + max(as.numeric(myeloid$seurat_clusters)) -1
unique(lymphoid$seurat_clusters)

# Combining the two datasets into a single Seurat object
combined_data <- merge(x = lymphoid, y = myeloid)
unique(combined_data@meta.data$seurat_clusters)
```

## Sample subsetting
```{r}
# Define femoral v.s. carotid
carotid_ids <- c("carotid 1 & 2", "carotid 3", "carotid 4")
femoral_ids <- setdiff(unique(combined_data@meta.data$orig.ident), carotid_ids)
# Subset the combined data for each group
carotid_data <- combined_data[, combined_data@meta.data$orig.ident %in% carotid_ids]
femoral_data <- combined_data[, combined_data@meta.data$orig.ident %in% femoral_ids]

# Get the unique cell_type labels in the order of the seurat_clusters
ordered_cell_types <- unique(combined_data$cell_type[order(as.numeric(combined_data$seurat_clusters))])
# Set the levels of the cell_type factor based on this order
carotid_data$cell_type <- factor(carotid_data$cell_type, levels = ordered_cell_types)
femoral_data$cell_type <- factor(femoral_data$cell_type, levels = ordered_cell_types)

# Additional cell marker checks
p1 <- VlnPlot(carotid_data, features = "MIF", group.by = "cell_type") + ggtitle("MIF Expression in Carotid Data")
print(p1)
p2 <- VlnPlot(femoral_data, features = "MIF", group.by = "cell_type") + ggtitle("MIF Expression in Femoral Data")
print(p2)
p3 <- VlnPlot(carotid_data, features = "SPP1", group.by = "cell_type") + ggtitle("SPP1 Expression in Carotid Data")
print(p3)
p4 <- VlnPlot(femoral_data, features = "SPP1", group.by = "cell_type") + ggtitle("SPP1 Expression in Femoral Data")
print(p4)

```

```{r}
# Plot for MIF in myeloid cells
p5 <- VlnPlot(myeloid, features = "MIF", group.by = "cell_type", split.by = "type") + ggtitle("MIF Expression in Myeloid Cells")

# Plot for SPP1 in lymphoid cells
p6 <- VlnPlot(myeloid, features = "SPP1", group.by = "cell_type", split.by = "type") + ggtitle("SPP1 Expression in Myeloid Cells")

p5+p6

# Plot for CD44 in lymphoid cells
p7 <- VlnPlot(lymphoid, features = "CD44", group.by = "cell_type", split.by = "type") + ggtitle("CD44 Expression in Lymphoid Cells")

# Plot for CD74 in lymphoid cells
p8 <- VlnPlot(lymphoid, features = "CD74", group.by = "cell_type", split.by = "type") + ggtitle("CD74 Expression in Lymphoid Cells")

# Plot for CXCR4 in lymphoid cells
p9 <- VlnPlot(lymphoid, features = "CXCR4", group.by = "cell_type", split.by = "type") + ggtitle("CXCR4 Expression in Lymphoid Cells")
p7+p8+p9
```

```{r}
library(dplyr)
library(tidyr)
library(gridExtra)
library(grid)
library(tibble)

# List of genes of interest
genes_of_interest <- c("MIF", "SPP1", "CD44", "CD74", "CXCR4")

# Use AverageExpression to get mean expression by cell type and type
avg_expression <- AverageExpression(combined_data, slot = "data", features = genes_of_interest, group.by = c("cell_type", "type"))
# trace(AverageExpression, edit = T)

# Convert the averaged expression values to a data frame and reshape for easy manipulation
df_avg_expression <- as.data.frame(avg_expression$RNA) %>%
  rownames_to_column("Gene") %>%
  t() %>%
  as.data.frame()

# Make the first row the header
colnames(df_avg_expression) <- as.character(df_avg_expression[1, ])

# Remove the first row
df_avg_expression <- df_avg_expression[-1, ]

# Set the working directory 
setwd("/Volumes/biologos/Bioinformatics/Using-CellChat-for-Ligand-Receptor-Network-Analysis/Visuals/Comparison")

# Save the grid object as PNG
write.csv(df_avg_expression, "table_plot_revised.csv", row.names = T)
```

AverageExpression Function Validation
```{r}
# Apply expm1 transformation to RNA@data 
transformed_data <- expm1(combined_data@assays$RNA@data[genes_of_interest,])

# Extract metadata to group by "cell_type" and "type":
metadata <- combined_data@meta.data[, c("cell_type", "type")]

# Function to calculate mean expressions for a given gene after transformation
calculate_mean <- function(gene_name) {
  gene_expression <- transformed_data[gene_name,]
  combined <- data.frame(Gene_Expression=gene_expression, metadata)
  
  mean_values <- combined %>%
    group_by(cell_type, type) %>%
    summarise(mean_expression = mean(Gene_Expression, na.rm = TRUE)) %>%
    mutate(Gene = gene_name)
  
  return(mean_values)
}

# Use lapply to get mean values for all genes and bind them into a single data frame
all_mean_values <- do.call(rbind, lapply(genes_of_interest, calculate_mean))

# Spread the table to have genes as columns and mean expression values as rows
clean_table <- all_mean_values %>%
  spread(key = Gene, value = mean_expression)

# View the clean summary table
print(clean_table)
```

```{r}
library(tidyverse)
library(gt)
# combined_data_sct <- SCTransform(combined_data, vst.flavor = "v2", verbose = FALSE)

table(Idents(combined_data), combined_data$seurat_clusters)
combined_data$seurat_clusters <- as.factor(combined_data$seurat_clusters)
Idents(object = combined_data) <- "seurat_clusters"

# Getting the unique cell types in your data
unique_cell_types <- unique(combined_data$cell_type)

# Initialize a list to store the results for each cell type
all_markers_list <- list()

# Loop through each unique cell type
for (cell_type in unique_cell_types) {
  # Subset data for the current cell type in loop
  current_data <- subset(combined_data, cell_type == !!cell_type)
  
  # Run FindMarkers for the current subset of data and store the result in the list
  all_markers_list[[cell_type]] <- FindMarkers(
    object = current_data,
    ident.1 = "carotid",
    ident.2 = "femoral",
    features = c("MIF", "SPP1", "CD44", "CD74", "CXCR4"),
    logfc.threshold = 0,
    only.pos = FALSE,
    group.by = "type"
  )

  # Optional: print a message to keep track of progress
  print(paste("Completed analysis for", cell_type))
}

all_markers <- bind_rows(all_markers_list, .id = "cell_type")
all_markers$gene <- rownames(all_markers)
rownames(all_markers) <- NULL

# Clean the data frame
cleaned_all_marker <- all_markers %>%
  # Separate gene names and numbers into two separate columns
  separate(gene, into = c("gene", "number"), sep = "\\.\\.\\.") %>%
  # Remove the "number" column
  select(-number) %>%
  # Rearrange the columns
  select(cell_type, gene, avg_log2FC, p_val, p_val_adj, pct.1, pct.2)

# Create a beautiful table
table <- cleaned_all_marker %>%
  gt() %>%
  tab_header(
    title = "Fold-Change Comparisons of SPP1 & MIF Pathway Genes"
  ) %>%
  cols_label(
    cell_type = "Cell Type",
    gene = "Gene",
    p_val = "P-Value",
    avg_log2FC = "Avg log2 Fold Change",
    pct.1 = "% of Carotid Cells",
    pct.2 = "% of Femoral Cells",
    p_val_adj = "Adjusted P-Value"
  ) %>%
  fmt_scientific(
    columns = vars(p_val, p_val_adj),
    scale_by = 1,
    decimals = 3
  ) %>%
  fmt_number(
    columns = vars(avg_log2FC, `pct.1`, `pct.2`),
    decimals = 3
  ) %>%
  data_color(
    columns = vars(p_val_adj),
    colors = function(x) {
      ifelse(x < 0.05, "lightgreen", "white")
    }
  )

# Set the working directory 
setwd("/Volumes/biologos/Bioinformatics/Using-CellChat-for-Ligand-Receptor-Network-Analysis/Visuals/Comparison")

# Save the table as a file
gtsave(table, filename = "fold-change_table.png")

```

```{r}
# Get the number of cells for each `cell_type` and `type`
# cell_counts <- metadata %>%
#   group_by(cell_type, type) %>%
#   summarise(count = n())

# Set the working directory 
# setwd("/Volumes/biologos/Bioinformatics/Using-CellChat-for-Ligand-Receptor-Network-Analysis/Visuals/Comparison")

# Save the cell counts as a CSV file
# write.csv(cell_counts, "cell_counts.csv", row.names = FALSE)

# Print the cell counts
# print(cell_counts)
```

```{r}
# Wilcox Test
library(ggpubr)
library(gridExtra)

p5 + stat_compare_means(# label = "p.signif",
                         # aes(label=..p.adj..),
                         symnum.args = list(cutpoints = c(0, 0.05, Inf), 
                                            symbols = c("*", "")),
                         method = "wilcox.test",
                         paired = FALSE,
                         hide.ns = F)

p6 + stat_compare_means(# label = "p.signif",
                         # aes(label=..p.adj..),
                         symnum.args = list(cutpoints = c(0, 0.05, Inf), 
                                            symbols = c("*", "")),
                         method = "wilcox.test",
                         paired = FALSE,
                         hide.ns = F)

p7 + stat_compare_means(# label = "p.signif",
                         # aes(label=..p.adj..),
                         symnum.args = list(cutpoints = c(0, 0.05, Inf), 
                                            symbols = c("*", "")),
                         method = "wilcox.test",
                         paired = FALSE,
                         hide.ns = F)

p8 + stat_compare_means(# label = "p.signif",
                        # aes(label=..p.adj..),
                         symnum.args = list(cutpoints = c(0, 0.05, Inf), 
                                            symbols = c("*", "")),
                         method = "wilcox.test",
                         paired = FALSE,
                         hide.ns = F)

p9 + stat_compare_means(# label = "p.signif",
                        # aes(label=..p.adj..),
                         symnum.args = list(cutpoints = c(0, 0.05, Inf), 
                                            symbols = c("*", "")),
                         method = "wilcox.test",
                         paired = FALSE,
                         hide.ns = F)

# add_stat_comparison <- function(plot) {
#   return(
#     plot + 
#       stat_compare_means(label = "p.signif",
#                          # aes(label=..p.adj..),
#                          symnum.args = list(cutpoints = c(0, 0.05, Inf), 
#                                             symbols = c("*", "")),
#                          method = "wilcox.test",
#                          paired = FALSE,
#                          hide.ns = F)
#   )
# }
# 
# # Using lapply to apply the function to each plot object
# plot_list <- list(p5, p6, p7, p8, p9)
# plot_list_updated <- lapply(plot_list, add_stat_comparison)
# 
# # Assigning the modified plots back to the original variables
# p5 <- plot_list_updated[[1]]
# p6 <- plot_list_updated[[2]]
# p7 <- plot_list_updated[[3]]
# p8 <- plot_list_updated[[4]]
# p9 <- plot_list_updated[[5]]
```

```{r}
# Set the working directory 
setwd("/Volumes/biologos/Bioinformatics/Using-CellChat-for-Ligand-Receptor-Network-Analysis/Visuals/Comparison")

# Combine and save p5 + p6
# combined_plot_1 <- arrangeGrob(p5, p6, ncol=2)
# ggsave("myeloid_violinplot.png", combined_plot_1, dpi=600, width=30, height=10)
library(showtext)
showtext_auto()
showtext_opts(dpi=600)
ggsave("myeloid_violinplot_MIF.png", p5, dpi=600, width=14, height=5)
ggsave("myeloid_violinplot_SPP1.png", p6, dpi=600, width=14, height=5)

# Combine and save p7 + p8 + p9
# combined_plot_2 <- arrangeGrob(p7, p8, p9, ncol=3)
# ggsave("lymphoid_violinplot.png", combined_plot_2, dpi=600, width=42, height=10)
ggsave("lymphoid_violinplot_CD44.png", p7, dpi=600, width=15, height=4.5)
ggsave("lymphoid_violinplot_CD74.png", p8, dpi=600, width=15, height=4.5)
ggsave("lymphoid_violinplot_CXCR4.png", p9, dpi=600, width=15, height=4.5)

```

## Carotid sample analysis
```{r}
# carotid
cellchat_carotid <- createCellChat(object = carotid_data, meta = carotid_data@meta.data, group.by = "cell_type")

# Set the ligand-receptor interaction database
CellChatDB.use <- CellChatDB.human

# Using a subset of CellChatDB for cell-cell communication analysis
CellChatDB.use <- subsetDB(CellChatDB.human, search = "Secreted Signaling")

# Add the Secreted Signaling database in the CellChat object
cellchat_carotid@DB <- CellChatDB.use

```

```{r}
# Pre-process the expression data
cellchat_carotid <- subsetData(cellchat_carotid)
cellchat_carotid <- identifyOverExpressedGenes(cellchat_carotid)
cellchat_carotid <- identifyOverExpressedInteractions(cellchat_carotid)

```

```{r}
# Project gene expression data onto protein-protein interaction (PPI)
cellchat_carotid <- projectData(cellchat_carotid, PPI.human)
```

```{r}
# Compute communication probabilities again
cellchat_carotid <- computeCommunProb(cellchat_carotid, raw.use = TRUE)
```

```{r}
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat_carotid <- filterCommunication(cellchat_carotid, min.cells = 10)

```

```{r}
# Infer the cell-cell communication at a signaling pathway level
cellchat_carotid <- computeCommunProbPathway(cellchat_carotid, thresh = 0.05)

```

```{r}
# Calculate the aggregated cell-cell communication network
cellchat_carotid <- aggregateNet(cellchat_carotid, remove.isolate = TRUE, thresh = 0.05)
cellchat_carotid@net$count
cellchat_carotid@net$weight

```

## Femoral sample analysis
```{r}
# femoral group
cellchat_femoral <- createCellChat(object = femoral_data, meta = femoral_data@meta.data, group.by = "cell_type")

# Set the ligand-receptor interaction database
CellChatDB.use <- CellChatDB.human

# Using a subset of CellChatDB for cell-cell communication analysis
CellChatDB.use <- subsetDB(CellChatDB.human, search = "Secreted Signaling")

# Add the Secreted Signaling database in the CellChat object
cellchat_femoral@DB <- CellChatDB.use

```

```{r}
# Pre-process the expression data
cellchat_femoral <- subsetData(cellchat_femoral)
cellchat_femoral <- identifyOverExpressedGenes(cellchat_femoral)
cellchat_femoral <- identifyOverExpressedInteractions(cellchat_femoral)

```

```{r}
# Project gene expression data onto protein-protein interaction (PPI)
cellchat_femoral <- projectData(cellchat_femoral, PPI.human)
```

```{r}
# Compute communication probabilities again
cellchat_femoral <- computeCommunProb(cellchat_femoral, raw.use = TRUE)

```

```{r}
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat_femoral <- filterCommunication(cellchat_femoral, min.cells = 10)

```

```{r}
# Infer the cell-cell communication at a signaling pathway level
cellchat_femoral <- computeCommunProbPathway(cellchat_femoral, thresh = 0.05)

```

```{r}
# Calculate the aggregated cell-cell communication network
cellchat_femoral <- aggregateNet(cellchat_femoral, remove.isolate = TRUE, thresh = 0.05)
cellchat_femoral@net$count
cellchat_femoral@net$weight

```
## Save CellChat Objects
```{r}
# set the working directory for retrieving CellChat objects
setwd("/Volumes/biologos/Bioinformatics/Using-CellChat-for-Ligand-Receptor-Network-Analysis/CellChatObjects")

# Save CellChat objects
saveRDS(cellchat_carotid, file = "fullcarotid_cellchat.rds")
saveRDS(cellchat_femoral, file = "fullfemoral_cellchat.rds")

# Recall CellChat objects
cellchat_carotid <- readRDS("fullcarotid_cellchat.rds")
cellchat_femoral <- readRDS("fullfemoral_cellchat.rds")
```

```{r}
# Compute network centrality
cellchat_carotid <- netAnalysis_computeCentrality(cellchat_carotid, slot.name = "netP")
cellchat_femoral <- netAnalysis_computeCentrality(cellchat_femoral, slot.name = "netP")

# Find available genes
available_genes <- cellchat_carotid@data.signaling@Dimnames[[1]]
print(available_genes)

# Signaling network roles
netAnalysis_signalingRole_network(cellchat_carotid, signaling = "MIF", width = 8, height = 4, font.size = 12)
netAnalysis_signalingRole_network(cellchat_carotid, signaling = "SPP1", width = 8, height = 4, font.size = 12)

netAnalysis_signalingRole_network(cellchat_femoral, signaling = "MIF", width = 8, height = 4, font.size = 12)
netAnalysis_signalingRole_network(cellchat_femoral, signaling = "SPP1", width = 8, height = 4, font.size = 12)
```

```{r}
library(NMF)
library(ggalluvial)

# selectK to infer the number of patterns, which is based on two metrics that have been implemented in the NMF R package, including Cophenetic and Silhouette. 
selectK(cellchat_carotid, pattern = "outgoing")
selectK(cellchat_femoral, pattern = "outgoing")

# Both Cophenetic and Silhouette values begin to drop suddenly when the number of outgoing patterns is 3.
nPatterns = 2
# Outgoing patterns reveal how the sender cells (i.e. cells as signal source) coordinate with each other as well as how they coordinate with certain signaling pathways to drive communication.
cellchat_carotid <- identifyCommunicationPatterns(cellchat_carotid, pattern = "outgoing", k = nPatterns)
cellchat_femoral <- identifyCommunicationPatterns(cellchat_femoral, pattern = "outgoing", k = nPatterns)
```

```{r}
# Unsuccessful dataset merging
# object.list <- list(carotid=cellchat_carotid, femoral=cellchat_femoral)
# cellchat <- mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE)
```

```{r}
gg1 <- netVisual_heatmap(cellchat_carotid, color.heatmap = c("#F2E9EA", "#2A9E00"))
gg2 <- netVisual_heatmap(cellchat_carotid, measure = "weight", color.heatmap = c("#F2E9EA", "#2A9E00"), title.name = "Carotid Interaction")

gg3 <- netVisual_heatmap(cellchat_femoral, color.heatmap = c("#F2E9EA", "#2A9E00"))
gg4 <- netVisual_heatmap(cellchat_femoral, measure = "weight", color.heatmap = c("#F2E9EA", "#2A9E00"), title.name = "Femoral Interaction")

gg2 + gg4

gg5 <- netVisual_heatmap(cellchat_carotid, signaling = "MIF", measure = "weight", color.heatmap = c("#F2E9EA", "#2A9E00"), title.name = "Carotid MIF Interaction")
gg6 <- netVisual_heatmap(cellchat_femoral, signaling = "MIF", measure = "weight", color.heatmap = c("#F2E9EA", "#2A9E00"), title.name = "Femoral MIF Interaction")

gg5 + gg6

gg7 <- netVisual_heatmap(cellchat_carotid, signaling = "SPP1", measure = "weight", color.heatmap = c("#F2E9EA", "#2A9E00"), title.name = "Carotid SPP1 Interaction")
gg8 <- netVisual_heatmap(cellchat_femoral, signaling = "SPP1", measure = "weight", color.heatmap = c("#F2E9EA", "#2A9E00"), title.name = "Femoral SPP1 Interaction")

gg7 + gg8
```

```{r}
# cellchat_carotid <- computeNetSimilarityPairwise(cellchat_carotid, type = "functional")
# cellchat_carotid <- netEmbedding(cellchat_carotid, type = "functional",
#                          umap.method = "uwot")
# cellchat_carotid <- netClustering(
#   cellchat_carotid, type = "functional")
# netVisual_embeddingPairwise(
#   cellchat_carotid, type = "functional", label.size = 3.5)
# rankSimilarity(cellchat_carotid, type = "functional")
# rankNet(cellchat_femoral, mode = "comparison", stacked = T, do.stat = TRUE)

# Compute and visualize the pathway distance in the learned joint manifold
# Y <- methods::slot(cellchat_carotid, slot.name)$similarity[[type]]$dr[[comparison.name]]
# rownames(Y) <- rownames(cellchat_carotid@netP[["similarity"]]$functional$matrix$1-2)
# group <- sub(".*--", "", rownames(Y))
# rankSimilarity(cellchat_carotid, slot.name = "netP", type = c("functional","structural"))
# rankNet(cellchat_carotid, mode = "comparison", stacked = T, do.stat = TRUE)
```

```{r}
# Set the working directory 
setwd("/Volumes/biologos/Bioinformatics/Using-CellChat-for-Ligand-Receptor-Network-Analysis/Visuals/Comparison")
# png("BubblePlot.png", width=12, height=5, units="in", res=600)
# bubble <- netVisual_bubble(cellchat, sources.use = c(6,7,9), targets.use = c(1,3,4,8,10), comparison = c(1, 2), remove.isolate = TRUE, angle.x=45, thresh = 0.05, font.size = 7)
# print(bubble)
# dev.off()

# png("BubblePlot-Refined.png", width=7, height=5, units="in", res=3000)
# bubble_refined <- netVisual_bubble(cellchat, sources.use = 7, targets.use = c(1,3,4,8,10), comparison = c(1, 2), remove.isolate = TRUE, angle.x=45, thresh = 0.05, font.size = 5.5)
# print(bubble_refined)
# dev.off()

# Bubble plot
# png("BubblePlot_carotid.png", width=10, height=5, units="in", res=600)
bubble1 <- netVisual_bubble(cellchat_carotid, targets.use = 1, remove.isolate = TRUE, angle.x=45, thresh = 0.05)
print(bubble1)
# dev.off()

# png("BubblePlot_femoral.png", width=10, height=5, units="in", res=600)
bubble2 <- netVisual_bubble(cellchat_femoral, targets.use = 1, remove.isolate = FALSE, angle.x=45, thresh = 0.05)
print(bubble2)
# dev.off()

bubble1+bubble2
```

```{r, echo=F}
## DO NOT DELETE THIS BLOCK!
Sys.info()
sessionInfo()
```