●Code
################
#1. Integration
################
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
library(hdf5r)
library(sctransform)


# 10x data install
WTgdata <- Read10X_h5("filtered_feature_bc_matrix.Wtg_CD45plus.h5")
ELdata <- Read10X_h5("filtered_feature_bc_matrix.EL_CD45plus.h5")

# to Seurat Object
WTg <- CreateSeuratObject(counts = WTgdata, project = "WTg", min.cells = 1, min.features = 200)
EL <- CreateSeuratObject(counts = ELdata, project = "EL", min.cells = 1, min.features = 200)

# Compute percent mito ratio
WTg$mitoRatio <- PercentageFeatureSet(object = WTg, pattern = "^mt-")
WTg$mitoRatio <- WTg@meta.data$mitoRatio / 100
EL$mitoRatio <- PercentageFeatureSet(object = EL, pattern = "^mt-")
EL$mitoRatio <- EL@meta.data$mitoRatio / 100

# Add number of genes per UMI for each cell to metadata
WTg$log10GenesPerUMI <- log10(WTg$nFeature_RNA) / log10(WTg$nCount_RNA)
EL$log10GenesPerUMI <- log10(EL$nFeature_RNA) / log10(EL$nCount_RNA)

# Filter out low quality cells using selected thresholds - these will change with experiment
WTg <- subset(x = WTg, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & mitoRatio < 0.15)
EL <- subset(x = EL, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & mitoRatio < 0.15)

# SCTranform
WTg <- SCTransform(WTg, vars.to.regress = c("mitoRatio"))
EL <- SCTransform(EL, vars.to.regress = c("mitoRatio"))

# Create a merged Seurat object
merged_seurat <- merge(x = WTg, y = EL, add.cell.id = c("WTg", "EL"))
View(merged_seurat@meta.data)

# Create metadata dataframe
metadata <- merged_seurat@meta.data

# Add cell IDs to metadata
metadata$cells <- rownames(metadata)

# Create sample column
metadata$sample <- NA
metadata$sample[which(str_detect(metadata$cells, "^WTg_"))] <- "WTg"
metadata$sample[which(str_detect(metadata$cells, "^EL_"))] <- "EL"

# Add metadata back to Seurat object
merged_seurat@meta.data <- metadata

## FigS6A
# Visualize the number of cell counts per sample
metadata %>% 
  	ggplot(aes(x=sample, fill=sample)) + 
  	geom_bar() +
  	theme_classic() +
  	theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  	theme(plot.title = element_text(hjust=0.5, face="bold")) +
  	ggtitle("NCells")

# Visualize the number UMIs/transcripts per cell
metadata %>% 
  	ggplot(aes(color=sample, x=nCount_RNA, fill= sample)) + 
  	geom_density(alpha = 0.2) + 
  	scale_x_log10() + 
  	theme_classic() +
  	ylab("Cell density") +
  	geom_vline(xintercept = 500)

# Visualize the distribution of genes detected per cell via histogram
metadata %>% 
  	ggplot(aes(color=sample, x=nFeature_RNA, fill= sample)) + 
  	geom_density(alpha = 0.2) + 
  	theme_classic() +
  	scale_x_log10() + 
  	geom_vline(xintercept = 300)

# Visualize the overall complexity of the gene expression by visualizing the genes detected per UMI (novelty score)
metadata %>%
  	ggplot(aes(x=log10GenesPerUMI, color = sample, fill=sample)) +
  	geom_density(alpha = 0.2) +
  	theme_classic() +
  	geom_vline(xintercept = 0.8)

# Visualize the distribution of mitochondrial gene expression detected per cell
metadata %>% 
  	ggplot(aes(color=sample, x=mitoRatio, fill=sample)) + 
  	geom_density(alpha = 0.2) + 
  	scale_x_log10() + 
  	theme_classic() +
  	geom_vline(xintercept = 0.2)

# Visualize the correlation between genes detected and number of UMIs and determine whether strong presence of cells with low numbers of genes/UMIs
metadata %>% 
  	ggplot(aes(x=nCount_RNA, y=nFeature_RNA, color=mitoRatio)) + 
  	geom_point() + 
	scale_colour_gradient(low = "gray90", high = "black") +
  	stat_smooth(method=lm) +
  	scale_x_log10() + 
  	scale_y_log10() + 
  	theme_classic() +
  	geom_vline(xintercept = 500) +
  	geom_hline(yintercept = 250) +
  	facet_wrap(~sample)


# Split seurat object by condition to perform SCT on all samples
split_seurat <- SplitObject(merged_seurat, split.by = "sample")
split_seurat <- split_seurat[c("WTg", "EL")]

# perform SCTransform normalization
split_seurat <- lapply(X = split_seurat, FUN = SCTransform)

## Integration using CCA
integ_features <- SelectIntegrationFeatures(object.list = split_seurat,nfeatures = 10000) 

# Prepare the SCT list object for integration
split_seurat <- PrepSCTIntegration(object.list = split_seurat,anchor.features = integ_features)

# Find best buddies - can take a while to run
integ_anchors <- FindIntegrationAnchors(object.list = split_seurat, normalization.method = "SCT", anchor.features = integ_features)

# Integrate across conditions
seurat_integrated <- IntegrateData(anchorset = integ_anchors, normalization.method = "SCT")

# Run PCA
seurat_integrated <- RunPCA(object = seurat_integrated, npcs = 30, verbose = FALSE)

# Plot PCA
PCAPlot(seurat_integrated,split.by = "sample")  

# Run UMAP
seurat_integrated <- RunUMAP(seurat_integrated,dims = 1:30,reduction = "pca")

# Plot UMAP                             
DimPlot(seurat_integrated)

# Plot UMAP split by sample
DimPlot(seurat_integrated,split.by = "sample")  

# Save integrated seurat object
saveRDS(seurat_integrated, "seurat_integrated.rds")


## Explore heatmap of PCs
DimHeatmap(seurat_integrated, dims = 1:9, cells = 500, balanced = TRUE)

