Harmony_cancer

1 Functions

2 Data

library(GSEABase)
GBMSeurat_cancer = readRDS("./Data/GSM3828673_10X_GBM_seurat_cancer.RDS")
metaModuleGenes <- read.table(file = "./Data/MetaModuleGenes_CellPaper.tsv", header = TRUE)
# Creating genesets for the meta module markers
MesLike2GeneSet <- GeneSet(metaModuleGenes[,1])
MesLike1GeneSet <- GeneSet(metaModuleGenes[,2])
ACLikeGeneSet <- GeneSet(metaModuleGenes[!is.na(metaModuleGenes[,3]), 3])
OPCLikeGeneSet <- GeneSet(metaModuleGenes[,4])
NPCLike1GeneSet <- GeneSet(metaModuleGenes[,5])
NPCLike2GeneSet <- GeneSet(metaModuleGenes[,6])

MesLike1GeneSet@setName <- "MesLike1"
MesLike2GeneSet@setName <- "MesLike2"
ACLikeGeneSet@setName <- "ACLike"
OPCLikeGeneSet@setName <- "OPCLike"
NPCLike1GeneSet@setName <- "NPCLike1"
NPCLike2GeneSet@setName <- "NPCLike2"

library(harmony)
GBMSeurat_cancer <- GBMSeurat_cancer %>% 
    RunHarmony("orig.ident", plot_convergence = TRUE)

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 7 iterations
Warning: Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.RNA.harmony; see ?make.names for more details on syntax validity

GBMSeurat_cancer_harmony <- GBMSeurat_cancer %>% 
    RunUMAP(reduction = "harmony", dims = 1:15) %>% 
    FindNeighbors(reduction = "harmony", dims = 1:15) %>% 
    FindClusters(resolution = 0.5) %>% 
    identity()

17:59:38 UMAP embedding parameters a = 0.9922 b = 1.112
17:59:38 Read 10000 rows and found 15 numeric columns
17:59:38 Using Annoy for neighbor search, n_neighbors = 30
17:59:38 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:59:39 Writing NN index file to temp file /tmp/RtmpbskuF3/file186136a8df094
17:59:39 Searching Annoy index using 1 thread, search_k = 3000
17:59:42 Annoy recall = 100%
17:59:43 Commencing smooth kNN distance calibration using 1 thread
17:59:45 Initializing from normalized Laplacian + noise
17:59:45 Commencing optimization for 500 epochs, with 429308 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
17:59:56 Optimization finished
Computing nearest neighbor graph
Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 10000
Number of edges: 375295

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9160
Number of communities: 14
Elapsed time: 1 seconds

3 UMAP after integration

 DimPlot(object = GBMSeurat_cancer_harmony, reduction = "umap", pt.size = .1, group.by = "orig.ident")

 DimPlot(object = GBMSeurat_cancer_harmony, reduction = "umap", pt.size = .1, group.by = "seurat_clusters",label = T)

cancer_markers = list(
  MesLike1 = MesLike1GeneSet@geneIds,
  MesLike2 = MesLike2GeneSet@geneIds,
  ACLike = ACLikeGeneSet@geneIds,
  OPCLike = OPCLikeGeneSet@geneIds,
  NPCLike1 = NPCLike1GeneSet@geneIds,
  NPCLike2 = NPCLike2GeneSet@geneIds
)

GBMSeurat_cancer <- AddModuleScore(
  object = GBMSeurat_cancer,
  features = cancer_markers,ctrl = 50,name = "cancer_markers")

Warning: The following features are not present in the object: ERO1L, not searching for symbol synonyms
Warning: The following features are not present in the object: PPAP2B, not searching for symbol synonyms
Warning: The following features are not present in the object: SOX2-OT, LPPR1, not searching for symbol synonyms
Warning: The following features are not present in the object: CD24, GPR56, not searching for symbol synonyms
Warning: The following features are not present in the object: CD24, HMP19, LOC150568, not searching for symbol synonyms

