models_cnmf_analysis_k5

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

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

3 K selection plot

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

k = 5
density_threshold = 0.1 
cnmf_obj.consensus(k=k, density_threshold=density_threshold,show_clustering=True)
usage_norm5, gep_scores5, _, _ = cnmf_obj.load_results(K=k, density_threshold=density_threshold)

usage_norm5 = py$usage_norm5
gep_scores5 = py$gep_scores5

4 NMF usage

5 Programs GSEA

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

gep1

gep2

gep3

gep4

gep5

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]

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

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

7 programs expression

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

FeaturePlot(object = xeno,features = colnames(xeno_5_metagenes),ncol = 3)

NA
NA

8 Programs dotplot

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

9 NMF programs regulation

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:4])

IFNa per patient

IFNa

immune_response per patient

immune_response

hypoxia per patient

hypoxia

cell_cycle per patient

cell_cycle

NA

10 Top program 2 genes expression correlation

top_ot = gep_scores5_xeno [order(gep_scores5_xeno [,2],decreasing = T),2,drop = F]%>% head(200) %>% rownames()

num_of_clusters = 7
annotation = plot_genes_cor(dataset = xeno,hallmark_name = NULL,num_of_clusters = num_of_clusters,geneIds = top_ot)

##   genes expression heatmap {.unnumbered }  

NA

11 program 2 all clusters expression

for (chosen_clusters in 1:num_of_clusters) {
  chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  # print(chosen_genes)
  hyp_obj <- hypeR(chosen_genes, genesets_env, test = "hypergeometric", fdr=1, plotting=F,background = rownames(xeno_5_gep_scores))

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

  
  print_tab(plt = 
              hyp_dots(hyp_obj,size_by = "none",title = paste0("cluster",chosen_clusters))+
              FeaturePlot(object = xeno,features = paste0("cluster",chosen_clusters)),
            title = chosen_clusters)
}

1

2

3

4

5

6

7

NA

12 Correlation of clusters

for (chosen_clusters in 1:num_of_clusters) {
  
  cor_res = cor(xeno$TNFa,xeno[[paste0("cluster",chosen_clusters)]])
print(paste("correlation of TNFa program to", paste0("cluster",chosen_clusters),":", cor_res))

}

[1] "correlation of TNFa program to cluster1 : 0.339424898015558"
[1] "correlation of TNFa program to cluster2 : 0.187545650255459"
[1] "correlation of TNFa program to cluster3 : 0.644457044249512"
[1] "correlation of TNFa program to cluster4 : 0.667805824488597"
[1] "correlation of TNFa program to cluster5 : 0.403323740001404"
[1] "correlation of TNFa program to cluster6 : 0.399837119655965"
[1] "correlation of TNFa program to cluster7 : -0.0285178052065907"

clusters_idents = c("cluster1", "KEGG_OXIDATIVE_PHOSPHORYLATION","HALLMARK_TNFA_SIGNALING_VIA_NFKB","KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION","HALLMARK_P53_PATHWAY","cluster6","cluster7")