# Printing out the most variable genes driving PCs
print(x = seurat_integrated[["pca"]], dims = 1:10, nfeatures = 5)

# Plot the elbow plot
ElbowPlot(object = seurat_integrated, ndims = 30)

# Determine the K-nearest neighbor graph
seurat_integrated <- FindNeighbors(object = seurat_integrated,dims = 1:30)

# Determine the clusters for various resolutions                                
seurat_integrated <- FindClusters(object = seurat_integrated,resolution = c(0.5, 0.75, 1.0))

# Assign identity of clusters
Idents(object = seurat_integrated) <- "integrated_snn_res.0.5"

# Plot the UMAP
DimPlot(seurat_integrated,
        reduction = "umap",
        label = TRUE,
        label.size = 6)

Idents(object = seurat_integrated) <- "integrated_snn_res.0.75"
DimPlot(seurat_integrated,
        reduction = "umap",
        label = TRUE,
        label.size = 6)

## Subsequent analyses were performed at "0.5"

# FigS6B
# UMAP of cells in each cluster by sample
DimPlot(seurat_integrated, 
        label = TRUE, 
        split.by = "sample")  + NoLegend()


### for Fig5C ###
## Extract identity and sample information from seurat object to determine the number of cells per cluster per sample
n_cells <- FetchData(seurat_integrated, 
                     vars = c("ident", "orig.ident")) %>%
        dplyr::count(ident, orig.ident) %>%
        tidyr::spread(ident, n)

# View table
View(n_cells)

# UMAP of cells in each cluster by sample
DimPlot(seurat_integrated, 
        label = TRUE, 
        split.by = "sample")  + NoLegend()

saveRDS(seurat_integrated, "seurat_integrated_0.5.rds")


##### heatmap ######
table(seurat_integrated@meta.data[,c("orig.ident","integrated_snn_res.0.5")])

temp.mat <- t(table(seurat_integrated@meta.data[,c("orig.ident","integrated_snn_res.0.5")]))
temp.mat.scale <- temp.mat/apply(temp.mat, 2 , sum)
temp.mat.scale2 <- temp.mat.scale/apply(temp.mat.scale, 1, max)

library(gplots)
heatmap.2(temp.mat.scale2, col=colorRampPalette(c("#0068b7","white","#f39800")), trace="none", density.info="none", Colv=F, cexCol=0.8, cexRow=0.8, keysize=1.2, margins=c(6,14))


######################
#2. Find marker genes
######################
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
library(AnnotationHub)
library(ensembldb)
library(multtest)
library(metap)

DefaultAssay(seurat_integrated) <- "SCT"


## https://doi.org/10.1172/jci.insight.129212.
FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("Amy1", "Amy2a2", "Pyy", "Sst", "Col1a1", "Col1a2", "Ins1", "Ins2", "Cdh5", "Pecam1", "Sox9", "Krt18"),
            label = TRUE, 
            order = TRUE,
            min.cutoff = 'q10',
	    repel = TRUE
            )

FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("Cd14", "Adgre1", "Cd52", "Cd3d", "Cd3e", "Cd19", "Cd79b", "Hba-a1", "Hba-a2"),
            label = TRUE, 
            order = TRUE,
            min.cutoff = 'q10',
	    repel = TRUE
            )

# FigS6C
markers.to.plot <- c("Amy1", "Amy2a2", "Pyy", "Sst", "Col1a1", "Col1a2", "Ins1", "Ins2", "Cdh5", "Pecam1", "Sox9", "Krt18", "Cd14", "Adgre1", "Cd52", "Cd3d", "Cd3e", "Cd19", "Cd79b", "Hba-a1", "Hba-a2")
VlnPlot(seurat_integrated, 
        features = markers.to.plot, stack = T, flip = T)


# Fig3B
## Plot genes highly expressed in ACC obtained from GSE48643
FeaturePlot(object = seurat_integrated, 
                        features = c("Cyp3a13", "Muc1", "Apob"),
                         order = TRUE, 
                         label = TRUE,
			 repel = TRUE, cols = c("lightgray","darkred"), ncol = 3)



seurat_rename[["cell.label"]] <- NA
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "0"] <- "ACC"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "1"] <- "ACC"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "2"] <- "Macrophages"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "3"] <- "ACC"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "4"] <- "ACC"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "5"] <- "Monocytes"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "6"] <- "T/NK cells"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "7"] <- "Duct cells"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "8"] <- "ACC"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "9"] <- "Macrophages"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "10"] <- "B cells"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "11"] <- "Macrophages"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "12"] <- "Endothelial cells"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "13"] <- "ACC"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "14"] <- "ACC"
seurat_rename$cell.label[seurat_rename$integrated_snn_res.0.5 == "15"] <- "Fibroblasts"

new.cluster.ids <- c("ACC",
                     "ACC",
                     "Macrophages",
                     "ACC",
                     "ACC",
                     "Monocytes",
                     "T/NK cells",
                     "Duct cells",
                     "ACC",
                     "Macrophages",
                     "B cells",
                     "Macrophages",
                     "Endothelial cells",
                     "ACC",
			   "ACC",
			   "Fibroblasts")
names(new.cluster.ids) <- levels(seurat_rename)
seurat_rename<- RenameIdents(seurat_rename, new.cluster.ids)


# Fig3A
# Plot the UMAP
DimPlot(object = seurat_rename, 
        reduction = "umap", 
        label = TRUE,
        label.size = 3,
        repel = TRUE)