GBMSeurat_cancer_harmony[["MesLike1"]] = GBMSeurat_cancer$cancer_markers1
GBMSeurat_cancer_harmony[["MesLike2"]] = GBMSeurat_cancer$cancer_markers2
GBMSeurat_cancer_harmony[["ACLike"]] = GBMSeurat_cancer$cancer_markers3
GBMSeurat_cancer_harmony[["OPCLike"]] = GBMSeurat_cancer$cancer_markers4
GBMSeurat_cancer_harmony[["NPCLike1"]] = GBMSeurat_cancer$cancer_markers5
GBMSeurat_cancer_harmony[["NPCLike2"]] = GBMSeurat_cancer$cancer_markers6


cancer_scores  = FetchData(object = GBMSeurat_cancer_harmony,vars = c("MesLike1","MesLike2","ACLike","OPCLike","NPCLike1","NPCLike2"))
cancer_scores$assignment<- colnames(cancer_scores)[max.col(cancer_scores, ties.method = "first")]

GBMSeurat_cancer_harmony %<>% AddMetaData(metadata = cancer_scores$assignment,col.name = "cancer_type")

DimPlot(GBMSeurat_cancer_harmony,group.by = "cancer_type",label = T)

4 Feature plots

# Testing for clusters enriched in markers of the different malignant meta-modules
print_tab(FeaturePlot(GBMSeurat_cancer_harmony, features = "MesLike1",reduction = "umap"),title = "MesLike1")

MesLike1

print_tab(FeaturePlot(GBMSeurat_cancer_harmony, features = "MesLike2",reduction = "umap"),title = "MesLike2")

MesLike2

print_tab(FeaturePlot(GBMSeurat_cancer_harmony, features = "ACLike",reduction = "umap"),title = "ACLike")

ACLike

print_tab(FeaturePlot(GBMSeurat_cancer_harmony, features = "OPCLike",reduction = "umap"),title = "OPCLike")

OPCLike

print_tab(FeaturePlot(GBMSeurat_cancer_harmony, features = "NPCLike1",reduction = "umap"),title = "NPCLike1")

NPCLike1

print_tab(FeaturePlot(GBMSeurat_cancer_harmony, features = "NPCLike2",reduction = "umap"),title = "NPCLike2")

NPCLike2

NA

5 Analysis

clusters_and_scores = FetchData(object = GBMSeurat_cancer_harmony,vars= c("cancer_type","seurat_clusters")) %>%  group_by(seurat_clusters,cancer_type) %>%  summarise(n_cells = n())%>% mutate(per =  100 *n_cells/sum(n_cells))

`summarise()` has grouped output by 'seurat_clusters'. You can override using the `.groups` argument.

integration_score = clusters_and_scores %>%  group_by(seurat_clusters) %>% filter(n_cells == max(n_cells)) %>% pull(per) %>% mean() %>% round(digits = 2)

v_factor_levels <-c( "MesLike1", "MesLike2", "NPCLike1", "NPCLike2", "OPCLike","ACLike")
colors = RColorBrewer::brewer.pal(6, "Paired"); colors[5] = "orange"
p2 = ggplot(data=clusters_and_scores, aes(x=seurat_clusters, y=per, fill=factor(cancer_type, levels = v_factor_levels))) +
  geom_bar(stat="identity")+theme_minimal() + scale_fill_manual(values = colors,name  = "Cancer type")+ labs(title = "Seurat",subtitle = "integration score=" %s+% integration_score %s+% "%")+ylab("% from cluster")
p2

6 silhouette score

library(cluster, quietly = TRUE)
dist.matrix <- dist(x = Embeddings(object = GBMSeurat_cancer_harmony[["umap"]]))
clusters <- GBMSeurat_cancer_harmony$cancer_type
sil <- silhouette(x = as.numeric(x = as.factor(x = clusters)), dist = dist.matrix)
summary(sil)

Silhouette of 10000 units in 6 clusters from silhouette.default(x = as.numeric(x = as.factor(x = clusters)),  from     dist = dist.matrix) :
 Cluster sizes and average silhouette widths:
       2842        2968        1479        1984         537         190 
 0.03627807  0.13422832 -0.16423981 -0.21607186  0.42678272 -0.12909305 
Individual silhouette widths:
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.881398 -0.315680  0.040535  0.003455  0.354209  0.573049 