metagenes_mean_compare <- function(dataset,time.point_var,prefix = "",patient.ident_var,pre_on = c("OSI","NT"),axis.text.x = 11,test = "t.test", programs = c("Hypoxia","TNFa","Cell_cycle"), with_split = T, without_split = T){
  
  for (metegene in programs) {
    #create data:
    genes_by_tp = FetchData(object = dataset,vars = metegene) %>% rowSums() %>% as.data.frame() #mean expression
    names(genes_by_tp)[1] = "Metagene_mean"
    genes_by_tp = cbind(genes_by_tp,FetchData(object = dataset,vars = c(patient.ident_var,time.point_var))) # add id and time points
    
    
    genes_by_tp_forPlot =  genes_by_tp %>% mutate(!!ensym(patient.ident_var) := paste(prefix,genes_by_tp[,patient.ident_var])) #add "model" before  each model/patient
    fm <- as.formula(paste("Metagene_mean", "~", time.point_var)) #make formula to plot
    
    #plot and split by patient:   
    stat.test = compare_means(formula = fm ,data = genes_by_tp_forPlot,method = test,group.by = patient.ident_var)%>% # Add pairwise comparisons p-value
      dplyr::filter(group1 == pre_on[1] & group2 == pre_on[2])  #filter for pre vs on treatment only
    
    plt = ggboxplot(genes_by_tp_forPlot, x = time.point_var, y = "Metagene_mean", color = time.point_var) + #plot
      stat_pvalue_manual(stat.test, label = "p = {p.adj}",  #add p value
                         y.position = max(genes_by_tp_forPlot$Metagene_mean))+ # set position at the top value
      grids()+  
      ylab(paste(metegene,"mean"))+
      theme(axis.text.x = element_text(size = axis.text.x))+
      ylim(0, max(genes_by_tp_forPlot$Metagene_mean)*1.2) # extend y axis to show p value
    
    plt = facet(plt, facet.by = patient.ident_var) #split by patients
    print_tab(plt = plt,title = c(metegene,"per patient")) 
    
    
    #plot = without split by patient:
    if(without_split){
          stat.test = compare_means(formula = fm ,data = genes_by_tp_forPlot,comparisons = my_comparisons,method = test)%>% 
      dplyr::filter(group1 == pre_on[1] & group2 == pre_on[2]) # Add pairwise comparisons p-value
    
    plt = ggboxplot(genes_by_tp_forPlot, x = time.point_var, y = "Metagene_mean", color = time.point_var) +
      stat_pvalue_manual(stat.test, label = "p = {p.adj}",  #add p value
                         y.position = max(genes_by_tp_forPlot$Metagene_mean))+ # set position at the top value
      grids()+  
      ylab(paste(metegene,"mean"))+
      ylim(0, max(genes_by_tp_forPlot$Metagene_mean)*1.2) # extend y axis to show p value
    
    
    print_tab(plt = plt,title = metegene)
    }

  }
  
  
}

13 program 2 intersected genes

programs_of_cluster = c()
for (chosen_clusters in 1:num_of_clusters) {
  chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  pathway_name = clusters_idents[chosen_clusters]
  if (!startsWith(x = pathway_name,prefix = "cluster")){
      chosen_genes  = (chosen_genes) %>% intersect(genesets[[pathway_name]])
      pathway_name = paste0(pathway_name,"_cluster")
  }
  programs_of_cluster = c(programs_of_cluster,pathway_name)
  print(pathway_name)
  print(chosen_genes)
  cat("\n")
  scoresAndIndices <- getPathwayScores(xeno@assays$RNA@data, chosen_genes)
  xeno=AddMetaData(xeno,scoresAndIndices$pathwayScores,pathway_name)

}

[1] "cluster1"
 [1] "STAC2"     "TACR1"     "TRPM4"     "SCPEP1"    "VIM"       "XBP1"      "GLYATL2"   "RCN3"      "TNIP3"     "PTN"      
[11] "RASSF9"    "PIK3R1"    "SOX4"      "SOX2"      "RUNX3"     "ZNF208"    "CLDN8"     "CACNB4"    "TSPAN8"    "PLEKHS1"  
[21] "PLA2G4A"   "LPAR3"     "MUCL1"     "TLE4"      "ALDH2"     "FSIP2"     "LINC00624" "MTRNR2L3"  "PDE4B"     "A1BG"     
[31] "DEFB1"     "ADGRV1"   

[1] "KEGG_OXIDATIVE_PHOSPHORYLATION_cluster"
[1] "MT-ND3"  "MT-ATP8" "MT-ND1"  "MT-ND2" 

[1] "HALLMARK_TNFA_SIGNALING_VIA_NFKB_cluster"
 [1] "ZFP36"   "JUNB"    "IER2"    "GADD45B" "CXCL2"   "BCL6"    "NR4A2"   "NR4A3"   "SOCS3"   "EGR2"    "NFKBIA"  "IL6"    