p1 <- DimPlot(seurat_rename, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(seurat_rename, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

saveRDS(seurat_rename, "seurat_renameR.rds")


###########
#3. ssGSEA
###########
library(escape)
library(dittoSeq)
library(SingleCellExperiment)
library(Seurat)
library(SeuratObject)
library(UCell)
library(tidyverse)

seurat_ACC <- subset(seurat_renameR, idents = "ACC")

# Getting Gene Sets
GS.hallmark <- getGeneSets(library = "H", species = "Mus musculus")

# Loading Processed Single-Cell Data
# Enrichment
ES.seurat <- enrichIt(obj = seurat_ACC, 
                      gene.sets = GS.hallmark, 
                      groups = 1000, cores = 2, 
                      min.size = 5)

seurat_ACC <- Seurat::AddMetaData(seurat_ACC, ES.seurat)

# Visualizations
colors <- colorRampPalette(c("#0D0887FF","#7E03A8FF","#CC4678FF","#F89441FF","#F0F921FF"))

## Heatmap
# Fig3C
dittoHeatmap(seurat_ACC, genes = NULL, metas = names(ES.seurat), 
             annot.by = "sample", 
             fontsize = 7, 
             cluster_cols = T,
             heatmap.colors = colors(50))

#FigS6D
dittoHeatmap(seurat_ACC, genes = NULL, 
             metas = c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", 
                       "HALLMARK_ANGIOGENESIS", 
                       "HALLMARK_HYPOXIA", 
                       "HALLMARK_MITOTIC_SPINDLE", 
                       "HALLMARK_E2F_TARGETS",
                       "HALLMARK_TGF_BETA_SIGNALING"), 
             annot.by = "sample", 
             fontsize = 7,
             heatmap.colors = colors(50))

# Fig3D
## Violin plots
multi_dittoPlot(seurat_ACC, vars = c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", 
                       "HALLMARK_ANGIOGENESIS", 
                       "HALLMARK_HYPOXIA", 
                       "HALLMARK_MITOTIC_SPINDLE", 
                       "HALLMARK_E2F_TARGETS",
                       "HALLMARK_TGF_BETA_SIGNALING"), 
                group.by = "sample", plots = c("jitter", "vlnplot", "boxplot"), 
                ylab = "Enrichment Scores", 
                theme = theme_classic() + theme(plot.title = element_text(size = 10)))


#############
#4. CellChat 
#############
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
library(hdf5r)
library(sctransform)
library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
library(ggalluvial)
library(ComplexHeatmap)
library(wordcloud)

### cellchat ###
cellchat <- createCellChat(object = seurat_renameR, group.by = "cell.label", assay = "SCT")
CellChatDB <- CellChatDB.mouse
CellChatDB.use <- CellChatDB
cellchat@DB <- CellChatDB.use

cellchat <- subsetData(cellchat)
future::plan("multisession", workers = 12)
options(future.globals.maxSize = 1000 * 1024^2)

cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
cellchat <- projectData(cellchat, PPI.mouse)

cellchat <- computeCommunProb(cellchat, raw.use = TRUE)
cellchat <- filterCommunication(cellchat, min.cells = 10)

cellchat <- computeCommunProbPathway(cellchat)

cellchat <- aggregateNet(cellchat)

groupSize <- as.numeric(table(cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)

# circular plot
netVisual_circle(cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
# Fig3E
netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")

# Fig.3F
netVisual_chord_gene(cellchat, 
                     sources.use = c(1,4,5,6), 
                     targets.use = c(1,4,5,6),  
                     color.use = brewer.pal(10, "Paired"), reduce = 0.003)


## Identify signaling roles (e.g., dominant senders, receivers) of cell groups as well as the major contributing signaling
# Compute the network centrality scores
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")

# Fig.3G
# Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways
# Heatmap incoming outgoing signals
ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing", font.size = 5, height = 15)
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming", font.size = 5, height = 15)
ht1 + ht2


# FigS7A
netAnalysis_signalingRole_network(cellchat, signaling = "COLLAGEN")
netAnalysis_signalingRole_network(cellchat, signaling = "THBS")
netAnalysis_signalingRole_network(cellchat, signaling = "LAMININ")
netAnalysis_signalingRole_network(cellchat, signaling = "FN1")
netAnalysis_signalingRole_network(cellchat, signaling = "MIF")
netAnalysis_signalingRole_network(cellchat, signaling = "APP")
netAnalysis_signalingRole_network(cellchat, signaling = "SPP1")
netAnalysis_signalingRole_network(cellchat, signaling = "PERIOSTIN")
netAnalysis_signalingRole_network(cellchat, signaling = "HSPG")
netAnalysis_signalingRole_network(cellchat, signaling = "TWEAK")
netAnalysis_signalingRole_network(cellchat, signaling = "MPZ")
netAnalysis_signalingRole_network(cellchat, signaling = "VISFATIN")
netAnalysis_signalingRole_network(cellchat, signaling = "IGF")
netAnalysis_signalingRole_network(cellchat, signaling = "PARs")
netAnalysis_signalingRole_network(cellchat, signaling = "JAM")
netAnalysis_signalingRole_network(cellchat, signaling = "NOTCH")

saveRDS(cellchat, file = "cellchat_renameR.rds")

####################
# isolate WTg and EL
#################### 
Idents(object = seurat_renameR) <- "orig.ident"
levels(seurat_renameR)

WTg <- subset(seurat_renameR, idents = "WTg")
WTg

EL <- subset(seurat_renameR, idents = "EL")
EL

cellchat.WTg <- createCellChat(object = WTg, group.by = "cell.label", assay = "SCT")
cellchat.EL <- createCellChat(object = EL, group.by = "cell.label", assay = "SCT")

cellchat.WTg <- updateCellChat(cellchat.WTg)
cellchat.EL <- updateCellChat(cellchat.EL)

### WTg ###
cellchat.WTg@DB <- CellChatDB.use
cellchat.WTg <- subsetData(cellchat.WTg)
future::plan("multisession", workers = 12) # do parallel

# Compute the communication probability and infer cellular communication network
cellchat.WTg <- identifyOverExpressedGenes(cellchat.WTg)
cellchat.WTg <- identifyOverExpressedInteractions(cellchat.WTg)
cellchat.WTg <- projectData(cellchat.WTg, PPI.mouse)

## Inference of cell-cell communication network
cellchat.WTg <- computeCommunProb(cellchat.WTg, raw.use = TRUE)
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat.WTg <- filterCommunication(cellchat.WTg, min.cells = 10)

cellchat.WTg <- computeCommunProbPathway(cellchat.WTg)

cellchat.WTg <- aggregateNet(cellchat.WTg)

groupSize <- as.numeric(table(cellchat.WTg@idents))
par(mfrow = c(1,2), xpd=TRUE)

### EL ###
cellchat.EL@DB <- CellChatDB.use
cellchat.EL <- subsetData(cellchat.EL)
future::plan("multisession", workers = 12) # do parallel

# Compute the communication probability and infer cellular communication network
cellchat.EL <- identifyOverExpressedGenes(cellchat.EL)
cellchat.EL <- identifyOverExpressedInteractions(cellchat.EL)
cellchat.EL <- projectData(cellchat.EL, PPI.mouse)

## Inference of cell-cell communication network
cellchat.EL <- computeCommunProb(cellchat.EL, raw.use = TRUE)
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat.EL <- filterCommunication(cellchat.EL, min.cells = 10)

cellchat.EL <- computeCommunProbPathway(cellchat.EL)

cellchat.EL <- aggregateNet(cellchat.EL)

groupSize <- as.numeric(table(cellchat.EL@idents))
par(mfrow = c(1,2), xpd=TRUE)


### merge cellchat objects ###
object.list <- list(EL = cellchat.EL, WTg = cellchat.WTg)
cellchat <- mergeCellChat(object.list, add.names = names(object.list))


saveRDS(cellchat, file = "cellchat_renameR_WTg_EL.rds")
saveRDS(cellchat.WTg, file = "cellchat_renameR_WTg.rds")
saveRDS(cellchat.EL, file = "cellchat_renameR_EL.rds")


# Compare the total number of interactions and interaction strength
gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2))
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "weight")
gg1 + gg2