GBMSeurat_cancer_harmony$sil = sil[,3]
VlnPlot(object = GBMSeurat_cancer_harmony,features = "sil",group.by = "cancer_type")

7 Combine cancer subtypes

GBMSeurat_cancer_harmony$cancer_type = GBMSeurat_cancer_harmony$cancer_type %>% gsub(pattern = "MesLike1|MesLike2",replacement = "MesLike")%>% gsub(pattern = "NPCLike1|NPCLike2",replacement = "NPCLike")

clusters_and_scores = FetchData(object = GBMSeurat_cancer_harmony,vars= c("cancer_type","seurat_clusters")) %>%  group_by(seurat_clusters,cancer_type) %>%  summarise(n_cells = n())%>% mutate(per =  100 *n_cells/sum(n_cells))

`summarise()` has grouped output by 'seurat_clusters'. You can override using the `.groups` argument.

integration_score = clusters_and_scores %>%  group_by(seurat_clusters) %>% filter(n_cells == max(n_cells)) %>% pull(per) %>% mean() %>% round(digits = 2)

v_factor_levels <-c( "MesLike", "NPCLike", "OPCLike","ACLike")
colors = RColorBrewer::brewer.pal(6, "Paired")[c(2,4,5,6)]; colors[3] = "orange"
p4 = ggplot(data=clusters_and_scores, aes(x=seurat_clusters, y=per, fill=factor(cancer_type, levels = v_factor_levels))) +
  geom_bar(stat="identity")+theme_minimal() + scale_fill_manual(values = colors,name  = "Cancer type")+ labs(title = "SiPSiC",subtitle = "integration score=" %s+% integration_score %s+% "%")+ylab("% from cluster")
p4

library(cluster, quietly = TRUE)
dist.matrix <- dist(x = Embeddings(object = GBMSeurat_cancer_harmony[["umap"]]))
clusters <- GBMSeurat_cancer_harmony$cancer_type
sil <- silhouette(x = as.numeric(x = as.factor(x = clusters)), dist = dist.matrix)
summary(sil)

Silhouette of 10000 units in 4 clusters from silhouette.default(x = as.numeric(x = as.factor(x = clusters)),  from     dist = dist.matrix) :
 Cluster sizes and average silhouette widths:
       2842        4447        2521         190 
 0.08042854  0.04773834  0.31742084 -0.09328014 
Individual silhouette widths:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.7576 -0.1140  0.2318  0.1223  0.3516  0.5490 

GBMSeurat_cancer_harmony$sil = sil[,3]
VlnPlot(object = GBMSeurat_cancer_harmony,features = "sil",group.by = "cancer_type")

pdf(file = "./Figures/Harmony/UMAP_orig.ident.pdf", onefile = FALSE)
DimPlot(GBMSeurat_cancer_harmony, reduction = "umap", group.by = "orig.ident")
dev.off()


pdf(file = "./Figures/Harmony/UMAP_clustering.pdf", onefile = FALSE)
DimPlot(GBMSeurat_cancer_harmony, reduction = "umap")
dev.off()


