cnmf_programs_in_ACC_V2

all_acc_cancer_cells = readRDS("./Data/acc_cancer_cells_V3.RDS")

  original_myo_genes = c("TP63", "TP73", "CAV1", "CDH3", "KRT5", "KRT14", "ACTA2", "TAGLN", "MYLK", "DKK3")
  original_lum_genes = c("KIT", "EHF", "ELF5", "KRT7", "CLDN3", "CLDN4", "CD24", "LGALS3", "LCN2", "SLPI")
  notch_targets_genes = c("NRARP", "NOTCH3", "HES4", "HEY1", "HEY2")
  notch_ligands = c("DLL1", "JAG1", "JAG2", "HEY1")

  orig_myoscore=apply(acc_cancerCells_noACC1@assays[["RNA"]][original_myo_genes,],2,mean)
  orig_lescore=apply(acc_cancerCells_noACC1@assays[["RNA"]][original_lum_genes,],2,mean)
  notch_score=apply(acc_cancerCells_noACC1@assays[["RNA"]][notch_targets_genes,],2,mean)
  notch_ligands_score=apply(acc_cancerCells_noACC1@assays[["RNA"]][notch_ligands,],2,mean)

  acc_cancerCells_noACC1 = AddMetaData(object = acc_cancerCells_noACC1,metadata = orig_lescore-orig_myoscore,col.name = "lum_over_myo")
    acc_cancerCells_noACC1 = AddMetaData(object = acc_cancerCells_noACC1,metadata = orig_myoscore,col.name = "myo_score")
    acc_cancerCells_noACC1 = AddMetaData(object = acc_cancerCells_noACC1,metadata = orig_lescore,col.name = "lum_score")
    acc_cancerCells_noACC1 = AddMetaData(object = acc_cancerCells_noACC1,metadata = notch_score,col.name = "notch_score")
    
  acc_cancerCells_noACC1 = AddMetaData(object = acc_cancerCells_noACC1,metadata = notch_ligands_score,col.name = "notch_ligands_score_score")

Load object

reticulate::repl_python()
from cnmf import cNMF
import pickle
nfeatures = "2K"
f = open('./Data/cNMF/HMSC_cNMF_' + nfeatures+ 'vargenes/cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

selected_k = 4
density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)

/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/lib/python3.7/site-packages/scanpy/preprocessing/_simple.py:843: UserWarning: Received a view of an AnnData. Making a copy.
  view_to_actual(adata)

usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

quit
gep_scores = py$gep_scores
gep_tpm = py$gep_tpm
all_metagenes= py$usage_norm

def get_norm_counts(counts, tpm,high_variance_genes_filter): #from cnmf.py
    import numpy as np
    import scipy.sparse as sp
    ## Subset out high-variance genes
    norm_counts = counts[:, high_variance_genes_filter].copy()
  
    ## Scale genes to unit variance
    if sp.issparse(tpm.X):
        sc.pp.scale(norm_counts, zero_center=False)
        if np.isnan(norm_counts.X.data).sum() > 0:
            print('Warning NaNs in normalized counts matrix')                       
    else:
        norm_counts.X /= norm_counts.X.std(axis=0, ddof=1)
        if np.isnan(norm_counts.X).sum().sum() > 0:
            print('Warning NaNs in normalized counts matrix')                    
    
  
    ## Check for any cells that have 0 counts of the overdispersed genes
    zerocells = norm_counts.X.sum(axis=1)==0
    if zerocells.sum()>0:
        examples = norm_counts.obs.index[zerocells]
        print('Warning: %d cells have zero counts of overdispersed genes. E.g. %s' % (zerocells.sum(), examples[0]))
        print('Consensus step may not run when this is the case')
    
    return(norm_counts)
  
def get_usage_from_score(counts,tpm, genes,cnmf_obj):
      import anndata as ad
      import scanpy as sc
      import numpy as np
      from sklearn.decomposition import non_negative_factorization
      import pandas as pd
      counts_adata = ad.AnnData(counts)
      tpm_adata = ad.AnnData(tpm)
      norm_counts = get_norm_counts(counts=counts_adata,tpm=tpm_adata,high_variance_genes_filter=np.array(genes)) #norm counts as cnmf
      spectra = cnmf_obj.get_median_spectra(k=4) #get score 
      spectra = spectra[spectra.columns.intersection(genes)] #remove genes not in @genes
      usage_by_calc,_,_ = non_negative_factorization(X=norm_counts.X, H = spectra.values, update_H=False,n_components = 4,max_iter=1000,init ="random")
      usage_by_calc = pd.DataFrame(usage_by_calc, index=counts.index, columns=spectra.index) #insert to df+add names
      usage_by_calc = usage_by_calc.div(usage_by_calc.sum(axis=1), axis=0) # sum rows to 1 
      return(usage_by_calc)

