find_new_clusters

1 Parameters

suffix = ""
data_to_read = "./Data/acc_tpm_nCount_mito_no146_15k_cancercells.rds"

2 functions

library(stringi)
source_from_github(repositoy = "DEG_functions",version = "0.2.24")

3 Data

acc = readRDS(file = data_to_read)
acc_primary = readRDS(file = "./Data/acc_cancer_no146_primaryonly15k_cancercells.rds")

(message("reading '" %>% paste0(data_to_read %>% basename()) %>% paste0("'")))

reading 'acc_tpm_nCount_mito_no146_15k_cancercells.rds'

NULL

pathways_scores = fread(file = "./Data/ACC_Canonical_Pathway_Scores.txt",sep = ",") %>% as.matrix(rownames=1) %>% t() %>%  as.data.frame()
hallmark_scores = fread(file = "./Data/ACC_Hallmark_Pathway_Scores.txt",sep = ",") %>% as.matrix(rownames=1) %>% t() %>%  as.data.frame()

ln_list = c("ACC22.LN.P11", "ACC22.P12.LN","ACC7.P13")
ln_plates = FetchData(object = acc,vars = "orig.ident") %>% mutate(
  tumor_type = if_else(condition = orig.ident %in% ln_list
                       ,true = "LN"
                       ,false = "primary"))

ln_plates["orig.ident"] <-NULL
acc= AddMetaData(object = acc,metadata = ln_plates)

pathways_scores = cbind(pathways_scores,hallmark_scores)
pathways_scores = pathways_scores[ , colSums(is.na(pathways_scores))==0] #remove cols with NA
pathways_scores = pathways_scores [rownames(pathways_scores) %in% colnames(acc),] #remove cells not in dataset
pathways_scores =  pathways_scores[order(row.names(pathways_scores)),] #order cells like dataset

4 Dim reduction

# run-dim-reduction on genes:
acc <- FindVariableFeatures(acc, selection.method = "vst", nfeatures = 2000)

acc <- ScaleData(acc)

acc <- RunPCA(acc,verbose = F)
ElbowPlot(acc)

pathway_scores_assay <- CreateAssayObject(counts = pathways_scores %>% t()) #create an assay

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

acc[["pathway_scores"]] = pathway_scores_assay

Warning: Keys should be one or more alphanumeric characters followed by an underscore, setting key from pathway_scores_ to pathwayscores_

# run-dim-reduction:
acc <- FindVariableFeatures(acc, selection.method = "vst", nfeatures = 2000,assay = "pathway_scores")

acc <- ScaleData(acc,assay = "pathway_scores",features = rownames(acc[["pathway_scores"]]))

acc <- RunPCA(acc, features = rownames(acc[["pathway_scores"]]),assay = "pathway_scores",reduction.name = "PCA_pathway_scores",verbose = F)
ElbowPlot(acc,reduction =  "PCA_pathway_scores")

acc <- RunUMAP(acc, dims = 1:5,reduction ="PCA_pathway_scores",reduction.name = "pathway_scores_umap",verbose = F)

Warning: Cannot add objects with duplicate keys (offending key: UMAP_), setting key to 'pathway_scores_umap_'

5 acc umaps

print_tab(plt = DimPlot(acc,group.by = "patient.ident"),title = "gene expression")

gene expression

print_tab(plt = DimPlot(acc,reduction = "pathway_scores_umap",group.by = "patient.ident"),title = "pathways scores")

pathways scores

NA

gs=acc@assays$RNA@var.features

myoscore=apply(acc@assays$RNA@scale.data[intersect(c("TP63","TP73","CAV1","CDH3","KRT5","KRT14","ACTA2","TAGLN","MYLK","DKK3"),gs),],2,mean)

lescore=apply(acc@assays$RNA@scale.data[intersect(c("KIT","EHF","ELF5","KRT7","CLDN3","CLDN4","CD24","LGALS3","LCN2","SLPI"),gs),],2,mean)
acc=AddMetaData(acc,lescore-myoscore,"luminal_over_myo")

#set lum_or_myo metadata
luminal_over_myo = FetchData(object = acc,vars = "luminal_over_myo")
luminal_over_myo$lum_or_myo  = case_when(luminal_over_myo$luminal_over_myo >1~"lum",luminal_over_myo$luminal_over_myo <(-1)~"myo",TRUE~"NA")
luminal_over_myo$luminal_over_myo <-NULL
acc=AddMetaData(object = acc,metadata = luminal_over_myo,col.name = "lum_or_myo")

print_tab(plt = FeaturePlot(object = acc,features = "luminal_over_myo",reduction = "pathway_scores_umap"),title = "luminal_over_myo")

luminal_over_myo

print_tab(plt = DimPlot(acc,group.by = "lum_or_myo",cols = c("red","green","grey"),reduction = "pathway_scores_umap"),title = "cell type")

cell type

print_tab(plt = FeaturePlot(object = acc,features = "luminal_over_myo"),title = "luminal_over_myo")

luminal_over_myo

print_tab(plt = DimPlot(acc,group.by = "lum_or_myo",cols = c("red","green","grey")),title = "cell type")