pdf(file = "./Figures/Harmony/macrophageScore.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "macrophageScore"))

dev.off()

pdf(file = "./Figures/Harmony/TCellScores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "TCellScores"))

dev.off()
pdf(file = "./Figures/Harmony/oligodendrocyteScores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "oligodendrocyteScores"))

dev.off()
pdf(file = "./Figures/Harmony/oligodendrocyteScores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "oligodendrocyteScores"))

dev.off()
pdf(file = "./Figures/Harmony/MESLike1Scores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "MESLike1Scores"))

dev.off()
pdf(file = "./Figures/Harmony/MESLike2Scores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "MESLike2Scores"))

dev.off()
pdf(file = "./Figures/Harmony/ACLikeScores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "ACLikeScores"))

dev.off()
pdf(file = "./Figures/Harmony/OPCLikeScores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "OPCLikeScores"))

dev.off()


pdf(file = "./Figures/Harmony/NPCLike1Scores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "NPCLike1Scores"))

dev.off()

pdf(file = "./Figures/Harmony/NPCLike2Scores.pdf", onefile = FALSE)
print(FeaturePlot(GBMSeurat_cancer_harmony,reduction = "umap", features = "NPCLike2Scores"))

dev.off()

# Measuring the percentage of cells of each malignant meta module that are in a single cluster
clusterAssignments <- Idents(GBMSeurat)
ScoresAndClusters <- cbind(clusterAssignments, MacrophageScores = macrophageScores$pathwayScores, OligodecdrocytesScores = oligodendrocyteScores$pathwayScores,
                           MESLike1Scores = MesLike1Scores$pathwayScores, MESLike2Scores = MesLike2Scores$pathwayScores, 
                           ACLikeScores = ACLikeScores$pathwayScores, OPCLikeScores = OPCLikeScores$pathwayScores, 
                           NPCLike1Scores = NPCLike1Scores$pathwayScores, NPCLike2Scores = NPCLike2Scores$pathwayScores, 
                           TCellScores = TCellScores$pathwayScores, as.data.frame(GBMSeurat@meta.data$orig.ident))
colnames(ScoresAndClusters)[colnames(ScoresAndClusters) == "GBMSeurat@meta.data$orig.ident"] <- "Patient"

# Checking the proportions of the different patients in each cluster
clusterByPatientIDs <- ScoresAndClusters[,c("clusterAssignments", "Patient")]
clusterByPatientIDs <- clusterByPatientIDs[order(clusterByPatientIDs$clusterAssignments),]
for(currCluster in 0:max(as.numeric(as.character(clusterByPatientIDs[,"clusterAssignments"]))))
{
  currClusterCells <- clusterByPatientIDs[clusterByPatientIDs[,"clusterAssignments"] == currCluster,"Patient"]
  cat("\n", "Cluster number: ", currCluster, " Total Cells: ", sum(clusterByPatientIDs[,"clusterAssignments"] == currCluster), "\n")
  print(table(currClusterCells))
}


NPCThreshold <- ACLikeThreshold <- OPCLikeThreshold <- MESLikesThreshold <- 0.13
OligodendrocyteThreshold <- 1.5
MacrophageThreshold <- TCellThreshold <- 0.05

macrophageOverThresh <- ScoresAndClusters[ScoresAndClusters[,"MacrophageScores"] > MacrophageThreshold,c("clusterAssignments", "MacrophageScores")]
print("Macrophage cells: Total over threshold = ")
print(nrow(macrophageOverThresh))
print(table(macrophageOverThresh[,"clusterAssignments"]))
cat("\n\n")

TCellsOverThresh <- ScoresAndClusters[ScoresAndClusters[,"TCellScores"] > TCellThreshold,c("clusterAssignments", "TCellScores")]
print("T cells: Total over threshold = ")
print(nrow(TCellsOverThresh))
print(table(TCellsOverThresh[,"clusterAssignments"]))
cat("\n\n")

OligodendrocytesOverThresh <- ScoresAndClusters[ScoresAndClusters[,"OligodecdrocytesScores"] > OligodendrocyteThreshold,c("clusterAssignments", "OligodecdrocytesScores")]
print("Oligodendrocyte cells: Total over threshold = ")
print(nrow(OligodendrocytesOverThresh))
print(table(OligodendrocytesOverThresh[,"clusterAssignments"]))
cat("\n\n")

NPCs1OverThresh <- ScoresAndClusters[ScoresAndClusters[,"NPCLike1Scores"] > NPCThreshold,c("clusterAssignments", "NPCLike1Scores")]
print("NPC1 cells: Total over threshold = ")
print(nrow(NPCs1OverThresh))
print(table(NPCs1OverThresh[,"clusterAssignments"]))
cat("\n\n")

NPCs2OverThresh <- ScoresAndClusters[ScoresAndClusters[,"NPCLike2Scores"] > NPCThreshold,c("clusterAssignments", "NPCLike2Scores")]
print("NPC2 cells: Total over threshold = ")
print(nrow(NPCs2OverThresh))
print(table(NPCs2OverThresh[,"clusterAssignments"]))
cat("\n\n")

ACLikesOverThresh <- ScoresAndClusters[ScoresAndClusters[,"ACLikeScores"] > ACLikeThreshold,c("clusterAssignments", "ACLikeScores")]
print("ACLike cells: Total over threshold = ")
print(nrow(ACLikesOverThresh))
print(table(ACLikesOverThresh[,"clusterAssignments"]))
cat("\n\n")

OPCLikesOverThresh <- ScoresAndClusters[ScoresAndClusters[,"OPCLikeScores"] > OPCLikeThreshold,c("clusterAssignments", "OPCLikeScores")]
print("OPCLike cells: Total over threshold = ")
print(nrow(OPCLikesOverThresh))
print(table(OPCLikesOverThresh[,"clusterAssignments"]))
cat("\n\n")

MES1OverThresh <- ScoresAndClusters[ScoresAndClusters[,"MESLike1Scores"] > MESLikesThreshold,c("clusterAssignments", "MESLike1Scores")]
print("MES1 cells: Total over threshold = ")
print(nrow(MES1OverThresh))
print(table(MES1OverThresh[,"clusterAssignments"]))
cat("\n\n")

MES2OverThresh <- ScoresAndClusters[ScoresAndClusters[,"MESLike2Scores"] > MESLikesThreshold,c("clusterAssignments", "MESLike2Scores")]
print("MES2 cells: Total over threshold = ")
print(nrow(MES2OverThresh))
print(table(MES2OverThresh[,"clusterAssignments"]))
cat("\n\n")

# Checking the distribution of each patient's malignant cells across clusters

isCellMalignant <- ScoresAndClusters[,"NPCLike1Scores"] > NPCThreshold |
                   ScoresAndClusters[,"NPCLike2Scores"] > NPCThreshold |
                   ScoresAndClusters[,"ACLikeScores"] > ACLikeThreshold |
                   ScoresAndClusters[,"OPCLikeScores"] > OPCLikeThreshold |
                   ScoresAndClusters[,"MESLike1Scores"] > MESLikesThreshold |
                   ScoresAndClusters[,"MESLike2Scores"] > MESLikesThreshold
malignantCellsOnly <- ScoresAndClusters[isCellMalignant,]

for (currPatient in levels(malignantCellsOnly$Patient))
{
  currPatientMaligCells <- malignantCellsOnly[malignantCellsOnly[,"Patient"] == currPatient,]
  cat("\n", "Patient number: ", currPatient, "\nTotal number of cells: ", nrow(currPatientMaligCells), "\n")
  print(table(currPatientMaligCells[,"clusterAssignments"]))
  
  # Finding the NPC-Like cells of this patient only
  currPatientNPCLike1 <- currPatientMaligCells[currPatientMaligCells[,"NPCLike1Scores"] > NPCThreshold,]
  currPatientNPCLike2 <- currPatientMaligCells[currPatientMaligCells[,"NPCLike2Scores"] > NPCThreshold,]
  
  cat("\n", "NPC-Like1 cells in this cluster:\n")
  cat("Total number of NPC-Like1: ", nrow(currPatientNPCLike1))
  print(table(currPatientNPCLike1[, "clusterAssignments"]))
  
  cat("\n\n", "NPC-Like2 cells in this cluster:\n")
  cat("Total number of NPC-Like2: ", nrow(currPatientNPCLike2))
  print(table(currPatientNPCLike2[, "clusterAssignments"]))
  
  cat("\n\n")
}

clusters_assignment_results = data.frame(cluster = 0:15)
clusters_assignment_results[,"Macrophages"] = table(macrophageOverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)
clusters_assignment_results[,"TCellS"] = table(TCellsOverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)
clusters_assignment_results[,"Oligodecdrocytes"] = table(OligodendrocytesOverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)
clusters_assignment_results[,"NPC1"] = table(NPCs1OverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)

clusters_assignment_results[,"NPC2"] = table(NPCs2OverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)
clusters_assignment_results[,"ACLike"] = table(ACLikesOverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)
clusters_assignment_results[,"OPCLike"] = table(OPCLikesOverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)

clusters_assignment_results[,"MES1"] = table(MES1OverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)
clusters_assignment_results[,"MES2"] = table(MES2OverThresh[,"clusterAssignments"]) %>% as.data.frame() %>% pull(Freq)

8 Session

session_info()