Get expression

nfeatures = 2000
acc1_cancer_cells = FindVariableFeatures(object = acc1_cancer_cells,nfeatures = nfeatures)
vargenes = VariableFeatures(object = acc1_cancer_cells)
hmsc_expression = t(as.matrix(GetAssayData(acc1_cancer_cells,slot='data')))
hmsc_expression = 2**hmsc_expression #convert from log2(tpm+1) to tpm
hmsc_expression = hmsc_expression-1
hmsc_expression = hmsc_expression[,!colSums(hmsc_expression==0, na.rm=TRUE)==nrow(hmsc_expression)] #delete rows that have all 0
hmsc_expression = hmsc_expression[,vargenes]
hmsc_expression = hmsc_expression  %>% as.data.frame()

Calculate usage

hmsc_expression = r.hmsc_expression
vargenes = r.vargenes
usage_by_calc = get_usage_from_score(counts=hmsc_expression,tpm=hmsc_expression,genes= vargenes,cnmf_obj=cnmf_obj)

all_metagenes= py$usage_by_calc
all_metagenes = all_metagenes[,c(1,4,3,2)]

Programs enrichments

#Hallmarks:
plt_list = list()
for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T)
   
  plt_list[[i]] = res$plt
}
p = ggarrange(plotlist  = plt_list)
print_tab(plt = p,title = "HALLMARK enrichment")

HALLMARK enrichment

#Canonical:
canonical_pathways = msigdbr(species = "Homo sapiens", category = "C2") %>% dplyr::filter(gs_subcat != "CGP") %>%  dplyr::distinct(gs_name, gene_symbol) 

plt_list = list()
for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T,custom_pathways = canonical_pathways)
   
  plt_list[[i]] = res$plt
}
p = ggarrange(plotlist  = plt_list)
print_tab(plt = p,title = "canonical pathway enrichment")

canonical pathway enrichment

NA

luminal and myo genes in score

# lum genes in metagenes
message("lum genes in metagenes:")

lum genes in metagenes:

for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  cat(paste0("metagene ",i,": "))
  print(original_lum_genes[original_lum_genes %in% top])

}

metagene 1: [1] "KRT7" "SLPI"
metagene 2: [1] "CLDN4"
metagene 3: character(0)
metagene 4: [1] "KRT7"   "LGALS3" "LCN2"   "SLPI"  

cat("\n")

# myo genes in metagenes
message("myo genes in metagenes:")

myo genes in metagenes:

for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  cat(paste0("metagene ",i,": "))
  print(original_myo_genes[original_myo_genes %in% top])

}

metagene 1: [1] "KRT14" "ACTA2" "DKK3" 
metagene 2: character(0)
metagene 3: character(0)
metagene 4: character(0)

cat("\n")

notch_genes = c("JAG1","JAG2","NOTCH3","NOTCH2","NOTCH1","DLL1","MYB","HES4","HEY1","HEY2","NRARP")
# notch genes in metagenes
message("notch genes in metagenes:")

notch genes in metagenes:

for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  cat(paste0("metagene ",i,": "))
  print(notch_genes[notch_genes %in% top])

}

metagene 1: character(0)
metagene 2: character(0)
metagene 3: character(0)
metagene 4: character(0)

cat("\n")

Metagenes on ACC

# Make metagene names
for (i in 1:ncol(all_metagenes)) {
  colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}


#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% select(i)
  acc_cancerCells_noACC1 = AddMetaData(object = acc_cancerCells_noACC1,metadata = metage_metadata)
}

print_tab(plt = 
            FeaturePlot(object = acc_cancerCells_noACC1,features = colnames(all_metagenes),max.cutoff = 100)
          ,title = "metagenes expression")

metagenes expression

all_metagenes= scale(all_metagenes) %>% as.data.frame()
for (i in 1:ncol(all_metagenes)) {
  colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}


#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% select(i)
  acc_cancerCells_noACC1 = AddMetaData(object = acc_cancerCells_noACC1,metadata = metage_metadata)
}

print_tab(plt = 
            FeaturePlot(object = acc_cancerCells_noACC1,features = colnames(all_metagenes),max.cutoff = 100)
          ,title = "metagenes z-score expression")

metagenes z-score expression

