patients_cnmf_analysis_V2

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_0-5sigma_2-7theta_100iter_26_9"

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.47")
source_from_github(repositoy = "cNMF_functions",version = "0.4.01",script_name = "cnmf_functions_V3.R")
source_from_github(repositoy = "sc_general_functions",version = "0.1.28",script_name = "functions.R")

genesets <- msigdb_download("Homo sapiens",category="H") %>% append( msigdb_download("Homo sapiens",category="C2",subcategory = "CP:KEGG"))
genesets[["HIF_targets"]] = hif_targets

genesets_go <- msigdb_download("Homo sapiens",category="C5",subcategory = "GO:BP")

3 K selection plot

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

cnmf_obj.consensus(k=3, density_threshold=0.1,show_clustering=True)
cnmf_obj.consensus(k=6, density_threshold=0.1,show_clustering=True)
cnmf_obj.consensus(k=7, density_threshold=0.1,show_clustering=True)
cnmf_obj.consensus(k=8, density_threshold=0.1,show_clustering=True)

4 gep scores for all NMF k’s

density_threshold = 0.1
usage_norm3, gep_scores3, gep_tpm, topgenes = cnmf_obj.load_results(K=3, density_threshold=density_threshold)
usage_norm6, gep_scores6, gep_tpm, topgenes = cnmf_obj.load_results(K=6, density_threshold=density_threshold)
usage_norm7, gep_scores7, gep_tpm, topgenes = cnmf_obj.load_results(K=7, density_threshold=density_threshold)
usage_norm8, gep_scores8, gep_tpm, topgenes = cnmf_obj.load_results(K=8, density_threshold=density_threshold)

5 Enrichment analysis by top 200 genes of each program

gep_scores3 = py$gep_scores3
gep_scores6 = py$gep_scores6
gep_scores7 = py$gep_scores7
gep_scores8 = py$gep_scores8

6 GSEA for every program

  for (col in seq_along(gep_scores8)) {
     ranked_vec = gep_scores8[,col] %>% setNames(rownames(gep_scores8)) %>% sort(decreasing = TRUE) 
     hyp_obj <- hypeR_fgsea(ranked_vec, genesets)
       print_tab(hyp_dots(hyp_obj,title = paste("program",col))+ aes(size=nes),title = paste0("gep",col))
  }

gep1

gep2

gep3

gep4

gep5

Warning in fgsea::fgseaMultilevel(stats = signature, pathways = gsets.obj$genesets, : There were 2 pathways for which P-values were not calculated properly due to unbalanced (positive and negative) gene-level statistic values. For such pathways pval, padj, NES, log2err are set to NA. You can try to increase the value of the argument nPermSimple (for example set it nPermSimple = 10000) ## gep6 {.unnumbered }

gep7

gep8

NA

gep_scores = py$gep_scores8
usage_norm = py$usage_norm8

7 progrmas from NMF

colnames(usage_norm) = paste0("gep",1:8)
#add each metagene to metadata
for (i  in 1:ncol(usage_norm )) {
  metagene_metadata = usage_norm [,i,drop=F]
  lung = AddMetaData(object = lung,metadata = metagene_metadata,col.name = colnames(usage_norm)[i])
}

FeaturePlot(object = lung,features = colnames(usage_norm),ncol = 3)

8 progrmas from NMF

metagenes_mean_compare(dataset = lung,time.point_var = "time.point",prefix = "patient",patient.ident_var = "patient.ident",pre_on = c("pre-treatment","on-treatment"),test = "wilcox.test",programs = colnames(usage_norm)[1:5],without_split = F)

gep1 per patient

gep2 per patient

gep3 per patient

gep4 per patient

gep5 per patient

NA

9 Calculate usage by counts before Harmony

# get expression with genes in cnmf input
genes = rownames(gep_scores)
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(verbose = F)

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

Ncells 13005635 694.6 21808952 1164.8 21808952 1164.8 Vcells 1505803258 11488.4 2197729618 16767.4 1767994410 13488.8

lung_expression = r.lung_expression
genes = r.genes
usage_by_calc = get_usage_from_score(counts=lung_expression,tpm=lung_expression,genes=genes,cnmf_obj=cnmf_obj,k=8,sumTo1=True)

/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

colnames(usage_by_calc) = paste0("gep",1:8)
#add each metagene to metadata
for (i  in 1:ncol(usage_by_calc )) {
  metagene_metadata = usage_by_calc [,i,drop=F]
  lung = AddMetaData(object = lung,metadata = metagene_metadata,col.name = colnames(usage_by_calc)[i])
}

FeaturePlot(object = lung,features = colnames(usage_by_calc),ncol = 3)

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

pht = pheatmap(cor_res,color = colors_for_plot,breaks = breaks)
print_tab(pht,title = "correlation")

groups_list = c(5,6)
patients_geps = union_programs(groups_list = groups_list,all_metagenes = patients_geps)

10 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"),test = "wilcox.test",programs = colnames(usage_by_calc)[1:5],without_split = F)

gep1 per patient

gep2 per patient

gep3 per patient

gep4 per patient

gep5 per patient

NA

programs_main_pathways = list(gep1 = 1:5, gep2 = 1:4,gep3 = 1:3)

11 programs LE genes

programs_main_pathways_names = list()

Warning message: 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

  for (col in seq_along(gep_scores)) {
     ranked_vec = gep_scores[,col] %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE) 
     hyp_obj <- hypeR_fgsea(ranked_vec, genesets)
     programs_main_pathways_names[[col]] =  hyp_obj$data[programs_main_pathways[[col]],"label",drop=T]
     for (pathway_num in programs_main_pathways[[col]]) {
        le_genes =  hyp_obj$data[pathway_num,,drop=F] %>% pull("le") %>% strsplit(",") %>% unlist()
        scoresAndIndices = getPathwayScores(dataMatrix = lung@assays$RNA@data,pathwayGenes = le_genes)
        pathway_name = paste0(hyp_obj$data[pathway_num,"label",drop=F],"_le")
        lung=AddMetaData(lung,scoresAndIndices$pathwayScores,col.name = pathway_name)
        print_tab(FeaturePlot(object = lung,features = pathway_name),title = pathway_name)
     }
  }