[13] "IRF1"    "ETS2"   

[1] "KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION_cluster"
[1] "CD74"     "HLA-DRA"  "HLA-DRB1" "HLA-DMA"  "HLA-DPB1" "HLA-DPA1" "HLA-DRB5"

[1] "HALLMARK_P53_PATHWAY_cluster"
[1] "FOS"     "ZFP36L1" "JUN"     "ATF3"    "INHBB"  

[1] "cluster6"
 [1] "SIX1"       "ANKRD36C"   "CP"         "PIGR"       "CNGA1"      "COLEC12"    "AC092683.1" "CHST9"      "NFIB"      
[10] "PIK3IP1"    "REL"        "LINC00342"  "MEF2C"      "BMP3"       "SV2B"       "SLC38A3"    "WFDC2"      "AP000851.1"

[1] "cluster7"
[1] "CFD"    "CASP14" "KRT13" 

14 program 2 intersected pathway regulation

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 = programs_of_cluster,without_split = F)

cluster1 per patient

KEGG_OXIDATIVE_PHOSPHORYLATION_cluster per patient

HALLMARK_TNFA_SIGNALING_VIA_NFKB_cluster per patient

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION_cluster per patient

HALLMARK_P53_PATHWAY_cluster per patient

cluster6 per patient

cluster7 per patient

NA

15 program 2 significant plot

16 Top program 3 genes expression correlation

top_hypoxia = gep_scores5_xeno [order(gep_scores5_xeno [,3],decreasing = T),2,drop = F]%>% head(200) %>% rownames()

num_of_clusters = 4
annotation = plot_genes_cor(dataset = xeno,hallmark_name = NULL,num_of_clusters = num_of_clusters,geneIds = top_hypoxia)

##   genes expression heatmap {.unnumbered }  

NA

17 program 3 all clusters expression

for (chosen_clusters in 1:num_of_clusters) {
  chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  # print(chosen_genes)
  hyp_obj <- hypeR(chosen_genes, genesets_env, test = "hypergeometric", fdr=1, plotting=F,background = rownames(xeno_5_gep_scores))

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

  
  print_tab(plt = 
              hyp_dots(hyp_obj,size_by = "none",title = paste0("cluster",chosen_clusters))+
              FeaturePlot(object = xeno,features = paste0("cluster",chosen_clusters)),
            title = chosen_clusters)


}

1

2

3

4

NA

18 Correlation of clusters

for (chosen_clusters in 1:num_of_clusters) {
  
  cor_res = cor(xeno$hypoxia,xeno[[paste0("cluster",chosen_clusters)]])
print(paste("correlation of hypoxia program to", paste0("cluster",chosen_clusters),":", cor_res))

}

[1] "correlation of hypoxia program to cluster1 : 0.826481794630301"
[1] "correlation of hypoxia program to cluster2 : 0.533339523425782"
[1] "correlation of hypoxia program to cluster3 : 0.419138131097487"
[1] "correlation of hypoxia program to cluster4 : 0.0551682337077848"

clusters_idents = c("HALLMARK_HYPOXIA", "HIF_targets","cluster3","cluster4")

programs_of_cluster = c()
for (chosen_clusters in 1:num_of_clusters) {
  chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  pathway_name = clusters_idents[chosen_clusters]
  if (!startsWith(x = pathway_name,prefix = "cluster")){
      chosen_genes  = (chosen_genes) %>% intersect(genesets[[pathway_name]])
      pathway_name = paste0(pathway_name,"_cluster")
  }
  programs_of_cluster = c(programs_of_cluster,pathway_name)
  print(pathway_name)
  print(chosen_genes)
  cat("\n")
  scoresAndIndices <- getPathwayScores(xeno@assays$RNA@data, chosen_genes)
  xeno=AddMetaData(xeno,scoresAndIndices$pathwayScores,pathway_name)

}