NA

Assign to programs by z-scored metagenes

acc_cancerCells_noACC1 = program_assignment(dataset = acc_cancerCells_noACC1,larger_by = 2,program_names = colnames(all_metagenes))

UMAPS

colors =  rainbow(acc_cancerCells_noACC1$program.assignment %>% unique() %>% length()-1)
colors = c(colors,"grey")
print_tab(plt = DimPlot(object = acc_cancerCells_noACC1,group.by = "program.assignment",cols =colors),title = "program.assignment")

program.assignment

print_tab(plt = DimPlot(object = acc_cancerCells_noACC1,group.by = "patient.ident"),title = "patient.ident")

patient.ident

print_tab(plt = FeaturePlot(object = acc_cancerCells_noACC1,features = "lum_score"),title = "lum_score")

lum_score

print_tab(plt = FeaturePlot(object = acc_cancerCells_noACC1,features = "myo_score"),title = "myo_score")

myo_score

NA

metagenes by lum_over_myo

acc_cancerCells_noACC1 = subset(acc_cancerCells_noACC1,subset = patient.ident!= "ACC1")

lumScore_vs_program = FetchData(object = acc_cancerCells_noACC1,vars = c("lum_over_myo","program.assignment"))
lumScore_vs_program$program.assignment <- factor(lumScore_vs_program$program.assignment, levels = c("metagene.1","metagene.2","metagene.4"))

lumScore_vs_program = lumScore_vs_program %>% dplyr::filter(program.assignment %in% c("metagene.1","metagene.2","metagene.4"))
ggboxplot(lumScore_vs_program, x = "program.assignment", y = "lum_over_myo",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("metagene.1","metagene.2"),c("metagene.2","metagene.4")))

Metagenes correlation with myo score

for (metagene in c("metagene.1","metagene.2","metagene.4")) {
   lum_score_vs_programScore = FetchData(object = acc_cancerCells_noACC1,vars = c(metagene,"myo_score"))
  print_tab(plt = 
    ggscatter(lum_score_vs_programScore, x = metagene, y = "myo_score", 
            add = "reg.line", conf.int = TRUE, 
            cor.coef = TRUE, cor.method = "pearson")
  ,title = metagene)
}

metagene.1

geom_smooth() using formula ‘y ~ x’

metagene.2

geom_smooth() using formula ‘y ~ x’

metagene.4

geom_smooth() using formula ‘y ~ x’

NA

NA

Metagenes correlation with lum score

for (metagene in c("metagene.1","metagene.2","metagene.4")) {
   lum_score_vs_programScore = FetchData(object = acc_cancerCells_noACC1,vars = c(metagene,"lum_score"))
  print_tab(plt = 
    ggscatter(lum_score_vs_programScore, x = metagene, y = "lum_score", 
            add = "reg.line", conf.int = TRUE, 
            cor.coef = TRUE, cor.method = "pearson")
  ,title = metagene)
}

metagene.1

geom_smooth() using formula ‘y ~ x’

metagene.2

geom_smooth() using formula ‘y ~ x’

metagene.4

geom_smooth() using formula ‘y ~ x’

NA

NA

Notch score

for (metagene in c("metagene.1","metagene.2","metagene.4")) {
   lum_score_vs_programScore = FetchData(object = acc_cancerCells_noACC1,vars = c(metagene,"notch_score"))
  print_tab(plt = 
    ggscatter(lum_score_vs_programScore, x = metagene, y = "notch_score", 
            add = "reg.line", conf.int = TRUE, 
            cor.coef = TRUE, cor.method = "pearson")
  ,title = metagene)
}

metagene.1

geom_smooth() using formula ‘y ~ x’

metagene.2

geom_smooth() using formula ‘y ~ x’

metagene.4

geom_smooth() using formula ‘y ~ x’

NA

NA

Notch ligands score

for (metagene in c("metagene.1","metagene.2","metagene.4")) {
   lum_score_vs_programScore = FetchData(object = acc_cancerCells_noACC1,vars = c(metagene,"notch_ligands_score_score"))
  print_tab(plt = 
    ggscatter(lum_score_vs_programScore, x = metagene, y = "notch_ligands_score_score", 
            add = "reg.line", conf.int = TRUE, 
            cor.coef = TRUE, cor.method = "pearson")
  ,title = metagene)
}

metagene.1

geom_smooth() using formula ‘y ~ x’

metagene.2

geom_smooth() using formula ‘y ~ x’

metagene.4

geom_smooth() using formula ‘y ~ x’

NA

NA