models_cnmf_analysis_k5_V2

1 Functions

library(stringi)

library(reticulate)
source_from_github(repositoy = "DEG_functions",version = "0.2.53")

ℹ SHA-1 hash of file is 8ed4e9017cead897a7d352226a5413f1b5fc0cb7
Welcome to enrichR
Checking connection ... 
Enrichr ... Connection is Live!
FlyEnrichr ... Connection is available!
WormEnrichr ... Connection is available!
YeastEnrichr ... Connection is available!
FishEnrichr ... Connection is available!

source_from_github(repositoy = "cNMF_functions",version = "0.4.04",script_name = "cnmf_functions_V3.R")

ℹ SHA-1 hash of file is a735d9bfb5bc55fcb203f83a7240b9713fc080d3

source_from_github(repositoy = "sc_general_functions",version = "0.1.34",script_name = "functions.R")

ℹ SHA-1 hash of file is 55839fb7db3fdb3e61d8d423b0b3129bfa1d777b

2 Data

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

DimPlot(xeno,group.by = "orig.ident")+ DimPlot(xeno,group.by = "treatment",shuffle = T)

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)

Warning in knitr::include_graphics(plot_path) :
  It is highly recommended to use relative paths for images. You had absolute paths: "/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"

k = 5
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_norm5_xeno  = py$usage_norm
gep_scores5_xeno = py$gep_scores

gep_scores = gep_scores5_xeno
usage_norm = usage_norm5_xeno

4 NMF usage

  for (i in 1:ncol(usage_norm)) {
    metage_metadata = usage_norm %>% dplyr::select(i)
    xeno = AddMetaData(object = xeno,metadata = metage_metadata,col.name = paste0("gep",i))
  }

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 = xeno,features = paste0("gep",1:ncol(usage_norm)),ncol = 3)

5 Programs GSEA

  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)
       print_tab(hyp_dots(hyp_obj,title = paste("program",col))+ aes(size=nes),title = paste0("gep",col))
  }

gep1

gep2

gep3

gep4

gep5

NA

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

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

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

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

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

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

NA
NA

8 Programs dotplot

DotPlot.2(object = xeno, features =  colnames(all_metagenes),group.by  = 'treatment',scale = T)+
  guides(size = guide_legend(title = "Cells expressing (%)"),color = guide_colorbar(title = "Average Score"))

all_pathways = list()
for (pathway in colnames(all_metagenes)) {
  data = FetchData(object = xeno,vars = c(pathway, "treatment", "orig.ident"))
  all_patients = list()
  for (patient in unique(xeno$orig.ident)) {
    mean1 = data %>% filter(orig.ident == patient, treatment == "NT") %>% pull(1) %>% mean()
    mean2 = data %>% filter(orig.ident == patient, treatment == "OSI") %>% pull(1) %>% mean()
    fc =(mean1+1) / (mean2+1)
    all_patients[[patient]] = fc
  }
  all_pathways[[pathway]] = all_patients
}


mat = as.data.frame(lapply(all_pathways, unlist))
mat = log2(t(mat) %>% as.data.frame())
breaks <- c(seq(-1,1,by=0.1))
colors_for_plot <- colorRampPalette(colors = c("blue", "white", "red"))(n = length(breaks))

pheatmap::pheatmap(mat,color = colors_for_plot,breaks = breaks,display_numbers = T,main = "log2(FC) pre/on")

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

IFNa per patient

immune_response per patient

hypoxia per patient

cell_cycle per patient

NA

10 known genes score

ranked_vec = gep_scores[,1] %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE) 
hyp_obj <- hypeR_fgsea(ranked_vec, genesets)
le_genes_ifna =  hyp_obj$data %>% filter(label == "HALLMARK_INTERFERON_ALPHA_RESPONSE") %>% pull("le") %>% strsplit(",") %>% unlist()

ranked_vec = gep_scores[,2] %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE) 
hyp_obj <- hypeR_fgsea(ranked_vec, genesets)
le_genes_tnfa =  hyp_obj$data %>% filter(label == "HALLMARK_TNFA_SIGNALING_VIA_NFKB") %>% pull("le") %>% strsplit(",") %>% unlist()

