models_cnmf_analysis_9programs

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)

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

k = 9
density_threshold = 0.1 
cnmf_obj.consensus(k=k, density_threshold=density_threshold,show_clustering=True)
usage_norm, gep_scores, _, _ = cnmf_obj.load_results(K=k, density_threshold=density_threshold)

usage_norm = py$usage_norm
gep_scores = py$gep_scores
gep_scores5_xeno = gep_scores #save 
usage_norm5_xeno = usage_norm #save

#load
usage_norm = usage_norm5_xeno
gep_scores = gep_scores5_xeno

3 NMF usage

4 Programs GSEA

  for (col in seq_along(gep_scores)) {
     ranked_vec = gep_scores[,col] %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE) 
     hyp_obj <- fgsea.wrapper(ranked_vec, genesets)
       print_tab(hyp_dots(hyp_obj),title = paste0("gep",col))
  }

gep1

gep2

gep3

gep4

gep5

gep6

gep7

gep8

gep9

NA

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

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

xeno_5_metagenes = py$usage_by_calc #save

all_metagenes = xeno_5_metagenes #load
colnames(all_metagenes) = c("IFNa","immune_response", "hypoxia","cell_cycle","unknown")

6 programs expression