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)

acc <- RunUMAP(acc, dims = 1:5)

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
11:02:06 UMAP embedding parameters a = 0.9922 b = 1.112
11:02:06 Read 951 rows and found 5 numeric columns
11:02:06 Using Annoy for neighbor search, n_neighbors = 30
11:02:06 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:02:06 Writing NN index file to temp file /tmp/Rtmpp1qPPG/filea27f57b560b6
11:02:06 Searching Annoy index using 1 thread, search_k = 3000
11:02:06 Annoy recall = 100%
11:02:07 Commencing smooth kNN distance calibration using 1 thread
11:02:09 Initializing from normalized Laplacian + noise
11:02:09 Commencing optimization for 500 epochs, with 34634 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:02:11 Optimization finished

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

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.1,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.9052
Number of communities: 2
Elapsed time: 0 seconds

DimPlot(acc,reduction = "pathway_scores_umap")

10 pathways DEG

print_tab(plt = datatable(pathway_markers_deg),title = "all deg")

all deg

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

all deg FC>1

NA

11 genes DEG

acc = SetIdent(object = acc,value = "pathway_scores_snn_res.0.1")
genes_markers = FindMarkers(object = acc,ident.1 = "0",ident.2 = "1",assay = "RNA",min.cells.feature = 10,logfc.threshold = 0,densify = T)
acc_deg  = genes_markers %>% mutate(fdr = p.adjust(p_val,method = "fdr"))%>% #add fdr
    dplyr::filter((avg_log2FC>1 & fdr<0.05) | (avg_log2FC< (-1) & fdr<0.05))  #filter significant

# enrichment_analysis(Differential_expression_genes = genes_markers,background = rownames(acc),fdr_Cutoff = 0.05,ident.1 = "0",ident.2 = "1")


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


down =  acc_deg %>% 
    dplyr::filter((avg_log2FC>1.5 & fdr<0.05)) %>% 
                    rownames()

genes_vec_enrichment(genes = down,background = rownames(acc),homer = T,title = "test",custom_pathways = cp_genes)

up =  acc_deg %>% 
    dplyr::filter((avg_log2FC<1.5 & fdr<0.05)) %>% 
                    rownames()

genes_vec_enrichment(genes = up,background = rownames(acc),homer = T,title = "test",custom_pathways = cp_genes)

# genes = msigdbr(species = "Homo sapiens", category = "C2") %>% dplyr::filter(gs_subcat != "CGP" & grepl('REACTOME_RHO_GTPASE_CYCLE', gs_name)) %>% pull("gene_symbol") %>% unique()


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")

all DEG

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

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")

NOTCH genes in DEG

NA

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")

NOTCH2 UMAP

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

RND3 UMAP

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

JAG1 UMAP

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

JAG2 UMAP

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

DLL1 UMAP

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

HES4 UMAP

NA

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")

UMAP

NA

all_markers = FindAllMarkers(object = acc,logfc.threshold = 0,densify = T,assay = "pathway_scores")

deg_markers=all_markers %>% 
  mutate(fdr = p.adjust(p_val,method = "fdr"))%>%  #add fdr 
  dplyr::filter(fdr<0.05) %>%   dplyr::filter(abs(avg_log2FC)>1)

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

cluster 0

cluster 1

cluster 2

cluster 3

cluster 4

NA