ranked_vec = gep_scores[,3] %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE) 
hyp_obj <- hypeR_fgsea(ranked_vec, genesets)
le_genes_hif =  hyp_obj$data %>% filter(label == "HIF_targets") %>% pull("le") %>% strsplit(",") %>% unlist()


gene_list = list(IFNa_genes = genesets$HALLMARK_INTERFERON_ALPHA_RESPONSE, TNFa_genes = genesets$HALLMARK_TNFA_SIGNALING_VIA_NFKB, hif_targets = hif_targets,E2F_genes = genesets$HALLMARK_E2F_TARGETS,gep2_top = gep_scores  %>%  arrange(desc(gep_scores[2])) %>% rownames() %>% head(200),TNFa_le = le_genes_tnfa,hif_le = le_genes_hif,IFNa_le  = le_genes_ifna)
# xeno = ScaleData(object = xeno,features = unlist(gene_list))
for (i in seq_along(gene_list)) {
  genes = gene_list[[i]]
  name = names(gene_list)[i]
  scores = FetchData(object = xeno,vars = genes,slot = "data")  %>%  expm1()  %>%  rowMeans() #use  TPM
    # scores = FetchData(object = xeno,vars = genes,slot = "data")  %>%  rowMeans() #use  TPM

  xeno %<>% AddMetaData(metadata = scores,col.name = name)

}

Warning in FetchData.Seurat(object = xeno, vars = genes, slot = "data") :
  The following requested variables were not found: WARS1
Warning in FetchData.Seurat(object = xeno, vars = genes, slot = "data") :
  The following requested variables were not found: CCN1
Warning in FetchData.Seurat(object = xeno, vars = genes, slot = "data") :
  The following requested variables were not found: AK4P1, BNIP3P1, LDHAP5, AL158201.1, MIR210, NLRP3P1, AL109946.1
Warning in FetchData.Seurat(object = xeno, vars = genes, slot = "data") :
  The following requested variables were not found: H2AX, H2AZ1

cor(xeno$immune_response, xeno$gep2_top)

[1] 0.6360318

cor(xeno$immune_response, xeno$TNFa_genes)

[1] 0.3350869

cor(xeno$immune_response, xeno$TNFa_le)

[1] 0.5579202

cat ("\n")

cor(xeno$hypoxia, xeno$hif_targets)

[1] 0.6395481

cor(xeno$hypoxia, xeno$hif_le)

[1] 0.4968675

p = DotPlot(object = xeno, features =  names(gene_list),group.by  = 'treatment',scale = T)+
  guides(size = guide_legend(title = "Cells expressing (%)"),color = guide_colorbar(title = "Average Score"))
 p+scale_radius()

all_pathways = list()
for (pathway in names(gene_list)) {
  data = FetchData(object = xeno,vars = c(pathway, "treatment", "orig.ident"))
  all_patients = list()
  for (patient in unique(xeno$orig.ident)) {
    mean1 = data %>% filter(orig.ident == patient, treatment == "NT") %>% pull(1) %>% mean()
    mean2 = data %>% filter(orig.ident == patient, treatment == "OSI") %>% pull(1) %>% mean()
    fc =(mean1+1) / (mean2+1)
    all_patients[[patient]] = fc
  }
  all_pathways[[pathway]] = all_patients
}


a = as.data.frame(lapply(all_pathways, unlist))
a = log2(t(a) %>% as.data.frame())
breaks <- c(seq(-1,1,by=0.05))
colors_for_plot <- colorRampPalette(colors = c("blue", "white", "red"))(n = length(breaks)-1); colors_for_plot[20:21] = "white"
pheatmap::pheatmap(a,color = colors_for_plot,breaks = breaks,display_numbers = T,main = "log2(FC) pre/on")

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

##   hif_le per patient {.unnumbered }  

NA

