HMSC_cNMF

0.1 Parameters

suffix = ""
data_to_read = ""

0.2 functions

source_from_github(repositoy = "DEG_functions",version = "0.2.24")
source_from_github(repositoy = "HMSC_functions",version = "0.1.12",script_name = "functions.R")
source_from_github(repositoy = "cNMF_functions",version = "0.3.72",script_name = "cnmf_function_Harmony.R")

no_neg <- function(x) {
  x = x + abs(min(x))
  x
}

sum_2_one <- function(x) {
  x =x/sum(x)
  x
}

0.3 Data

0.4 Load object

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)
usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

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

1 Hallmark Enrichment analysis by top 200 genes of each program

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
}
gridExtra::grid.arrange(grobs = plt_list)

2 Hallmark + canonical Enrichment analysis by top 200 genes of each program

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
}
gridExtra::grid.arrange(grobs = plt_list)

3 Expression by scores from cnmf

# 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)
  acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = metage_metadata)
}

print_tab(plt = 
            FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes))
,title = "expression",subtitle_num = 2)

expression

print_tab(
  plt = FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes),max.cutoff = 0.05)
          ,title = "max.cutoff = 0.05",subtitle_num = 2)

max.cutoff = 0.05

NA

4 Expression by scores from cnmf + Z score

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

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

print_tab(plt = 
            FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes))
,title = "expression",subtitle_num = 2)

expression

print_tab(
  plt = FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes),max.cutoff = 0.5)
          ,title = "max.cutoff = 0.5",subtitle_num = 2)

max.cutoff = 0.5

NA

5 Expression by coefficients inversion * expression

all_metagenes_inversed  = expression_inversion(gep_scores = gep_scores %>% as.matrix() ,dataset = acc1_cancer_cells) %>% t() %>% as.data.frame()

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

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

FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes_inversed))

6 Expression by (coefficients inversion sum to 1) * expression

gep_scores_norm = apply(gep_scores, 2, no_neg)
gep_scores_norm = apply(gep_scores_norm, 2, sum_2_one)


all_metagenes_inversed  = expression_inversion(gep_scores = gep_scores_norm ,dataset = acc1_cancer_cells) %>% t() %>% as.data.frame()

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

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

print_tab(plt = 
            FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes))
,title = "expression",subtitle_num = 2)

expression

print_tab(
  plt = FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes),max.cutoff = 750)
          ,title = "max.cutoff = 750",subtitle_num = 2)

max.cutoff = 750

NA

7 Expression by coefficients sum to 1 * expression

gep_scores_norm = apply(gep_scores, 2, no_neg)
gep_scores_norm = apply(gep_scores_norm, 2, sum_2_one) %>% as.data.frame()


all_metagenes_mult  = expression_mult(gep_scores = gep_scores_norm ,dataset = acc1_cancer_cells,top_genes = T)

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

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

print_tab(plt = 
            FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes))
,title = "expression",subtitle_num = 2)

expression

print_tab(
  plt = FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes),max.cutoff = 50)
          ,title = "max.cutoff = 50",subtitle_num = 2)

max.cutoff = 50

NA