KEGG_AUTOIMMUNE_THYROID_DISEASE_le

KEGG_TYPE_I_DIABETES_MELLITUS_le

KEGG_INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION_le

KEGG_SYSTEMIC_LUPUS_ERYTHEMATOSUS_le

KEGG_ALLOGRAFT_REJECTION_le

HALLMARK_TNFA_SIGNALING_VIA_NFKB_le

HALLMARK_P53_PATHWAY_le

HALLMARK_UV_RESPONSE_DN_le

KEGG_TGF_BETA_SIGNALING_PATHWAY_le

HALLMARK_HYPOXIA_le

HALLMARK_GLYCOLYSIS_le

HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION_le

Error in programs_main_pathways[[col]] : subscript out of bounds

12 programs LE genes 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"),test = "wilcox.test",programs = programs_main_pathways_names %>% unlist() %>% paste0("_le"),without_split = F)

KEGG_AUTOIMMUNE_THYROID_DISEASE_le per patient

KEGG_TYPE_I_DIABETES_MELLITUS_le per patient

KEGG_INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION_le per patient

KEGG_SYSTEMIC_LUPUS_ERYTHEMATOSUS_le per patient

KEGG_ALLOGRAFT_REJECTION_le per patient

HALLMARK_TNFA_SIGNALING_VIA_NFKB_le per patient

HALLMARK_P53_PATHWAY_le per patient

HALLMARK_UV_RESPONSE_DN_le per patient

KEGG_TGF_BETA_SIGNALING_PATHWAY_le per patient

HALLMARK_HYPOXIA_le per patient

HALLMARK_GLYCOLYSIS_le per patient

HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION_le per patient

NA

13 Top program 2 genes expression correlation

annotation = plot_genes_cor(dataset = xeno,geneIds = top_ot,height = 2.8)

##   genes expression heatmap {.unnumbered }  

NA

14 cluster program expression

for (chosen_clusters in 1:length(unique(annotation$cluster))) {
  chosen_genes = annotation %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  # print(chosen_genes)
  hyp_obj <- hypeR(chosen_genes, genesets, test = "hypergeometric", fdr=1, plotting=F,background = rownames(patients_geps))

   scoresAndIndices <- getPathwayScores(lung@assays$RNA@data, chosen_genes)
  lung=AddMetaData(lung,scoresAndIndices$pathwayScores,paste0("cluster",chosen_clusters))

  
  print_tab(plt = 
              hyp_dots(hyp_obj,size_by = "none",title = paste0("cluster",chosen_clusters))+
              FeaturePlot(object = lung,features = paste0("cluster",chosen_clusters)),
            title = chosen_clusters)
  
  cor_res = cor(lung$interferon_like,lung[[paste0("cluster",chosen_clusters)]])
# print(paste("correlation of TNFa program to", paste0("cluster",chosen_clusters),":", cor_res))

}

