R Notebook

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

Functions

library(stringi)
library(reticulate)
source_from_github(repositoy = "DEG_functions",version = "0.2.24")
source_from_github(repositoy = "cNMF_functions",version = "0.3.84",script_name = "cnmf_function_Harmony.R") 

K selection plot

selected_k = 8
suffix = paste(suffix,paste0(selected_k,"nmfK"),sep="_")
print(suffix)

[1] "xeno_genes_normalized_0-5sigma_2-7theta_8nmfK"

selected_k = int(r.selected_k)
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)

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)

Enrichment analysis by top 200 genes of each program

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_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_scores5 {.unnumbered }

|————————————————–| |==================================================| ## gep_scores6 {.unnumbered }

gep_scores7

|————————————————–| |==================================================| ## gep_scores8 {.unnumbered }

gep_scores9

NA

Correlation of programs

gep_scores = py$gep_scores
cor_res = cor(gep_scores)
pheatmap(cor_res)

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)

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)

Calculate usage

lung = FindVariableFeatures(object = lung,nfeatures = 2000)
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()

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
if not (gep_scores is None):
  spectra = gep_scores
  print ("a")
else:
  print ("b")
  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
usage_by_calc = usage_by_calc.div(usage_by_calc.sum(axis=1), axis=0) # sum rows to 1 

def get_usage_from_score(counts,tpm, genes,cnmf_obj,k):
    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
    usage_by_calc = usage_by_calc.div(usage_by_calc.sum(axis=1), axis=0) # sum rows to 1 
    reorder = usage_by_calc.sum(axis=0).sort_values(ascending=False)
    usage_by_calc = usage_by_calc.loc[:, reorder.index]
    return(usage_by_calc)

usage_by_calc = get_usage_from_score(counts=lung_expression,tpm=lung_expression,genes=genes,cnmf_obj=cnmf_obj,k=9)

/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_calc = py$usage_by_calc

Usgage UMAP

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

#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)
}

Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(i)` instead of `i` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.

FeaturePlot(object = lung,features = colnames(all_metagenes),max.cutoff = 70)

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

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

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

by timepoint

NA

program.assignment

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 = "metagene.3")

metagene.3 per patient

metagene.3

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)

Registered S3 method overwritten by 'parameters':
  method                         from      
  format.parameters_distribution datawizard
You can cite this package as:
     Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
     Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167

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

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

patients_vector = lung$patient.ident %>% unique()
for (patient_name in patients_vector) {
  patient_data = subset(x = lung, subset = patient.ident == patient_name)
  cell_percentage(dataset = patient_data,time.point_var = "time.point")

}