patients_cnmf_analysis

1 Data

lung = readRDS("./Data/lung_cancercells_withTP_onlyPatients.rds")
lung_patients = lung$patient.ident %>% unique() %>% as.character()
lung_patients_filtered = lung_patients[!(lung_patients %in% c("X1055new","X1099"))] # remove patients with less than 100 malignant cells
lung = subset(x = lung,subset = patient.ident %in% lung_patients_filtered)

suffix ="xeno_genes_normalized_0-5sigma_2-7theta"
by_expression = T #calculate metagenes by multiplie expression in genes coef, or by cnmf usage

from cnmf import cNMF
suffix = r.suffix
import pickle
f = open('./Data/cnmf/cnmf_objects/patients_' + suffix + '_cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

2 Functions

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

ℹ SHA-1 hash of file is a183df1c565702ecd8ed338bb2abfb0e13415d8e

source_from_github(repositoy = "cNMF_functions",version = "0.3.88",script_name = "cnmf_function_Harmony.R") 

ℹ SHA-1 hash of file is 3e4ca171702a700577edf6399c789fdf505397c2

3 K selection plot

4 gep scores for all NMF k’s

density_threshold = 0.1
usage_norm, gep_scores3, gep_tpm, topgenes = cnmf_obj.load_results(K=3, density_threshold=density_threshold)
usage_norm, gep_scores4, gep_tpm, topgenes = cnmf_obj.load_results(K=4, density_threshold=density_threshold)
usage_norm, gep_scores5, gep_tpm, topgenes = cnmf_obj.load_results(K=5, density_threshold=density_threshold)
usage_norm, gep_scores6, gep_tpm, topgenes = cnmf_obj.load_results(K=6, density_threshold=density_threshold)
usage_norm, gep_scores7, gep_tpm, topgenes = cnmf_obj.load_results(K=7, density_threshold=density_threshold)
usage_norm, gep_scores8, gep_tpm, topgenes = cnmf_obj.load_results(K=8, density_threshold=density_threshold)
usage_norm, gep_scores9, gep_tpm, topgenes = cnmf_obj.load_results(K=9, density_threshold=density_threshold)

5 Enrichment analysis by top 200 genes of each program

gep_scores3 = py$gep_scores3
gep_scores4 = py$gep_scores4
gep_scores5 = py$gep_scores5
gep_scores6 = py$gep_scores6
gep_scores7 = py$gep_scores7
gep_scores8 = py$gep_scores8
gep_scores9 = py$gep_scores9

all_gep_scores =  list(gep_scores3 = gep_scores3, gep_scores4 = gep_scores4, gep_scores5 = gep_scores5, gep_scores6 = gep_scores6, gep_scores7= gep_scores7, gep_scores8 = gep_scores8, gep_scores9 = gep_scores9)
# canonical_pathways = msigdbr(species = "Homo sapiens", category = "C2") %>% dplyr::filter(gs_subcat != "CGP") %>%  dplyr::distinct(gs_name, gene_symbol) 
for (gep_name in names(all_gep_scores)) {
  gep_scores = all_gep_scores[[gep_name]]
  top_genes_num = 200
  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),top_genes_num) #take top top_genes_num
    res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                      names(gep_scores)[i],silent = T,return_all = T,custom_pathways = NULL )
     
    plt_list[[i]] = res$plt
  }
  print_tab(plt =   ggarrange(plotlist = plt_list),
            title = gep_name)
}

gep_scores3

Warning in grSoftVersion() : unable to load shared object ‘/usr/local/lib/R/modules//R_X11.so’: libXt.so.6: cannot open shared object file: No such file or directory

gep_scores4

gep_scores5

gep_scores6

gep_scores7

gep_scores8

gep_scores9

NA

6 Chosen K

selected_k = 8
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

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

top_genes_num = 200
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),top_genes_num) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    names(gep_scores)[i],silent = T,return_all = T,custom_pathways = NULL )
   
  plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)

7 Correlation of programs

gep_scores = py$gep_scores
cor_res = cor(gep_scores)
breaks <- c(seq(-1,1,by=0.01))
colors <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                 (n = length(colors)-3), "red")