# Differential number of interactions or interaction strength among different cell populations
# circularplot
par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")


# FigS7B
# circularplot2
weight.max <- getMaxWeight(object.list, attribute = c("idents","weight"))
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_circle(object.list[[i]]@net$weight, weight.scale = T, label.edge= F, edge.weight.max = weight.max[2], edge.width.max = 12, title.name = paste0("Interaction weights/strength - ", names(object.list)[i]))
}

# FigS7C
## Compare the major sources and targets in 2D space
cellchat.WTg <- netAnalysis_computeCentrality(cellchat.WTg, slot.name = "netP")
cellchat.EL <- netAnalysis_computeCentrality(cellchat.EL, slot.name = "netP")
object.list <- list(EL = cellchat.EL, WTg = cellchat.WTg)
cellchat <- mergeCellChat(object.list, add.names = names(object.list))

num.link <- sapply(object.list, function(x) {rowSums(x@net$count) + colSums(x@net$count)-diag(x@net$count)})
weight.MinMax <- c(min(num.link), max(num.link)) # control the dot size in the different datasets
gg <- list()
for (i in 1:length(object.list)) {
  gg[[i]] <- netAnalysis_signalingRole_scatter(object.list[[i]], title = names(object.list)[i], weight.MinMax = weight.MinMax)
}
patchwork::wrap_plots(plots = gg)

# Compare the overall information flow of each signaling pathway
gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
gg1 + gg2

# FigS7D
# Compare outgoing (or incoming) signaling associated with each cell population
library(ComplexHeatmap)
cellchat.WTg <- netAnalysis_computeCentrality(cellchat.WTg, slot.name = "netP")
cellchat.EL <- netAnalysis_computeCentrality(cellchat.EL, slot.name = "netP")
object.list <- list(EL = cellchat.EL, WTg = cellchat.WTg)

i = 1  
pathway.union <- union(object.list[[i]]@netP$pathways, object.list[[i+1]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 15)
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i+1], width = 5, height = 15)

ht3 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 15, color.heatmap = "GnBu")
ht4 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i+1], width = 5, height = 15, color.heatmap = "GnBu")

draw(ht1 + ht2 + ht3 + ht4, ht_gap = unit(0.5, "cm"))


########################################
#5. re-clustering of TME 
# Only this section was reanalyzed in December 2024, so Seurat V5 was used. Therefore, the integration method is different from the others.
########################################
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
library(hdf5r)
library(sctransform)
library(viridis)

options(Seurat.object.assay.version = "v3")

WTgdata <- Read10X_h5("filtered_feature_bc_matrix.Wtg_CD45plus.h5")
ELdata <- Read10X_h5("filtered_feature_bc_matrix.EL_CD45plus.h5")

TME <- rownames(seurat_renameR@meta.data[which(seurat_renameR@meta.data[,"cell.label"] %in% c("Endothelial cells","Fibroblasts", "Macrophages")),])

TME.el <- gsub("EL_", "",  TME[grep("EL_", TME)])
TME.wtg <- gsub("WTg_", "",  TME[grep("WTg_", TME)])

WTg <- CreateSeuratObject(counts = WTgdata[,TME.wtg], project = "WTg", min.cells = 1, min.features = 200)
EL <- CreateSeuratObject(counts = ELdata[,TME.el], project = "EL", min.cells = 1, min.features = 200)

# Compute percent mito ratio
WTg$mitoRatio <- PercentageFeatureSet(object = WTg, pattern = "^mt-")
WTg$mitoRatio <- WTg@meta.data$mitoRatio / 100
EL$mitoRatio <- PercentageFeatureSet(object = EL, pattern = "^mt-")
EL$mitoRatio <- EL@meta.data$mitoRatio / 100

# Add number of genes per UMI for each cell to metadata
WTg$log10GenesPerUMI <- log10(WTg$nFeature_RNA) / log10(WTg$nCount_RNA)
EL$log10GenesPerUMI <- log10(EL$nFeature_RNA) / log10(EL$nCount_RNA)

# Filter out low quality cells using selected thresholds - these will change with experiment
WTg <- subset(x = WTg, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & mitoRatio < 0.15)
EL <- subset(x = EL, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & mitoRatio < 0.15)