[1] "HALLMARK_HYPOXIA_cluster"
 [1] "VEGFA"    "NDRG1"    "IGFBP3"   "P4HA1"    "ADM"      "DDIT4"    "HK2"      "STC2"     "HSPA5"    "SERPINE1"
[11] "AKAP12"   "PDGFB"    "STC1"     "CAV1"     "TNFAIP3"  "COL5A1"   "PLAUR"    "SLC2A3"   "TGFBI"   

[1] "HIF_targets_cluster"
[1] "ERO1A"   "PGK1"    "SLC16A3" "ENO1"   

[1] "cluster3"
 [1] "ENO2"       "AL133453.1" "PLIN2"      "CA12"       "BHLHE40"    "OSMR"       "SPAG4"      "ZNF395"     "P4HB"      
[10] "PIM3"       "G0S2"       "HIF1A-AS2"  "ATP1B1"     "KATNBL1"    "SAMD4A"     "LRP1"       "SLC6A8"     "IL1RAP"    
[19] "HIST1H2BC"  "MT1F"       "IER3"       "FLNA"       "FAM69C"     "MKNK2"      "HIST1H2AE"  "ELL2"       "HIST1H1C"  
[28] "CA2"        "LPCAT1"     "SLITRK6"    "BMP2"       "FZD8"       "MIF"        "ATP8B3"     "BIK"        "TTYH3"     
[37] "MEIOB"      "APOL2"      "IFNGR2"     "HIST1H4H"   "SPINK1"     "GRIN2A"    

[1] "cluster4"
 [1] "ADAM8"      "MAF"        "DDIT3"      "NUPR1"      "PDCD1"      "MIR155HG"   "FAM13A"     "ETS2"       "DNAJB9"    
[10] "PKP4-AS1"   "GDF15"      "SNHG12"     "TNS1"       "HERPUD1"    "SEMA5B"     "ARRDC3"     "WFIKKN1"    "CLEC2B"    
[19] "XIST"       "CDH7"       "FLNC"       "HLA-DQB1"   "SEC61G"     "AC068672.2" "GPNMB"      "CP"         "PNPLA5"    
[28] "PROX1"      "CASC15"     "RARRES1"    "NME8"       "FBXO32"     "MIAT"       "PGF"        "DEPP1"      "THEMIS2"   
[37] "SMAD3"      "CDH2"       "SPOCK3"     "SOD2"       "IFFO1"      "HLX"        "SYNJ2"      "CHI3L1"     "AC034213.1"
[46] "HTR3C"      "HMOX1"      "PLTP"      

19 program 3 intersected pathway regulation

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 = programs_of_cluster,without_split = F)

HALLMARK_HYPOXIA_cluster per patient

HIF_targets_cluster per patient

cluster3 per patient

cluster4 per patient

NA

20 Top program 3 genes expression correlation

top_cc = gep_scores5_xeno [order(gep_scores5_xeno [,4],decreasing = T),2,drop = F]%>% head(200) %>% rownames()

num_of_clusters = 4
annotation = plot_genes_cor(dataset = xeno,hallmark_name = NULL,num_of_clusters = num_of_clusters,geneIds = top_cc)

##   genes expression heatmap {.unnumbered }  

NA

21 program 3 all clusters expression

for (chosen_clusters in 1:num_of_clusters) {
  chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  # print(chosen_genes)
  hyp_obj <- hypeR(chosen_genes, genesets_env, test = "hypergeometric", fdr=1, plotting=F,background = rownames(xeno_5_gep_scores))

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

  
  print_tab(plt = 
              hyp_dots(hyp_obj,size_by = "none",title = paste0("cluster",chosen_clusters))+
              FeaturePlot(object = xeno,features = paste0("cluster",chosen_clusters)),
            title = chosen_clusters)
  
}

1

2

3

4

NA