cell type

NA

6 Genes

print_tab(plt = 
            FeaturePlot(object = acc,features = c("FGF1","FGF2","FGF11","FGF12"),reduction = "pathway_scores_umap")
          ,title = "FGF")

FGF

print_tab(plt = 
FeaturePlot(object = acc,features = c("FGF18","FGF20","FGF15","FGF23"),reduction = "pathway_scores_umap")
          ,title = "FGF")

FGF

Warning in FetchData.Seurat(object = object, vars = c(dims, “ident”, features), : The following requested variables were not found: FGF15

print_tab(plt = 
FeaturePlot(object = acc,features = c("EGF"),reduction = "pathway_scores_umap")
          ,title = "EGF")

EGF

print_tab(plt = 
FeaturePlot(object = acc,features = c("NOTCH1","NOTCH2","NOTCH3","NOTCH4"),reduction = "pathway_scores_umap")
          ,title = "NOTCH genes")

NOTCH genes

print_tab(plt = 
FeaturePlot(object = acc,features = c("HES4","HEY1","HEY2","NRARP"),reduction = "pathway_scores_umap")
          ,title = "NOTCH targets")

NOTCH targets

NA

7 All PC’s

for (i in 1:8) {
print_tab(plt = VizDimLoadings(acc, dims = i, reduction = "PCA_pathway_scores"),title = paste("PC", i))
}

PC 1

PC 2

PC 3

PC 4

PC 5

PC 6

PC 7

PC 8

NA

8 cycling cells clustring

hallmark_name = "HALLMARK_G2M_CHECKPOINT"
genesets  =GSEABase::getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc@assays$RNA@data[intersect(geneIds,var_features),],2,mean)
acc=AddMetaData(acc,score,hallmark_name)

print_tab(plt = FeaturePlot(acc, reduction = "umap",features = "HALLMARK_G2M_CHECKPOINT"),title = "by genes")

by genes

print_tab(plt = FeaturePlot(acc, reduction = "pathway_scores_umap",features = "HALLMARK_G2M_CHECKPOINT"),title = "by genes")

by genes

NA

9 UMAP clusters

acc <- FindNeighbors(acc, dims = 1:10,reduction = "PCA_pathway_scores")

Computing nearest neighbor graph
Computing SNN

acc <- FindClusters(acc, resolution = 0.15,graph.name = "pathway_scores_snn")

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

Number of nodes: 951
Number of edges: 33164

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.8799
Number of communities: 2
Elapsed time: 0 seconds

DimPlot(acc,reduction = "pathway_scores_umap")

FeaturePlot(object = acc,features = "nFeature_RNA",reduction = "pathway_scores_umap")