# Create a merged Seurat object
merged_seurat <- merge(x = WTg, y = EL, add.cell.id = c("WTg", "EL"))
View(merged_seurat@meta.data)

# Create metadata dataframe
metadata <- merged_seurat@meta.data

# Add cell IDs to metadata
metadata$cells <- rownames(metadata)

# Create sample column
metadata$sample <- NA
metadata$sample[which(str_detect(metadata$cells, "^WTg_"))] <- "WTg"
metadata$sample[which(str_detect(metadata$cells, "^EL_"))] <- "EL"


# Add metadata back to Seurat object
merged_seurat@meta.data <- metadata


# Split seurat object by condition to perform SCT on all samples
split_seurat <- SplitObject(merged_seurat, split.by = "sample")
split_seurat <- split_seurat[c("WTg", "EL")]

# perform SCTransform normalization
split_seurat <- lapply(X = split_seurat, FUN = SCTransform)

## Integration using CCA
integ_features <- SelectIntegrationFeatures(object.list = split_seurat,nfeatures = 10000) 

# Prepare the SCT list object for integration
split_seurat <- PrepSCTIntegration(object.list = split_seurat,anchor.features = integ_features)

# Find best buddies - can take a while to run
integ_anchors <- FindIntegrationAnchors(object.list = split_seurat, normalization.method = "SCT", anchor.features = integ_features)

# Integrate across conditions
seurat_integrated <- IntegrateData(anchorset = integ_anchors, normalization.method = "SCT")

# Run PCA
seurat_integrated <- RunPCA(object = seurat_integrated, npcs = 30, verbose = FALSE)

# Plot PCA
PCAPlot(seurat_integrated,split.by = "sample")  

# Run UMAP
seurat_integrated <- RunUMAP(seurat_integrated,dims = 1:30,reduction = "pca")

# Determine the K-nearest neighbor graph
seurat_integrated <- FindNeighbors(object = seurat_integrated,dims = 1:30)

# Determine the clusters for various resolutions                                
seurat_integrated <- FindClusters(object = seurat_integrated,resolution = c(0.5, 0.75, 1.0))

# Assign identity of clusters
Idents(object = seurat_integrated) <- "integrated_snn_res.0.5"


# Fig.S8A
DimPlot(object = seurat_integrated, 
        reduction = "umap", 
        label = TRUE,
        label.size = 3,
        repel = TRUE)

# Fig.S8B
markers.to.plot <- c("Col1a1", "Col1a2", "Cdh5", "Pecam1", "Cd14", "Adgre1", "Cd52")
VlnPlot(seurat_integrated, 
        features = markers.to.plot, stack = T, flip = T)

saveRDS(seurat_integrated, "seurat_integrated_TME.rds")

## TAM marker
# Cell 2021 Vol. 184 Issue 11 Pages 2988-3005 e16
TAM <- c("Cd14", "Fcgr3", "Adgre1", "Apoe", "Trem2", "C1qa")

## M1/M2 macrophage marker
## JCI Insight 2019:4:e126556
M1  <- c("Azin1", "Cd38", "Cd86", "Cxcl10", "Cxcl9", "Fpr2", "Gpr18", "Il12b", "Il18", "Irf5", "Nfkbiz", "Nos2", "Ptgs2", "Tlr4", "Tnf")
M2 <- c("Alox15", "Arg1", "Chil3", "Chil4", "Egr2", "Il10", "Irf4", "Klf4", "Myc", "Socs2", "Tgm2")

## CAF marker
## PLoS One 9:e91808. doi: 10.1371/journal.pone.0091808
iCAF <- c("Col1a1", "Vim", "Dcn", "Pdgfra", "Cxcl12")
myCAF <- c("Col1a1", "Vim","Tagln", "Acta2", "Fap")

## Tumor Endothelial cells (TEC)
## Front. Immunol., 08 February 2019, Sec. Cytokines and Soluble Mediators in Immunity, Volume 10 - 2019 | https://doi.org/10.3389/fimmu.2019.00147
TEC <- c("Cdh5", "Pecam1","Itgb3","Cxcr4")


seurat_integrated <- AddModuleScore(seurat_integrated, features = list(TAM), name = "TAM_score")
seurat_integrated <- AddModuleScore(seurat_integrated, features = list(M1), name = "M1_score")
seurat_integrated <- AddModuleScore(seurat_integrated, features = list(M2), name = "M2_score")
seurat_integrated <- AddModuleScore(seurat_integrated, features = list(iCAF), name = "iCAF_score")
seurat_integrated <- AddModuleScore(seurat_integrated, features = list(myCAF), name = "myCAF_score")
seurat_integrated <- AddModuleScore(seurat_integrated, features = list(TEC), name = "TEC_score")


# Fig4A
FeaturePlot(seurat_integrated, features = c("TEC_score1", "TAM_score1", "M1_score1", "M2_score1","iCAF_score1","myCAF_score1"), order = T, pt.size = 2) & scale_color_viridis(option = "H")

# Fig.S7C
P1 <- VlnPlot(seurat_integrated, features = c("TAM_score1"), flip = T) & geom_hline(yintercept = median(seurat_integrated@meta.data$TAM_score1))
P2 <- VlnPlot(seurat_integrated, features = c("M1_score1"), flip = T) & geom_hline(yintercept = median(seurat_integrated@meta.data$M1_score1))
P3 <- VlnPlot(seurat_integrated, features = c("M2_score1"), flip = T) & geom_hline(yintercept = median(seurat_integrated@meta.data$M2_score1))
P4 <- VlnPlot(seurat_integrated, features = c("iCAF_score1"), flip = T) & geom_hline(yintercept = median(seurat_integrated@meta.data$iCAF_score1))
P5 <- VlnPlot(seurat_integrated, features = c("myCAF_score1"), flip = T) & geom_hline(yintercept = median(seurat_integrated@meta.data$myCAF_score1))
P6 <- VlnPlot(seurat_integrated, features = c("TEC_score1"), flip = T) & geom_hline(yintercept = median(seurat_integrated@meta.data$TEC_score1))
P1 + P2 + P3 + P4 + P5 + P6 + NoLegend()