for (model in unique(xeno$orig.ident)) {
  cell_name = FetchData(xeno,vars = "orig.ident") %>% filter(orig.ident == model) %>% rownames()
  data = FetchData(object = xeno,vars = colnames(all_metagenes),cells = cell_name)%>% scale() %>% t() 
  annotation = FetchData(object = xeno,vars = c("treatment","orig.ident"),cells = cell_name) %>% arrange(treatment)
  col = list(treatment = c("NT" = "red", "OSI" = "green", "res" = "blue"),orig.ident = c(model = "grey"));names(col[[2]]) = model
  column_ha = HeatmapAnnotation(df = annotation,col = col) 
  data = data[,rownames(annotation)]
  print(ComplexHeatmap::Heatmap(data,show_column_names = F,show_row_names = T,cluster_rows = F,name = "Z-score expression",use_raster = T,cluster_columns = F,
                          row_names_gp = grid::gpar(fontsize = 10), top_annotation = column_ha))

}

# data = FetchData(object = xeno,vars = colnames(all_metagenes))%>% scale() %>% t() 
# annotation = FetchData(object = xeno,vars = c("treatment","orig.ident")) %>% arrange(treatment)
# column_ha = HeatmapAnnotation(df = annotation) 
# data = data[,rownames(annotation)]
# ComplexHeatmap::Heatmap(data,show_column_names = F,show_row_names = T,cluster_rows = F,name = "Z-score expression",use_raster = T,cluster_columns = F,
#                         row_names_gp = grid::gpar(fontsize = 10), top_annotation = column_ha)

library(ComplexHeatmap)
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),50) #take top top_genes_num
  data = FetchData(object = xeno,vars = top)%>% scale() %>% t() 
  metagene_data = FetchData(object = xeno,vars = colnames(all_metagenes)[i],cells = colnames(xeno),slot = "data")
  metagene_data$names = rownames(metagene_data)
  col_list = list(circlize::colorRamp2(c(0, 1), c("white", "red"))); names(col_list) = colnames(all_metagenes)[i]
  column_ha = HeatmapAnnotation(df = metagene_data,col = col_list)
  print(ComplexHeatmap::Heatmap(data,show_column_names = F,show_row_names = T,cluster_rows = F,name = "Z-score expression",use_raster = T,cluster_columns = F))
    
  
  a = DoHeatmap(object = xeno,features = top,group.by = "treatment",combine = F,cells = colnames(xeno),slot = "data")
  print(a[[1]]+
          geom_xsidetile(data = metagene_data,aes(x = names,y = "immune_response",xfill = `immune_response`,width=5, height=15)))
  


}

11 LE genes programs UMAP

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

12 LE genes 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 = programs_main_pathways %>% unlist() %>% paste0("_le"),without_split = F)

col=2
ranked_vec = gep_scores[, col] %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE)
print (paste("running gep",col))
hyp_obj <-hypeR_fgsea(ranked_vec, genesets_go, up_only = T)

print(hyp_dots(hyp_obj, title = paste("program", col), abrv = 70) + aes(size =nes))
  

13 Top program 2 genes expression correlation

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

num_of_clusters = 7
annotation = plot_genes_cor(dataset = xeno,num_of_clusters = num_of_clusters,geneIds = top_ot,height = 3)

14 program 2 all clusters expression

for (chosen_clusters in 1:num_of_clusters) {
  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(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)
}

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

}

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

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

}

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

18 program 2 significant plot

signf_plot_pre_vs_on<- function(dataset,programs,patient.ident_var,prefix,pre_on,test,time.point_var) {
    final_df = data.frame()
    for (metegene in programs) {
      genes_by_tp = FetchData(object = dataset,vars = metegene) %>% rowSums() %>% as.data.frame() #mean expression
      names(genes_by_tp)[1] = metegene
      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(metegene, "~", 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)%>% 
              dplyr::filter(group1 == pre_on[1] & group2 == pre_on[2])  #filter for pre vs on treatment only
      final_df = rbind(final_df,stat.test)
    }
    return(final_df)
}

undebug(signf_plot_pre_vs_on)
final_df = signf_plot_pre_vs_on(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 )
final_df = reshape2::dcast(final_df, orig.ident  ~.y.,value.var = "p.adj") %>% column_to_rownames("orig.ident")

sig_heatmap(all_patients_result = final_df,title = "ad")

19 Top program 3 genes expression correlation

top_hypoxia = gep_scores [order(gep_scores [,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)

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


}

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

}

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)

}

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

23 Top program 3 genes expression correlation

top_cc = gep_scores [order(gep_scores [,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)

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