SiPSiC_CP

1 Functions

2 Data

calculatePathwayScores <- function(pathwayToScore, countsMatrix)
{
  pathwayName <- pathwayToScore@setName
  print(pathwayName)
  pathwayGenes <- pathwayToScore@geneIds
  pathwayScoresObject <- getPathwayScores(countsMatrix, pathwayGenes) 
  suppressWarnings(
    if(is.na(pathwayScoresObject)){return (NA)}
  )
  return(pathwayScoresObject$pathwayScores)
}
GBMSeurat_cancer = readRDS("./Data/GSM3828673_10X_GBM_seurat_cancer.RDS")
genesets_cp <- getGmt("./Data/c2.cp.v2023.2.Hs.symbols.gmt")
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"

pathwayScoreLists2 <- pathwayScoreLists[!na.rm(pathwayScoreLists)]

Error in na.rm(pathwayScoreLists) : could not find function "na.rm"

## Adding the pathway name to each list of pathway scores, as the names will be used during the clustering pre-processing steps
for(currPathwayIndex in 1:length(pathwayScoreLists)){
  rownames(pathwayScoresMatrix)[currPathwayIndex] <- genesets_cp@.Data[[currPathwayIndex]]@setName
}

GBMSeurat_cancer_sipsic_cp <- CreateSeuratObject(counts = pathwayScoresMatrix)

Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')

saveRDS(object = GBMSeurat_cancer_sipsic_cp,file = "./Data/GBMSeurat_cancer_sipsic_cp.RDS")

GBMSeurat_cancer_sipsic_cp = readRDS(file = "./Data/GBMSeurat_cancer_sipsic_cp.RDS")

# calculate cancer_markers on original data
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_sipsic_cp[["MesLike1"]] = GBMSeurat_cancer$cancer_markers1
GBMSeurat_cancer_sipsic_cp[["MesLike2"]] = GBMSeurat_cancer$cancer_markers2
GBMSeurat_cancer_sipsic_cp[["ACLike"]] = GBMSeurat_cancer$cancer_markers3
GBMSeurat_cancer_sipsic_cp[["OPCLike"]] = GBMSeurat_cancer$cancer_markers4
GBMSeurat_cancer_sipsic_cp[["NPCLike1"]] = GBMSeurat_cancer$cancer_markers5
GBMSeurat_cancer_sipsic_cp[["NPCLike2"]] = GBMSeurat_cancer$cancer_markers6

# further pre-processing steps before moving on to clustering
# GBMSeurat_cancer_sipsic_cp <- FindVariableFeatures(GBMSeurat_cancer_sipsic_cp, selection.method = "vst", nfeatures = 100)
GBMSeurat_cancer_sipsic_cp <- ScaleData(GBMSeurat_cancer_sipsic_cp, rownames(GBMSeurat_cancer_sipsic_cp))
GBMSeurat_cancer_sipsic_cp <- RunPCA(GBMSeurat_cancer_sipsic_cp,features =  rownames(GBMSeurat_cancer_sipsic_cp))
# print(GBMSeurat_cancer_sipsic_cp[["pca"]], dims=1:5, nfeatures = 5)
print(ElbowPlot(GBMSeurat_cancer_sipsic_cp))

# Clustering based on the top principal components
GBMSeurat_cancer_sipsic_cp <- FindNeighbors(GBMSeurat_cancer_sipsic_cp, dims = 1:15)
GBMSeurat_cancer_sipsic_cp <- FindClusters(GBMSeurat_cancer_sipsic_cp, resolution = 0.65)
GBMSeurat_cancer_sipsic_cp <- RunUMAP(GBMSeurat_cancer_sipsic_cp, dims = 1:15)

3 Feature plots

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

MESLike1Scores

print_tab(FeaturePlot(GBMSeurat_cancer_sipsic_cp, features = "MesLike2"),title = "MESLike2Scores")

MESLike2Scores

print_tab(FeaturePlot(GBMSeurat_cancer_sipsic_cp, features = "ACLike"),title = "ACLikeScores")

ACLikeScores

print_tab(FeaturePlot(GBMSeurat_cancer_sipsic_cp, features = "OPCLike"),title = "OPCLikeScores")

OPCLikeScores

print_tab(FeaturePlot(GBMSeurat_cancer_sipsic_cp, features = "NPCLike1"),title = "NPCLike1Scores")

NPCLike1Scores

print_tab(FeaturePlot(GBMSeurat_cancer_sipsic_cp, features = "NPCLike2"),title = "NPCLike2Scores")

NPCLike2Scores

NA

4 UMAPS

print_tab(DimPlot(GBMSeurat_cancer_sipsic_cp, reduction = "umap"),title = "clusters")

clusters

print_tab(DimPlot(GBMSeurat_cancer_sipsic_cp, reduction = "umap",group.by = "orig.ident"),title = "by patient")

by patient

print_tab(DimPlot(GBMSeurat_cancer_sipsic_cp, reduction = "umap",group.by = "cancer_type"),title = "by cell types")

by cell types

NA

5 stacked batplot

clusters_and_scores = FetchData(object = GBMSeurat_cancer_sipsic_cp,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"
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 CP",subtitle = "integration score=" %s+% integration_score %s+% "%")+ylab("% from cluster")
p4

6 silhouette score

library(cluster, quietly = TRUE)
dist.matrix <- dist(x = Embeddings(object = GBMSeurat_cancer_sipsic_cp[["umap"]]))
clusters <- GBMSeurat_cancer_sipsic_cp$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.05610227  0.04080160 -0.15085952 -0.23096790  0.42317404 -0.16351472 
Individual silhouette widths:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.81349 -0.31760 -0.01108 -0.05235  0.18586  0.55330 

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

7 Combine cancer subtypes

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

clusters_and_scores = FetchData(object = GBMSeurat_cancer_sipsic_cp,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 = "original",subtitle = "integration score=" %s+% integration_score %s+% "%")+ylab("% from cluster")
p4

library(cluster, quietly = TRUE)
dist.matrix <- dist(x = Embeddings(object = GBMSeurat_cancer_sipsic_cp[["umap"]]))
clusters <- GBMSeurat_cancer_sipsic_cp$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.03131561 -0.03240623  0.31179352 -0.04210191 
Individual silhouette widths:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.68597 -0.15797  0.15271  0.07229  0.27370  0.50481 

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

print(DimPlot(GBMSeurat_cancer_sipsic_cp, reduction = "umap",group.by = "cancer_type"))