seurat_integrated[["cell.label"]] <- NA
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "0"] <- "TAM"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "1"] <- "Macrophages"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "2"] <- "TAM"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "3"] <- "TAM"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "4"] <- "Macrophages"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "5"] <- "TEC"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "6"] <- "Macrophages"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "7"] <- "TAM"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "8"] <- "Macrophages"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "9"] <- "Macrophages"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "10"] <- "myCAF"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "11"] <- "TEC"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "12"] <- "iCAF"
seurat_integrated$cell.label[seurat_integrated$integrated_snn_res.0.5 == "13"] <- "Macrophages"

Idents(object = seurat_integrated) <- "cell.label"

# Fig.S8D
# Plot the UMAP
DimPlot(object = seurat_integrated, 
        reduction = "umap", 
        label = TRUE,
        label.size = 3,
        repel = TRUE)

saveRDS(seurat_integrated, "seurat_TME.rds")


#############
#6. CellChat for re-clustaering 
#############
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
library(hdf5r)
library(sctransform)
library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
library(ggalluvial)
library(ComplexHeatmap)
library(wordcloud)


seurat_ACC <- subset(seurat_renameR, idents = "ACC")
seurat_TME_subset <- subset(seurat_TME, idents = c("iCAF", "myCAF", "TEC", "TAM"))

ACC_TME <- merge(seurat_ACC, seurat_TME_subset)
levels(ACC_TME)
saveRDS(ACC_TME, "seurat_ACC_TME.rds")

####################
##### CellChat #####
####################
cellchat <- createCellChat(object = ACC_TME, group.by = "cell.label", assay = "SCT")
CellChatDB <- CellChatDB.mouse

CellChatDB.use <- CellChatDB
cellchat@DB <- CellChatDB.use

cellchat <- subsetData(cellchat)
future::plan("multisession", workers = 12)

cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
cellchat <- projectData(cellchat, PPI.mouse)

cellchat <- computeCommunProb(cellchat, raw.use = TRUE)
cellchat <- filterCommunication(cellchat, min.cells = 10)

cellchat <- computeCommunProbPathway(cellchat)

cellchat <- aggregateNet(cellchat)

groupSize <- as.numeric(table(cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)

# Fig4B
# circular plot
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")

netVisual_chord_gene(cellchat, 
                     sources.use = c(1), 
                     targets.use = c(1:5),  
                     color.use = brewer.pal(10, "Paired"), reduce = 0.005)

## Identify signaling roles (e.g., dominant senders, receivers) of cell groups as well as the major contributing signaling
# Compute the network centrality scores
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")

# Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways
# Heatmap incoming outgoing signals
ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing", font.size = 5, height = 15)
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming", font.size = 5, height = 15)
ht1 + ht2


saveRDS(cellchat, file = "cellchat_ACC_TME.rds")



# isolate WTg and EL 
Idents(object = ACC_TME) <- "orig.ident"
levels(ACC_TME)

WTg <- subset(ACC_TME, idents = "WTg")
WTg

EL <- subset(ACC_TME, idents = "EL")
EL

cellchat.WTg <- createCellChat(object = WTg, group.by = "cell.label", assay = "SCT")
cellchat.EL <- createCellChat(object = EL, group.by = "cell.label", assay = "SCT")

cellchat.WTg <- updateCellChat(cellchat.WTg)
cellchat.EL <- updateCellChat(cellchat.EL)

### WTg ###
cellchat.WTg@DB <- CellChatDB.use
cellchat.WTg <- subsetData(cellchat.WTg)
future::plan("multisession", workers = 12) # do parallel

# Compute the communication probability and infer cellular communication network
cellchat.WTg <- identifyOverExpressedGenes(cellchat.WTg)
cellchat.WTg <- identifyOverExpressedInteractions(cellchat.WTg)
cellchat.WTg <- projectData(cellchat.WTg, PPI.mouse)

## Inference of cell-cell communication network
cellchat.WTg <- computeCommunProb(cellchat.WTg, raw.use = TRUE)
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat.WTg <- filterCommunication(cellchat.WTg, min.cells = 10)

cellchat.WTg <- computeCommunProbPathway(cellchat.WTg)

cellchat.WTg <- aggregateNet(cellchat.WTg)

groupSize <- as.numeric(table(cellchat.WTg@idents))
par(mfrow = c(1,2), xpd=TRUE)

### EL ###
cellchat.EL@DB <- CellChatDB.use
cellchat.EL <- subsetData(cellchat.EL)
future::plan("multisession", workers = 12) # do parallel

# Compute the communication probability and infer cellular communication network
cellchat.EL <- identifyOverExpressedGenes(cellchat.EL)
cellchat.EL <- identifyOverExpressedInteractions(cellchat.EL)
cellchat.EL <- projectData(cellchat.EL, PPI.mouse)

## Inference of cell-cell communication network
cellchat.EL <- computeCommunProb(cellchat.EL, raw.use = TRUE)
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat.EL <- filterCommunication(cellchat.EL, min.cells = 10)

cellchat.EL <- computeCommunProbPathway(cellchat.EL)

cellchat.EL <- aggregateNet(cellchat.EL)

groupSize <- as.numeric(table(cellchat.EL@idents))
par(mfrow = c(1,2), xpd=TRUE)


### merge cellchat objects ###
object.list <- list(EL = cellchat.EL, WTg = cellchat.WTg)
cellchat <- mergeCellChat(object.list, add.names = names(object.list))


saveRDS(cellchat, file = "cellchat_ACC_TME_WTg_EL.rds")
saveRDS(cellchat.WTg, file = "cellchat_ACC_TME_WTg.rds")
saveRDS(cellchat.EL, file = "cellchat_ACC_TME_EL.rds")



# Compare the total number of interactions and interaction strength
gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2))
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "weight")
gg1 + gg2


# Differential number of interactions or interaction strength among different cell populations
# circularplot
par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")