1

2

3

4

5

6

7

NA

15 cluster program 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"),test = "wilcox.test",programs = paste0("cluster",1:length(unique(annotation$cluster))),without_split = F)

cluster1 per patient

cluster2 per patient

cluster3 per patient

cluster4 per patient

cluster5 per patient

cluster6 per patient

cluster7 per patient

NA

16 cluster program enrichment

top_ot = patients_geps[order(patients_geps[,1],decreasing = T),2,drop = F]%>% head(10) %>% rownames()

chosen_genes = top_ot
 scoresAndIndices <- getPathwayScores(lung@assays$RNA@data, chosen_genes)
  lung=AddMetaData(lung,scoresAndIndices$pathwayScores,"immune_gep")

  
  
  cor_res = cor(lung$interferon_like,lung[["immune_gep"]])
print(paste("correlation of TNFa program to", "immune_gep",":", cor_res))

for (chosen_clusters in 1:length(unique(annotation$cluster))) {
  chosen_genes = annotation %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  # print(chosen_genes)
  hyp_obj <- hypeR(chosen_genes, genesets, test = "hypergeometric", fdr=1, plotting=F,background = rownames(patients_geps))

   scoresAndIndices <- getPathwayScores(lung@assays$RNA@data, chosen_genes)
  lung=AddMetaData(lung,scoresAndIndices$pathwayScores,paste0("cluster",chosen_clusters))

  
  print_tab(plt = 
              hyp_dots(hyp_obj,size_by = "none",title = paste0("cluster",chosen_clusters))+
              FeaturePlot(object = lung,features = paste0("cluster",chosen_clusters)),
            title = chosen_clusters)
  
  cor_res = cor(lung$interferon_like,lung[[paste0("cluster",chosen_clusters)]])
print(paste("correlation of TNFa program to", paste0("cluster",chosen_clusters),":", cor_res))

}

##   1 {.unnumbered }  

 

[1] "correlation of TNFa program to cluster1 : 0.155039646876199"
##   2 {.unnumbered }  

 

[1] "correlation of TNFa program to cluster2 : 0.35203084026905"
##   3 {.unnumbered }  

 

[1] "correlation of TNFa program to cluster3 : 0.163348283357048"
##   4 {.unnumbered }  

 

[1] "correlation of TNFa program to cluster4 : -0.0422938672700064"
##   5 {.unnumbered }  

 

[1] "correlation of TNFa program to cluster5 : 0.430856560234176"
##   6 {.unnumbered }  

 

[1] "correlation of TNFa program to cluster6 : -0.107526561142319"
##   7 {.unnumbered }  

 

[1] "correlation of TNFa program to cluster7 : -0.0617622603084714"

mult = as.matrix(t(gep_scores[,1,drop=F])) %*% (lung@assays$RNA@data[rownames(gep_scores),] %>% as.matrix())
names(mult) = colnames(lung)
mult = as.data.frame(t(mult))
 lung=AddMetaData(lung,mult,"immune_gep_mult")
  cor_res = cor(lung$interferon_like,lung[["immune_gep_mult"]])
print(paste("correlation of TNFa program to", "immune_gep_mult",":", cor_res))

metagenes_mean_compare(dataset = lung,time.point_var = "time.point",prefix = "patient",patient.ident_var = "patient.ident",pre_on = c("pre-treatment","on-treatment"),test = "wilcox.test",programs = "immune_gep_mult",without_split = F)

# get expression with genes in cnmf input
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(verbose = F)

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

usage_by_calc = py$usage_by_calc
groups_list = c(4,3,6)
usage_by_calc = union_programs(groups_list = groups_list,all_metagenes = usage_by_calc)
# usage_by_calc = apply(usage_by_calc, MARGIN = 1, sum_2_one) %>% t() %>% as.data.frame()
usage_by_calc =usage_by_calc %>% rename(cell_cycle = gep4.3.6, hypoxia_like = gep2, interferon_like = gep1, TNFa =  gep5, INF2 = gep7)

lung=AddMetaData(lung,usage_by_calc[,1,drop=F],"immune_gep_no_sum2one")
metagenes_mean_compare(dataset = lung,time.point_var = "time.point",prefix = "patient",patient.ident_var = "patient.ident",pre_on = c("pre-treatment","on-treatment"),test = "wilcox.test",programs = "immune_gep_no_sum2one",without_split = F)