models_cnmf_analysis_5programs

1 Functions

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

2 Data

xeno = readRDS("./Data/10x_xeno_1000.Rds")
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)

3 Models 2K vargenes

import pickle
from cnmf import cNMF
f = open('./Data/cnmf/cnmf_objects/models_2Kvargenes_cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

4 K selection plot

plot_path = paste0("/sci/labs/yotamd/lab_share/avishai.wizel/R_projects/EGFR/Data/cnmf/cNMF_models_Varnorm_Harmony_2Kvargenes/cNMF_models_Varnorm_Harmony_2Kvargenes.k_selection.png")
knitr::include_graphics(plot_path)

selected_k = 5
density_threshold = 0.1
# cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold)
usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

xeno_5_gep_scores = py$gep_scores

genesets <- msigdb_download("Homo sapiens",category="H") %>% append( msigdb_download("Homo sapiens",category="C2",subcategory = "CP:KEGG"))
genesets[["HIF_targets"]] = hif_targets
genesets = gsets$new(genesets, name="my genesets", version="v1.0")

ranked_list = list()
  for (col in (xeno_5_gep_scores)) {
   lst = col %>% setNames(rownames(xeno_5_gep_scores)) %>% sort(decreasing = TRUE) 
   ranked_list  %<>% append(list(lst))
  }
  names(ranked_list) = paste0("gep",1:ncol(xeno_5_gep_scores))
  
  
  hyp_obj <- hypeR(ranked_list, genesets, test="kstest", fdr=0.05, plotting=F,background = rownames(xeno_5_gep_scores))
  print(
    hyp_dots(hyp_obj,size_by = "significance",abrv = 100,merge = F)
  )

$gep1

$gep2

$gep3

$gep4

$gep5

xeno = FindVariableFeatures(object = xeno,nfeatures = 2000)
xeno_vargenes = VariableFeatures(object = xeno)

xeno_expression = FetchData(object = xeno,vars = xeno_vargenes,slot='counts')
all_0_genes = colnames(xeno_expression)[colSums(xeno_expression==0, na.rm=TRUE)==nrow(xeno_expression)] #delete rows that have all 0
xeno_vargenes = xeno_vargenes[!xeno_vargenes %in% all_0_genes]

5 calculate score for Xeno

import numpy as np
import scanpy as sc
xeno_expression = r.xeno_expression
xeno_vargenes = r.xeno_vargenes
tpm =  compute_tpm(xeno_expression)
usage_by_calc = get_usage_from_score(counts=xeno_expression,tpm=tpm,genes=xeno_vargenes, cnmf_obj=cnmf_obj,k=5)

/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.

all_metagenes = py$usage_by_calc
colnames(all_metagenes) = c("IFNa","immune_response", "hypoxia","cell_cycle","unknown")

6 programs expression

#add each metagene to metadata
for (i  in 1:ncol(all_metagenes)) {
  metagene_metadata = all_metagenes[,i,drop=F]
  xeno = AddMetaData(object = xeno,metadata = metagene_metadata,col.name = names(all_metagenes)[i])
}

print_tab(plt = FeaturePlot(object = xeno,features = colnames(all_metagenes),ncol = 3),title = "umap expression")

umap expression

NA

DotPlot(object = xeno, features =  colnames(all_metagenes),group.by  = 'treatment')

metagenes_mean_compare(dataset = xeno,time.point_var = "treatment",prefix = "model",patient.ident_var = "orig.ident",pre_on = c("NT","OSI"),test = "wilcox.test",programs = colnames(all_metagenes)[1:3])