VlnPlot(acc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

10 pathways DEG

print_tab(plt = datatable(pathway_markers_deg),title = "all deg")
print_tab(plt = datatable(pathway_markers_deg %>% dplyr::filter(abs(avg_log2FC)>1)),title = "all deg FC>1")

11 genes DEG

rho_genes = c("CDC42",
"RHOQ",
"RHOJ",
"RHOUV",
'RHOU',
"RHOV",
"RAC1",
"RAC2",
"RAC3",
"RHOG",
"RHOBTB1",
'RHOBTB2',
'RHOBTB3',
'RHOH',
'RHOA',
'RHOB',
'RHOC',
'RND1',
'RND2',
'RND3',
'RHOF',
'RHOD',
'RHOF')

notch_genes = c("JAG1","JAG2","NOTCH3","NOTCH2","NOTCH1","DLL1","MYB","HES4","HEY1","HEY2","NRARP")

print_tab(plt = acc_deg %>% datatable(class = 'cell-border stripe') ,title = "all DEG")
print_tab(plt = acc_deg[rownames(acc_deg) %in%  rho_genes,]%>% datatable(class = 'cell-border stripe')   ,title = "Rho genes in DEG")
print_tab(plt = acc_deg[rownames(acc_deg) %in%  notch_genes,] %>% datatable(class = 'cell-border stripe') ,title = "NOTCH genes in DEG")

DefaultAssay(object = acc)<- "pathway_scores"
genes_expression = FetchData(object = acc,vars = c("REACTOME-RHO-GTPASE-CYCLE","WP-HEAD-AND-NECK-SQUAMOUS-CELL-CARCINOMA"))
cor(genes_expression)

top_correlated <- function(dataset, genes, threshold,anti_cor = F) {
  markers_expression = FetchData(object = dataset,vars = genes,slot = "data") #get genes expression
  markers_average = rowMeans(markers_expression) %>% as.data.frame() %>% rename("average" = 1) #average them
  cor_mat = cor(expression %>% t(), markers_average)%>% as.data.frame() #cor with all genes
  cor_mat = cor_mat[complete.cases(cor_mat),,drop=F]  %>% as.data.frame %>%  rename("corr" = 1) #remove rows with NA in at least one col
  if (threshold<1){ #if threshold is based on pearson correlation 
      if(anti_cor == T){top_genes =   cor_mat %>% as.data.frame %>% select(1) %>% dplyr::filter(.< threshold) %>% rownames()}else{
          top_genes =   cor_mat %>% as.data.frame %>% select(1) %>% dplyr::filter(.> threshold) %>% rownames()
      }
  }else{ #if threshold is based on top correlated genes 
      if(anti_cor == T){threshold  = threshold*(-1)}
      top_genes =   cor_mat %>%  top_n(threshold,corr) %>% rownames()
      }
  return(top_genes)
}
expression = GetAssayData(object = acc,assay = "RNA",slot = "data") %>% as.data.frame()

top_correlated(dataset = acc,genes = "RND3",threshold = 20)

#UMAPS

print_tab(plt = 
            FeaturePlot(object = acc,features = "NOTCH2",reduction = "pathway_scores_umap")
 ,title = "NOTCH2 UMAP")

print_tab(plt = 
            FeaturePlot(object = acc,features = "RND3",reduction = "pathway_scores_umap")
 ,title = "RND3 UMAP")

print_tab(plt = 
            FeaturePlot(object = acc,features = "JAG1",reduction = "pathway_scores_umap")
 ,title = "JAG1 UMAP")

print_tab(plt = 
            FeaturePlot(object = acc,features = "JAG2",reduction = "pathway_scores_umap")
 ,title = "JAG2 UMAP")

print_tab(plt = 
            FeaturePlot(object = acc,features = "DLL1",reduction = "pathway_scores_umap")
 ,title = "DLL1 UMAP")

print_tab(plt = 
            FeaturePlot(object = acc,features = "HES4",reduction = "pathway_scores_umap")
 ,title = "HES4 UMAP")

12 HEAD-AND-NECK-SQUAMOUS

FeaturePlot(object = acc,features = "WP-HEAD-AND-NECK-SQUAMOUS-CELL-CARCINOMA",reduction = "pathway_scores_umap")

13 ACC1/2

ACC1_genes = c("MYC", "SOX6", "SOX8", "CTNND2", "NOTCH3","BCL2")
ACC2_genes = c("TP63","COL17A1","PDGFA", "DKK3","EGFR", "AXL","PDGFRA")

gs=acc@assays$RNA@var.features

acc1_score=apply(acc@assays$RNA@data[ACC1_genes,],2,mean)

acc2_score=apply(acc@assays$RNA@data[ACC2_genes,],2,mean)
acc=AddMetaData(acc,acc1_score-acc2_score,"acc1_over_acc2")

FeaturePlot(object = acc,features = "acc1_over_acc2",reduction = "pathway_scores_umap")

14 More clusters

acc <- FindNeighbors(acc, dims = 1:10,reduction = "PCA_pathway_scores")
acc <- FindClusters(acc, resolution = 0.5,graph.name = "pathway_scores_snn")

print_tab(plt = DimPlot(acc,reduction = "pathway_scores_umap"),title = "UMAP")
DimPlot(acc,group.by = "lum_or_myo",cols = c("red","green","grey"),reduction = "pathway_scores_umap")
FeaturePlot(object = acc,features = "nFeature_RNA",reduction = "pathway_scores_umap")

acc <- FindNeighbors(acc, dims = 1:10)
acc <- FindClusters(acc, resolution = 0.5)

print_tab(plt = DimPlot(acc),title = "UMAP")
DimPlot(acc,group.by = "lum_or_myo",cols = c("red","green","grey"))

for (i in 0:4) {
print_tab(plt = deg_markers %>% dplyr::filter(cluster == i) %>% datatable(),title = paste("cluster",i))
  }

14.1 Genes in DEG

print_tab(plt = acc_deg %>% dplyr::filter(gene %in% rho_genes),title = "Rho genes",subtitle_num = 3)
print_tab(plt = acc_deg %>% dplyr::filter(gene %in% notch_genes),title = "Notch genes",subtitle_num = 3)

15 MET pathway score

FeaturePlot(object = acc,features = "BIOCARTA-MET-PATHWAY",reduction = "pathway_scores_umap",max.cutoff = 0.3)

Warning: Could not find BIOCARTA-MET-PATHWAY in the default search locations, found in pathway_scores assay instead

16 MET pathway score based on gene expression

met_genes = msigdbr(species = "Homo sapiens",category = "C2") %>% dplyr::filter(gs_subcat != "CGP") %>%dplyr::filter(gs_name == "BIOCARTA_MET_PATHWAY") %>%  dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame() %>% pull(gene_symbol)

geneIds= met_genes
score <- apply(acc@assays$RNA@data[geneIds,],2,mean)
acc=AddMetaData(acc,score,"BIOCARTA_MET_PATHWAY")
FeaturePlot(object = acc,features = "BIOCARTA_MET_PATHWAY" ,reduction =  "pathway_scores_umap")

pathways_scores[,"BIOCARTA_MET_PATHWAY"]
pathways_scores[,"BIOCARTA_MET_PATHWAY"] %>% sort()

17 nFeatures of met genes

acc <- subset(acc, subset = nFeature_RNA > 2000 & percent.mt < 20)