pheatmap(cor_res,color = colors,breaks = breaks)

8 correlation by top 150 combined

all_top= c()
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),150) #take top top_genes_num
all_top = c(all_top,top)
}
gep_scores_top = gep_scores[all_top,]
cor_res = cor(gep_scores_top)
pheatmap(cor_res)

9 Combine similar programs

gep_scores = py$gep_scores
# groups_list = list(c(4,5,3),c(1),c(2),c(6))
# groups_list = list(c(4,5,6),c(1),c(2),c(3),c(7),c(9),c(8))
# groups_list = list(c(4,5,7),c(1),c(2),c(3),c(7))
groups_list = list(c(4,3,6),c(1),c(2),c(5),c(7),c(8))

gep_scores = union_programs(groups_list = groups_list,all_metagenes = gep_scores)

top_genes_num = 200
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),top_genes_num) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                     names(gep_scores)[i],silent = T,return_all = T)
   
  plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)

10 Calculate usage

# get expression with genes in cnmf input
lung = FindVariableFeatures(object = lung,nfeatures = 2000)

Calculating gene variances 0% 10 20 30 40 50 60 70 80 90 100% [—-|—-|—-|—-|—-|—-|—-|—-|—-|—-| **************************************************| Calculating feature variances of standardized and clipped values 0% 10 20 30 40 50 60 70 80 90 100% [—-|—-|—-|—-|—-|—-|—-|—-|—-|—-| **************************************************|

genes = rownames(lung)[rownames(lung) %in% VariableFeatures(object = xeno)[1:2000]]

lung_expression = t(as.matrix(GetAssayData(lung,slot='data'))) 
lung_expression = 2**lung_expression #convert from log2(tpm+1) to tpm
lung_expression = lung_expression-1
lung_expression = lung_expression[,genes] %>% as.data.frame()

all_0_genes = colnames(lung_expression)[colSums(lung_expression==0, na.rm=TRUE)==nrow(lung_expression)] #delete rows that have all 0
genes = genes[!genes %in% all_0_genes]
lung_expression = lung_expression[,!colnames(lung_expression) %in% all_0_genes]
gc()

         used    (Mb) gc trigger    (Mb)   max used    (Mb)

Ncells 11889314 635.0 20979452 1120.5 17004190 908.2 Vcells 1641096627 12520.6 2627951469 20049.7 1997200248 15237.5

def get_usage_from_score(counts,tpm, genes,cnmf_obj,k, sumTo1 = True):
      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 like cnmf
      spectra = cnmf_obj.get_median_spectra(k=k) #get score 
      spectra = spectra[spectra.columns.intersection(genes)] #remove genes not in @genes
      spectra = spectra.T.reindex(norm_counts.to_df().columns).T #reorder spectra genes like norm_counts
      
      usage_by_calc,_,_ = non_negative_factorization(X=norm_counts.X, H = spectra.values, update_H=False,n_components = k,max_iter=1000,init ='random')
      usage_by_calc = pd.DataFrame(usage_by_calc, index=counts.index, columns=spectra.index) #insert to df+add names
      if(sumTo1):
          usage_by_calc = usage_by_calc.div(usage_by_calc.sum(axis=1), axis=0) # sum rows to 1 and assign to main df
      usage_by_calc_sumTo1 = usage_by_calc.div(usage_by_calc.sum(axis=1), axis=0) # sum rows to 1
      reorder = usage_by_calc_sumTo1.sum(axis=0).sort_values(ascending=False) #reorder after sum to 1
      usage_by_calc = usage_by_calc.loc[:, reorder.index]
      return(usage_by_calc)

lung_expression = r.lung_expression
genes = r.genes
gep_scores = r.gep_scores
usage_by_calc = get_usage_from_score(counts=lung_expression,tpm=lung_expression,genes=genes,cnmf_obj=cnmf_obj,k=selected_k,sumTo1=False)

/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/bin/python3.7:7: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/bin/python3.7:8: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.

usage_by_calc2 =usage_by_calc %>% rename(gep4.3.6 = cell_cycle)

Error in `chr_as_locations()`:
! Can't rename columns that don't exist.
✖ Column `cell_cycle` doesn't exist.
Backtrace:
  1. usage_by_calc %>% rename(gep4.3.6 = cell_cycle)
  3. dplyr:::rename.data.frame(., gep4.3.6 = cell_cycle)
  4. tidyselect::eval_rename(expr(c(...)), .data)
  5. tidyselect:::rename_impl(...)
  6. tidyselect:::eval_select_impl(...)
     ...
 18. tidyselect:::reduce_sels(node, data_mask, context_mask, init = init)
 19. tidyselect:::walk_data_tree(new, data_mask, context_mask)
 20. tidyselect:::as_indices_sel_impl(...)
 21. tidyselect:::as_indices_impl(x, vars, call = call, strict = strict)
 22. tidyselect:::chr_as_locations(x, vars, call = call)

11 Usgage UMAP

all_metagenes= usage_by_calc

#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% dplyr::select(i)
  lung = AddMetaData(object = lung,metadata = metage_metadata)
}
FeaturePlot(object = lung,features = colnames(all_metagenes),max.cutoff = 70)

12 Assignment

larger_by = 1.25
lung = program_assignment(dataset = lung,larger_by = larger_by,program_names = colnames(all_metagenes))

p = cell_percentage(dataset = lung,time.point_var = "time.point",by_program = T)
print_tab(plt = p,title = "by program")
p = cell_percentage(dataset = lung,time.point_var = "time.point",by_tp = T,x_order = NULL)
print_tab(plt = p,title = "by timepoint")

12.1 program.assignment

# colors =  rainbow(ncol(all_metagenes))
# fc <- colorRampPalette(c("lightgreen", "darkgreen"))
# greens = fc(4)
# colors[1] = "blue"
# colors[2:3] = greens[1:2]
# colors[4] = "red"
# colors[5:6] = greens[3:4]
# colors = c(colors,"grey")
# DimPlot(lung,group.by = "program.assignment",pt.size = 0.5,cols =colors)


# colors =  rainbow(ncol(all_metagenes))
# colors = c(colors,"grey")
# DimPlot(lung,group.by = "program.assignment",pt.size = 0.5,cols =colors)

DimPlot(lung,group.by = "program.assignment",pt.size = 0.5)
DimPlot(lung,group.by = "patient.ident",pt.size = 0.5)
DimPlot(lung,group.by = "time.point",pt.size = 0.5)

13 Score regulation

metagenes_mean_compare(dataset = lung, time.point_var = "time.point",prefix = "patient",patient.ident_var = "patient.ident",pre_on = c("pre-treatment","on-treatment"),axis.text.x = 8,programs = c("hypoxia_like","interferon_like","cell_cycle"))

hypoxia_like per patient

hypoxia_like

interferon_like per patient

interferon_like

cell_cycle per patient

cell_cycle

NA

 df  = FetchData(object = lung,vars = c("program.assignment","time.point")) %>% 
    mutate(program.assignment = if_else(!(program.assignment %in% "metagene.3"), 
                       "not_hypoxia", 
                       program.assignment)) %>% 
    filter (time.point %in% c("pre-treatment","on-treatment")) %>% 
    droplevels() 
  test = fisher.test(table(df))
    
  library(ggstatsplot)
print(
    ggbarstats(
    df, program.assignment, time.point,
    results.subtitle = FALSE,
    subtitle = paste0(
      "Fisher's exact test", ", p-value = ",
      ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
    )
  )
)

14 per patient fisher test

patients_vector = lung$patient.ident %>% unique()
for (patient_name in patients_vector) {
  df  = FetchData(object = lung,vars = c("program.assignment","patient.ident","time.point")) %>% 
    filter (patient.ident == patient_name) %>% 
    filter (program.assignment %in% c("metagene.1","metagene.2")) %>% 
    filter (time.point %in% c("pre-treatment","on-treatment")) %>% 
    select(-patient.ident) %>% 
    droplevels() 
  test = fisher.test(table(df))
    
  library(ggstatsplot)
print(
    ggbarstats(
    df, program.assignment, time.point,
    results.subtitle = FALSE,
    subtitle = paste0(
      "Fisher's exact test", ", p-value = ",
      ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
    ),title = patient_name
  )
)
}