#Fig.S9A
# circularplot2
weight.max <- getMaxWeight(object.list, attribute = c("idents","weight"))
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_circle(object.list[[i]]@net$weight, weight.scale = T, label.edge= F, edge.weight.max = weight.max[2], edge.width.max = 12, title.name = paste0("Interaction weights/strength - ", names(object.list)[i]))
}

## Compare the major sources and targets in 2D space
cellchat.WTg <- netAnalysis_computeCentrality(cellchat.WTg, slot.name = "netP")
cellchat.EL <- netAnalysis_computeCentrality(cellchat.EL, slot.name = "netP")
object.list <- list(EL = cellchat.EL, WTg = cellchat.WTg)
cellchat <- mergeCellChat(object.list, add.names = names(object.list))

num.link <- sapply(object.list, function(x) {rowSums(x@net$count) + colSums(x@net$count)-diag(x@net$count)})
weight.MinMax <- c(min(num.link), max(num.link)) # control the dot size in the different datasets
gg <- list()
for (i in 1:length(object.list)) {
  gg[[i]] <- netAnalysis_signalingRole_scatter(object.list[[i]], title = names(object.list)[i], weight.MinMax = weight.MinMax)
+ scale_x_continuous(limits = c(0, 10)) + scale_y_continuous(limits = c(0, 7.5))
}
patchwork::wrap_plots(plots = gg)

gg1 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "ACC")
gg2 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "iCAF")
gg3 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "myCAF")
gg4 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "TAM")
gg5 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "TEC")
patchwork::wrap_plots(plots = list(gg1,gg2,gg3,gg4,gg5))


# Compare the overall information flow of each signaling pathway
gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
gg1 + gg2


#Fig.S9B
# Compare outgoing (or incoming) signaling associated with each cell population
library(ComplexHeatmap)
cellchat.WTg <- netAnalysis_computeCentrality(cellchat.WTg, slot.name = "netP")
cellchat.EL <- netAnalysis_computeCentrality(cellchat.EL, slot.name = "netP")
object.list <- list(EL = cellchat.EL, WTg = cellchat.WTg)

i = 1
pathway.union <- union(object.list[[i]]@netP$pathways, object.list[[i+1]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 15)
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i+1], width = 5, height = 15)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 15, color.heatmap = "GnBu")
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i+1], width = 5, height = 15, color.heatmap = "GnBu")
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "all", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 15, color.heatmap = "OrRd")
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "all", signaling = pathway.union, title = names(object.list)[i+1], width = 5, height = 15, color.heatmap = "OrRd")
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))


i = 1
pathway.union <- union(object.list[[i]]@netP$pathways, object.list[[i+1]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 15)
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i+1], width = 5, height = 15)

ht3 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 15, color.heatmap = "GnBu")
ht4 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i+1], width = 5, height = 15, color.heatmap = "GnBu")

draw(ht1 + ht2 + ht3 + ht4, ht_gap = unit(0.5, "cm"))


## 1=ACC, 2=iCAF, 3=myCAF, 4=TAM, 5=TEC

## Identify dysfunctional signaling by using differential expression analysis
# define a positive dataset, i.e., the dataset with positive fold change against the other dataset
object.list <- list(EL = cellchat.EL, WTg = cellchat.WTg)
pos.dataset = "WTg"
# define a char name used for storing the results of differential expression analysis
features.name = pos.dataset
# perform differential expression analysis
cellchat <- identifyOverExpressedGenes(cellchat, group.dataset = "datasets", pos.dataset = pos.dataset, features.name = features.name, only.pos = FALSE, thresh.pc = 0.1, thresh.fc = 0.1, thresh.p = 1)
#> Use the joint cell labels from the merged CellChat object
# map the results of differential expression analysis onto the inferred cell-cell communications to easily manage/subset the ligand-receptor pairs of interest
net <- netMappingDEG(cellchat, features.name = features.name)
# extract the ligand-receptor pairs with upregulated ligands in WTg
net.up <- subsetCommunication(cellchat, net = net, datasets = "WTg",ligand.logFC = 0.2, receptor.logFC = 0.2)
# extract the ligand-receptor pairs with upregulated ligands and upregulated recetptors in EL, i.e.,downregulated in WTg
net.down <- subsetCommunication(cellchat, net = net, datasets = "EL",ligand.logFC = -0.2, receptor.logFC = -0.2)
gene.up <- extractGeneSubsetFromPair(net.up, cellchat)
gene.down <- extractGeneSubsetFromPair(net.down, cellchat)

# Fig.4D
pairLR.use.up = net.up[, "interaction_name", drop = F]
gg1 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.up, sources.use = 1, targets.use = c(1:5), comparison = c(1, 2),  angle.x = 90, remove.isolate = T,title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
#> Comparing communications on a merged object
pairLR.use.down = net.down[, "interaction_name", drop = F]
gg2 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.down, sources.use = 1, targets.use = c(1:5), comparison = c(1, 2),  angle.x = 90, remove.isolate = T,title.name = paste0("Up-regulated signaling in ", names(object.list)[1]))
#> Comparing communications on a merged object
gg1 + gg2

# Fig.S9D
pairLR.use.up = net.up[, "interaction_name", drop = F]
gg1 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.up, sources.use = c(1:5), targets.use = c(1:5), comparison = c(1, 2),  angle.x = 90, remove.isolate = T,title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
#> Comparing communications on a merged object
pairLR.use.down = net.down[, "interaction_name", drop = F]
gg2 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.down, sources.use = c(1:5), targets.use = c(1:5), comparison = c(1, 2),  angle.x = 90, remove.isolate = T,title.name = paste0("Up-regulated signaling in ", names(object.list)[1]))
#> Comparing communications on a merged object
gg1 + gg2

# Chord diagram
# Fig.4C
par(mfrow = c(1,2), xpd=TRUE)
netVisual_chord_gene(object.list[[2]], sources.use = 1, targets.use = c(1:5), slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = 1, targets.use = c(1:5), slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[1]))

