SiPSiC

1 Functions

2 Data

calculatePathwayScores <- function(pathwayToScore, countsMatrix)
{
  pathwayName <- pathwayToScore@setName
  pathwayGenes <- pathwayToScore@geneIds
  pathwayScoresObject <- getPathwayScores(countsMatrix, pathwayGenes) 
  return(pathwayScoresObject$pathwayScores)
}
GBMSeurat_cancer = readRDS("./Data/GSM3828673_10X_GBM_seurat_cancer.RDS")
genesets <- getGmt("./Data/h.all.v7.0.symbols.pluscc.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"

TPMCounts = GBMSeurat_cancer@assays$RNA@counts %>% as.matrix() %>% Matrix(sparse = T)

pathwayScoreLists <- lapply(X = genesets@.Data, calculatePathwayScores, TPMCounts)
pathwayScoresMatrix <- as.data.frame(do.call("rbind", pathwayScoreLists))

## 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@.Data[[currPathwayIndex]]@setName
}

GBMSeurat_cancer_sipsic <- CreateSeuratObject(counts = pathwayScoresMatrix)

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

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[["MesLike1"]] = GBMSeurat_cancer$cancer_markers1
GBMSeurat_cancer_sipsic[["MesLike2"]] = GBMSeurat_cancer$cancer_markers2
GBMSeurat_cancer_sipsic[["ACLike"]] = GBMSeurat_cancer$cancer_markers3
GBMSeurat_cancer_sipsic[["OPCLike"]] = GBMSeurat_cancer$cancer_markers4
GBMSeurat_cancer_sipsic[["NPCLike1"]] = GBMSeurat_cancer$cancer_markers5
GBMSeurat_cancer_sipsic[["NPCLike2"]] = GBMSeurat_cancer$cancer_markers6

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

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

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

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

3 Feature plots

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

MESLike1Scores

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

MESLike2Scores

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

ACLikeScores

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

OPCLikeScores

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

NPCLike1Scores

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

NPCLike2Scores

NA

4 UMAPS

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

clusters

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

by patient

print_tab(DimPlot(GBMSeurat_cancer_sipsic, 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,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"
p1 = 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")
p1

6 silhouette score

library(cluster, quietly = TRUE)
dist.matrix <- dist(x = Embeddings(object = GBMSeurat_cancer_sipsic[["umap"]]))
clusters <- GBMSeurat_cancer_sipsic$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.05078936  0.03399211 -0.11382633 -0.17191359  0.21005554 -0.23310521 
Individual silhouette widths:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.76769 -0.24510 -0.01811 -0.04844  0.16534  0.33000 

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

7 Combine cancer subtypes

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

clusters_and_scores = FetchData(object = GBMSeurat_cancer_sipsic,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_sipsic[["umap"]]))
clusters <- GBMSeurat_cancer_sipsic$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.03706650 -0.01902861  0.33220367 -0.13660378 
Individual silhouette widths:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.68356 -0.10855  0.11759  0.08323  0.30599  0.51275 

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