netVisual_chord_gene(object.list[[2]], sources.use = 2, targets.use = c(1:5), slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = 2, targets.use = c(1:5), slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[1]))

netVisual_chord_gene(object.list[[2]], sources.use = 3, targets.use = c(1:5), slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = 3, targets.use = c(1:5), slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[1]))

netVisual_chord_gene(object.list[[2]], sources.use = 4, targets.use = c(1:5), slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = 4, targets.use = c(1:5), slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[1]))

netVisual_chord_gene(object.list[[2]], sources.use = 5, targets.use = c(1:5), slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = 5, targets.use = c(1:5), slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[1]))

# Fig.S9C
netVisual_chord_gene(object.list[[2]], sources.use = c(1:5), targets.use = 2, slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = c(1:5), targets.use = 2, slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[1]))

## Visually compare cell-cell communication using Hierarchy plot, Circle plot or Chord diagram
# circle
pathways.show <- c("PARs") 
weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) # control the edge weights across different datasets
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "circle", edge.weight.max = weight.max[1], edge.width.max = 10, signaling.name = paste(pathways.show, names(object.list)[i]))
}

# Fig.S10B
# heatmap
par(mfrow = c(1,2), xpd=TRUE)
ht <- list()
for (i in 1:length(object.list)) {
  ht[[i]] <- netVisual_heatmap(object.list[[i]], signaling = pathways.show, color.heatmap = "Reds",title.name = paste(pathways.show, "signaling ",names(object.list)[i]))
}
ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))

# Fig.S10A
# Chord diagram
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "chord", signaling.name = paste(pathways.show, names(object.list)[i]))
}

# Fig.S10C
## Compare the signaling gene expression distribution between different datasets
cellchat@meta$datasets = factor(cellchat@meta$datasets, levels = c("EL", "WTg")) # set factor level
plotGeneExpression(cellchat, signaling = "PARs", split.by = "datasets", colors.ggplot = T)


####################
##### NicheNetR ####
####################

library(nichenetr)
library(Seurat)
library(tidyverse)
library(circlize)
library(pheatmap)

seuratObj <- readRDS("seurat_ACC_TME.rds")
DefaultAssay(seuratObj) <- "SCT"
Idents(object = seuratObj) <- "cell.label"

ligand_target_matrix <- readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds"))
lr_network <- readRDS(url("https://zenodo.org/record/7074291/files/lr_network_mouse_21122021.rds"))
weighted_networks <- readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final_mouse.rds"))

sender_celltypes = c("ACC", "iCAF", "myCAF", "TAM", "TEC")

seuratObj <- UpdateSeuratObject(seuratObj)
seuratObj <- PrepSCTFindMarkers(seuratObj)

DefaultAssay(seuratObj) <- "integrated"

nichenet_output = nichenet_seuratobj_aggregate(
    seurat_obj = seuratObj, 
    receiver = "iCAF", 
    condition_colname = "sample",
    condition_oi = "WTg",
    condition_reference = "EL", 
    sender = sender_celltypes, 
    ligand_target_matrix = ligand_target_matrix,
    lr_network = lr_network,
    weighted_networks = weighted_networks
)

DefaultAssay(seuratObj) <- "SCT"

DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()

# Fig.S10D
# WTg vs. EL
DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), split.by = "sample", cols = "RdYlBu") + RotatedAxis()

# Inferred active ligand-target links
nichenet_output$ligand_target_heatmap

top_ligands <- 10
top_targets <- 50

reduced_heatmap <- nichenet_output$ligand_target_matrix[1:top_ligands, 1:top_targets]

# Fig.4E
pheatmap(reduced_heatmap, colorRampPalette(c("whitesmoke","royalblue"))(100),  
         treeheight_row = 0, 
         treeheight_col = 0)


#############
#7. Velocity analysis
#############
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
library(hdf5r)
library(sctransform)
library(viridis)
library(pheatmap)
library(RColorBrewer)

DefaultAssay(seurat_renameR) <- "SCT"

seurat_renameR <- AddMetaData(seurat_renameR, seurat_renameR@assays$SCT@data["Prss1",], col.name="Prss1")
seurat_renameR <- AddMetaData(seurat_renameR, seurat_renameR@assays$SCT@data["Spink1",], col.name="Spink1")
seurat_renameR <- AddMetaData(seurat_renameR, seurat_renameR$Prss1/seurat_renameR$Spink1, col.name="Prss1_Spink1_ratio")

# Fig.5A
FeaturePlot(seurat_renameR, features = c("Prss1", "Spink1", "Prss1_Spink1_ratio"), order = T, pt.size = 2) & scale_color_viridis(option = "H")

saveRDS(seurat_renameR, file = "seurat_renameR_Prss1_Spink1_ratio.rds")

seurat_ACC <- subset(seurat_renameR, idents = "ACC")
Idents(object = seurat_ACC) <- "integrated_snn_res.0.5"

# Fig.5B
VlnPlot(seurat_ACC, features = c("Prss1_Spink1_ratio"), flip = T) + 
     scale_fill_brewer(palette = "Spectral")

# Fig.5C
n_cells <- FetchData(seurat_ACC, 
                     vars = c("ident", "orig.ident")) %>%
        dplyr::count(ident, orig.ident) %>%
        tidyr::spread(ident, n)  %>%
        rename_with(~paste0("C", .), -orig.ident)

n_cells <- column_to_rownames(n_cells, var = "orig.ident")
n_cells_ratio <- data.frame(n_cells/rowSums(n_cells))

n_cells_ratio <- t(n_cells_ratio)
pheatmap(n_cells_ratio, show_rownames = T,color = colorRampPalette(c("blue", "gray", "orange"))(100))


DefaultAssay(seurat_ACC) <- "RNA"
markers <- FindAllMarkers(object = seurat_ACC, 
                          only.pos = TRUE,
                          logfc.threshold = 0.25) 

write.table(markers, "FindAllMarkers_ACC.txt", sep = "\t", quote=F, col.names=T, append=T, row.names=F)