models_cnmf_analysis

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

no_neg <- function(x) {
  x = x + abs(min(x))
  x
}

sum_2_one <- function(x) {
  x =x/sum(x)
  x
}

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)

5 gep scores for all NMF k’s

density_threshold = 0.1
usage_norm3, gep_scores3, _, _ = cnmf_obj.load_results(K=3, density_threshold=density_threshold)
usage_norm4, gep_scores4, _, _ = cnmf_obj.load_results(K=4, density_threshold=density_threshold)
usage_norm5, gep_scores5, _, _ = cnmf_obj.load_results(K=5, density_threshold=density_threshold)

gep_scores3 = py$gep_scores3
gep_scores4 = py$gep_scores4
gep_scores5 = py$gep_scores5

usage_norm3 = py$usage_norm3
usage_norm4 = py$usage_norm4
usage_norm5 = py$usage_norm5

all_usage_norm = list(usage_norm3 = usage_norm3, usage_norm4 = usage_norm4, usage_norm5 = usage_norm5)
for (usage_num in seq_along(all_usage_norm)) {
  usage = all_usage_norm[[usage_num]]
  #add each metagene to metadata
  for (i in 1:ncol(usage)) {
    metage_metadata = usage %>% dplyr::select(i)
    xeno = AddMetaData(object = xeno,metadata = metage_metadata,col.name = paste0("gep",i))
  }
  if (usage_num==1) {
      print_tab(FeaturePlot(object = xeno,features = paste0("gep",1:ncol(usage)),ncol = 2),title = usage_num)

  }else{ print_tab(FeaturePlot(object = xeno,features = paste0("gep",1:ncol(usage)),ncol = 3),title = usage_num)}
}

1

2

3

NA

6 all NMF k’s enrichment

all_gep_scores =  list(gep_scores3 = gep_scores3, gep_scores4 = gep_scores4, gep_scores5 = gep_scores5)

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

for (gep_num in seq_along(all_gep_scores)) {
  gep = all_gep_scores[[gep_num]]
  ranked_list = list()
  for (col in (gep)) {
   lst = col %>% setNames(rownames(gep)) %>% sort(decreasing = TRUE) 
   ranked_list  %<>% append(list(lst))
  }
  names(ranked_list) = paste0("gep",1:ncol(gep))
  
  
  
  hyp_obj <- hypeR(ranked_list, genesets_env, test="kstest", fdr=0.05, plotting=F,background = rownames(gep_scores3))
  print_tab(
    hyp_dots(hyp_obj,size_by = "significance",abrv = 100,merge = T)+
      scale_color_continuous(low = "red", high = "black", guide = guide_colorbar(reverse = TRUE))
    ,title = gep_num)
  
}

1

Scale for ‘colour’ is already present. Adding another scale for ‘colour’, which will replace the existing scale.

2

Scale for ‘colour’ is already present. Adding another scale for ‘colour’, which will replace the existing scale.

3

Scale for ‘colour’ is already present. Adding another scale for ‘colour’, which will replace the existing scale.

NA

7 programs enrichment

gep_scores = py$gep_scores3
usage_norm= py$usage_norm3

# or:
gep_scores = readRDS("/sci/labs/yotamd/lab_share/avishai.wizel/R_projects/EGFR/Data/cnmf/harmony_models_gep_scores.rds")

names (gep_scores) = c("Hypoxia","TNFa","Cell_cycle")
plt_list = list()

for (program  in names (gep_scores)) {
 p = ggplot(gep_scores, aes(x=!!ensym(program))) +
  geom_density()+xlab(program)+
   geom_vline(
    aes(xintercept=sort(gep_scores[,program],TRUE)[200]  ,color="top200"),
          linetype="dashed", size=1)+
   geom_vline(
    aes(xintercept=sort(gep_scores[,program],TRUE)[100]  ,color="top100"),
          linetype="dashed", size=1)+
      geom_vline(
    aes(xintercept=sort(gep_scores[,program],TRUE)[50]  ,color="top50"),
          linetype="dashed", size=1)+
         geom_vline(
    aes(xintercept=sort(gep_scores[,program],TRUE)[150]  ,color="top150"),
          linetype="dashed", size=1)+
   scale_color_manual(name = "top n genes", values = c(top200 = "blue",top100 = "red",top150 = "yellow",top50 = "green"))
   plt_list[[program]] <- p

}
 
ggarrange(plotlist = plt_list)

8 Programs hallmark enrichment

ntop = 150
plt_list = list()
hif_targets_set = data.frame(gs_name = "hif_targets",gene_symbol = hif_targets)

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),ntop) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T,custom_pathways  = hif_targets_set)
   
  plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)

9 program 2 all cp

ntop = 150
i = 2 

  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),ntop) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T,custom_pathways  = canonical_pathways)
   
 print(res$plt) 

NA
NA

10 Program 2 all msigdb

ntop = 150
i = 2 

  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),ntop) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T,custom_pathways  = all_pathways)
   
 print(res$plt) 

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]

11 Test with expr after harmony

import numpy as np

import scanpy as sc

expr_after_harmony = sc.read_h5ad('./Data/cnmf/xeno_Harmony_NoNeg_2Kvargenes.h5ad').to_df()
tpm =  compute_tpm(expr_after_harmony)
cnmf_genes = expr_after_harmony.keys().to_list()
usage_by_calc = get_usage_from_score(counts=expr_after_harmony,tpm=tpm,genes=cnmf_genes,cnmf_obj=cnmf_obj,k=3)

12 Check if original cNMF score is like the calculated score

usage_by_calc = py$usage_by_calc
usage_norm = py$usage_norm3
cor(usage_by_calc,usage_norm)

13 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=3)

/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_usage = py$usage_by_calc

14 programs expression

names (xeno_usage) = c("Hypoxia","TNFa","Cell_cycle")
#add each metagene to metadata
for (i  in 1:ncol(xeno_usage)) {
  metage_metadata = xeno_usage %>% dplyr::select(i)
  xeno = AddMetaData(object = xeno,metadata = metage_metadata,col.name = names(xeno_usage)[i])
}

print_tab(plt = FeaturePlot(object = xeno,features = colnames(xeno_usage)),title = "umap expression")

umap expression

NA

15 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 = c("Hypoxia","TNFa","Cell_cycle"))

Hypoxia per patient

Hypoxia

TNFa per patient

TNFa

Cell_cycle per patient

Cell_cycle

NA

DotPlot(object = xeno, features =  c("Hypoxia","TNFa","Cell_cycle"),group.by  = 'treatment')

16 CC signature

# get genes and plot cor heatmap

plot_genes_cor <- function(dataset,hallmark_name,num_of_clusters,geneIds = NULL) {
if(!is_null(geneIds)){
  geneIds = geneIds
}
else{
  geneIds= genesets[[hallmark_name]]@geneIds
}
hallmars_exp = FetchData(object = dataset,vars = c(geneIds))
hallmars_exp = hallmars_exp[,colSums(hallmars_exp[])>0] #remove no expression genes
hallmark_cor = cor(hallmars_exp)
pht1 = pheatmap(mat = hallmark_cor,silent = T)

# make annotations
num_of_clusters = num_of_clusters
clustering_distance = "euclidean"
myannotation = as.data.frame(cutree(pht1[["tree_row"]], k = num_of_clusters)) #split into k clusters
names(myannotation)[1] = "cluster"
myannotation$cluster = as.factor(myannotation$cluster)
palette1 <-brewer.pal(num_of_clusters, "Paired")
names(palette1) = unique(myannotation$cluster)
ann_colors = list (cluster = palette1)
annotation = list(ann_colors = ann_colors, myannotation = myannotation)

#choose colors
colors <- c(seq(-1,1,by=0.01))
my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
                                                 (n = length(colors)-3), "red")


print_tab(plt = 
              pheatmap(mat = hallmark_cor,annotation_col =  annotation[["myannotation"]], annotation_colors = annotation[["ann_colors"]], clustering_distance_rows = clustering_distance,clustering_distance_cols = clustering_distance,color = my_palette,breaks = colors,show_rownames = F,show_colnames = F)
            ,title = "genes expression heatmap")
return(annotation)

}

hallmark_name = "HALLMARK_G2M_CHECKPOINT"
annotation = plot_genes_cor(dataset = xeno,hallmark_name = hallmark_name,num_of_clusters = 4)

Warning in FetchData.Seurat(object = dataset, vars = c(geneIds)) : The following requested variables were not found: AC027237.1, AC091021.1, PTTG3P ## genes expression heatmap {.unnumbered }

NA

#choose clusters
chosen_clusters = c(1,2)
print (paste("chosen_clusters = ", chosen_clusters %>% toString()))

1 “chosen_clusters = 1, 2”

#UMAP expression of signature
chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster %in% chosen_clusters) %>% rownames() #take relevant genes
score <- apply(xeno@assays$RNA@data[chosen_genes,],2,mean)
xeno=AddMetaData(xeno,score,hallmark_name)
scoresAndIndices <- getPathwayScores(xeno@assays$RNA@data, chosen_genes)
xeno=AddMetaData(xeno,scoresAndIndices$pathwayScores,paste0(hallmark_name,"_sipsic"))
print_tab(FeaturePlot(object = xeno, features = paste0(hallmark_name,"_sipsic")),title = "Expression")

Expression

print_tab(FeaturePlot(object = xeno, features = hallmark_name),title = "Expression")

Expression

NA

#plot signature distribution
cc_scores = FetchData(object = xeno,vars = "HALLMARK_G2M_CHECKPOINT")
plt  =  ggplot(cc_scores, aes(x=HALLMARK_G2M_CHECKPOINT)) +
  geom_density()+
   geom_vline(
    aes(xintercept=mean(cc_scores$HALLMARK_G2M_CHECKPOINT) + sd(cc_scores$HALLMARK_G2M_CHECKPOINT) ,color="1 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(cc_scores$HALLMARK_G2M_CHECKPOINT) + 2*sd(cc_scores$HALLMARK_G2M_CHECKPOINT) ,color="2 SD"),
          linetype="dashed", size=1)

print_tab(plt = plt,title = "dist")

# Plot assignment 

cc_scores = cc_scores %>% mutate(is_cycling = if_else(condition = 
                                          HALLMARK_G2M_CHECKPOINT > mean(HALLMARK_G2M_CHECKPOINT) 
                                          +sd(HALLMARK_G2M_CHECKPOINT),
                                        true = "cycling",
                                        false = "non_cycling"))
xeno = AddMetaData(object = xeno,metadata = cc_scores$is_cycling,col.name = "is_cycling")
print_tab(plt = DimPlot(object = xeno,group.by = "is_cycling") , title = "assignment 1 sd")


# cc_scores = cc_scores %>% mutate(is_cycling = if_else(condition = 
#                                           HALLMARK_G2M_CHECKPOINT > mean(HALLMARK_G2M_CHECKPOINT) + 2*sd(cc_scores$HALLMARK_G2M_CHECKPOINT),
#                                         true = "cycling",
#                                         false = "non_cycling"))
# xeno = AddMetaData(object = xeno,metadata = cc_scores$is_cycling,col.name = "is_cycling")
# print_tab(plt = DimPlot(object = xeno,group.by = "is_cycling") , title = "assignment 2 sd")

 df  = FetchData(object = xeno,vars = c("is_cycling","treatment")) %>% 
    filter (treatment %in% c("NT","OSI")) %>% 
    droplevels() 
  test = fisher.test(table(df))
    
  library(ggstatsplot)

    plt = ggbarstats(
    df, is_cycling, treatment,
    results.subtitle = FALSE,
    subtitle = paste0(
      "Fisher's exact test", ", p-value = ",
       round(test$p.value,13))
    )
  
print_tab(plt = plt,title = "fisher")

16.1 correlation to nmf program

cor_res = cor(xeno$Cell_cycle,xeno[[hallmark_name]])
print(paste("correlation of TNFa program to", hallmark_name,":", cor_res))

[1] "correlation of TNFa program to HALLMARK_G2M_CHECKPOINT : 0.926264346933643"

cor_res = cor(xeno$Cell_cycle,xeno[[paste0(hallmark_name,"_sipsic")]])
print(paste("correlation of Cell_cycle program to", paste0(hallmark_name,"_sipsic"),":", cor_res))

[1] "correlation of Cell_cycle program to HALLMARK_G2M_CHECKPOINT_sipsic : 0.929443708586434"

17 TNFa signature

hallmark_name = "HALLMARK_TNFA_SIGNALING_VIA_NFKB"
annotation = plot_genes_cor(dataset = xeno,hallmark_name = hallmark_name,num_of_clusters = 6)

Warning in FetchData.Seurat(object = dataset, vars = c(geneIds)) : The following requested variables were not found: CCN1 ## genes expression heatmap {.unnumbered }

NA

chosen_clusters = c(1,5)
print (paste("chosen_clusters = ", chosen_clusters %>% toString()))

1 “chosen_clusters = 1, 5”

#UMAP expression of signature
chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster %in% chosen_clusters) %>% rownames() #take relevant genes
score <- apply(xeno@assays$RNA@data[chosen_genes,],2,mean)
xeno=AddMetaData(xeno,score,hallmark_name)
scoresAndIndices <- getPathwayScores(xeno@assays$RNA@data, chosen_genes)
xeno=AddMetaData(xeno,scoresAndIndices$pathwayScores,paste0(hallmark_name,"_sipsic"))
print_tab(FeaturePlot(object = xeno, features = paste0(hallmark_name,"_sipsic")),title = "Expression")

Expression

print_tab(FeaturePlot(object = xeno, features = hallmark_name),title = "Expression")

Expression

NA

17.1 correlation to nmf program

cor_res = cor(xeno$TNFa,xeno[[hallmark_name]])
print(paste("correlation of TNFa program to", hallmark_name,":", cor_res))

[1] "correlation of TNFa program to HALLMARK_TNFA_SIGNALING_VIA_NFKB : 0.471118461766633"

cor_res = cor(xeno$TNFa,xeno[[paste0(hallmark_name,"_sipsic")]])
print(paste("correlation of TNFa program to", paste0(hallmark_name,"_sipsic"),":", cor_res))

[1] "correlation of TNFa program to HALLMARK_TNFA_SIGNALING_VIA_NFKB_sipsic : 0.441336272152405"

cc_scores = FetchData(object = xeno,vars = hallmark_name)

plt  =  ggplot(cc_scores, aes(x=HALLMARK_TNFA_SIGNALING_VIA_NFKB)) +
  geom_density()+
   geom_vline(
    aes(xintercept=mean(HALLMARK_TNFA_SIGNALING_VIA_NFKB) + sd(HALLMARK_TNFA_SIGNALING_VIA_NFKB) ,color="1 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(HALLMARK_TNFA_SIGNALING_VIA_NFKB) + 2*sd(HALLMARK_TNFA_SIGNALING_VIA_NFKB) ,color="2 SD"),
          linetype="dashed", size=1)+
    geom_vline(
    aes(xintercept=mean(HALLMARK_TNFA_SIGNALING_VIA_NFKB) ,color="mean"),
          linetype="dashed", size=1)

print_tab(plt = plt,title = "dist")

18 weighted ks test

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 geneset", version="v1.0")
  
ranked_list1 = gep_scores %>% select(1) %>% as_vector() %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE) 
ranked_list2 = gep_scores %>% select(2) %>% as_vector() %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE) 
ranked_list3 = gep_scores %>% select(3) %>% as_vector() %>% setNames(rownames(gep_scores)) %>% sort(decreasing = TRUE) 



ranked_list = list("gep1" = ranked_list1, "gep2" = ranked_list2,"gep3" = ranked_list3)
hyp_obj <- hypeR(ranked_list, genesets, test="kstest", fdr=0.05, plotting=F,background = rownames(gep_scores))
plots = hyp_dots(hyp_obj,size_by = "none")
print_tab(plt = plots[[1]],title = "GEP 1")

GEP 1

print_tab(plt = plots[[2]],title = "GEP 2")

GEP 2

print_tab(plt = plots[[3]],title = "GEP 3")

GEP 3

NA

19 APC signature

hallmark_name = "KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION"
annotation = plot_genes_cor(dataset = xeno,hallmark_name = hallmark_name,num_of_clusters = 4,canocical = T)

genes expression heatmap

NA

chosen_clusters = c(1,3)
print (paste("chosen_clusters = ", chosen_clusters %>% toString()))

1 “chosen_clusters = 1, 3”

#UMAP expression of signature
chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster %in% chosen_clusters) %>% rownames() #take relevant genes
score <- apply(xeno@assays$RNA@data[chosen_genes,],2,mean)
xeno=AddMetaData(xeno,score,hallmark_name)
scoresAndIndices <- getPathwayScores(xeno@assays$RNA@data, chosen_genes)
xeno=AddMetaData(xeno,scoresAndIndices$pathwayScores,paste0(hallmark_name,"_sipsic"))
print_tab(FeaturePlot(object = xeno, features = paste0(hallmark_name,"_sipsic")),title = "Expression")

Expression

print_tab(FeaturePlot(object = xeno, features = hallmark_name),title = "Expression")

Expression

NA

19.1 correlatino to nmf

cor_res = cor(xeno$TNFa,xeno[[hallmark_name]])
print(paste("correlation of TNFa program to", hallmark_name,":", cor_res))

[1] "correlation of TNFa program to KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION : 0.396759488432559"

cor_res = cor(xeno$TNFa,xeno[[paste0(hallmark_name,"_sipsic")]])
print(paste("correlation of TNFa program to", paste0(hallmark_name,"_sipsic"),":", cor_res))

[1] "correlation of TNFa program to KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION_sipsic : 0.333142851072643"

20 Signatures regulation

metagenes_mean_compare(dataset = xeno,time.point_var = "treatment",prefix = "model",patient.ident_var = "orig.ident",pre_on = c("NT","OSI","res"),programs = c("TNFa","HALLMARK_TNFA_SIGNALING_VIA_NFKB","KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION"))

TNFa per patient

TNFa

HALLMARK_TNFA_SIGNALING_VIA_NFKB per patient

HALLMARK_TNFA_SIGNALING_VIA_NFKB

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION per patient

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION

NA

21 intersection

top_ot = gep_scores3[order(gep_scores3[,2],decreasing = T),2,drop = F]%>% head(200) %>% rownames()
hallmark_name = "KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION"
geneIds= genesets[[hallmark_name]]
chosen_genes = intersect(geneIds,top_ot)
print(chosen_genes)

1 “CD74” “HLA-DMA” “HLA-DPA1” “HLA-DPB1” “HLA-DRA” “HLA-DRB1” “HLA-DRB5”

score <- apply(xeno@assays$RNA@data[chosen_genes,],2,mean)
xeno=AddMetaData(xeno,score,paste0(hallmark_name,"intersected"))
scoresAndIndices <- getPathwayScores(xeno@assays$RNA@data, chosen_genes)
xeno=AddMetaData(xeno,scoresAndIndices$pathwayScores,paste0(hallmark_name,"_intersected_sipsic"))
print_tab(FeaturePlot(object = xeno, features = paste0(hallmark_name,"_intersected_sipsic")),title = "Expression")

Expression

Warning 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

print_tab(FeaturePlot(object = xeno, features = paste0(hallmark_name,"intersected")),title = "Expression")

Expression

NA

cor_res = cor(xeno$TNFa,xeno[[paste0(hallmark_name,"intersected")]])
print(paste("correlation of TNFa program to", paste0(hallmark_name,"intersected"),":", cor_res))

[1] "correlation of TNFa program to KEGG_ANTIGEN_PROCESSING_AND_PRESENTATIONintersected : 0.536876281844422"

top_ot = gep_scores3[order(gep_scores3[,2],decreasing = T),2,drop = F]%>% head(200) %>% rownames()
hallmark_name = "HALLMARK_TNFA_SIGNALING_VIA_NFKB"
geneIds= genesets[[hallmark_name]]
chosen_genes = intersect(geneIds,top_ot)
print(chosen_genes)

 [1] "ATF3"    "BCL6"    "BTG2"    "CEBPD"   "CXCL2"   "EGR1"    "EGR2"    "FOS"     "FOSB"    "GADD45B" "IER2"    "JUN"    
[13] "JUNB"    "KLF2"    "NR4A1"   "NR4A2"   "NR4A3"   "REL"     "RHOB"    "SAT1"    "SOCS3"   "ZFP36"  

score <- apply(xeno@assays$RNA@data[chosen_genes,],2,mean)
xeno=AddMetaData(xeno,score,paste0(hallmark_name,"intersected"))
scoresAndIndices <- getPathwayScores(xeno@assays$RNA@data, chosen_genes)
xeno=AddMetaData(xeno,scoresAndIndices$pathwayScores,paste0(hallmark_name,"_intersected_sipsic"))
print_tab(FeaturePlot(object = xeno, features = paste0(hallmark_name,"_intersected_sipsic")),title = "Expression")

##   Expression {.unnumbered }  

print_tab(FeaturePlot(object = xeno, features = paste0(hallmark_name,"intersected")),title = "Expression")

##   Expression {.unnumbered }  

NA

cor_res = cor(xeno$TNFa,xeno[[paste0(hallmark_name,"intersected")]])
print(paste("correlation of TNFa program to", paste0(hallmark_name,"intersected"),":", cor_res))

[1] "correlation of TNFa program to HALLMARK_TNFA_SIGNALING_VIA_NFKBintersected : 0.488904752032433"

22 Norm pathways

sum_scores= xeno$HIF_targets+xeno$KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION+xeno$HALLMARK_G2M_CHECKPOINT

new_score = xeno$KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION/ sum_scores
xeno %<>% AddMetaData(metadata = new_score,col.name = "KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION_norm")

new_score = (xeno$HIF_targets/ sum_scores)%>% replace(is.na(.), 0)
xeno %<>% AddMetaData(metadata = new_score,col.name = "HIF_targets_norm")

new_score = xeno$HALLMARK_G2M_CHECKPOINT/ sum_scores
xeno %<>% AddMetaData(metadata = new_score,col.name = "HALLMARK_G2M_norm")
FeaturePlot(object = xeno,features = c("HALLMARK_G2M_norm", "KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION_norm","HIF_targets_norm"))

metagenes_mean_compare(dataset = xeno,time.point_var = "treatment",prefix = "model",patient.ident_var = "orig.ident",pre_on = c("NT","OSI","res"),programs = c("TNFa","HALLMARK_TNFA_SIGNALING_VIA_NFKB","KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION","KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION_norm"))

TNFa per patient

TNFa

HALLMARK_TNFA_SIGNALING_VIA_NFKB per patient

HALLMARK_TNFA_SIGNALING_VIA_NFKB

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION per patient

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION_norm per patient

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION_norm

NA

23 program 2 clustering

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

24 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(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)
  
  cor_res = cor(xeno$TNFa,xeno[[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

25 correlation of all clusters

for (chosen_clusters in 1:num_of_clusters) {
  chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  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.196505816390455"
[1] "correlation of TNFa program to cluster2 : 0.366147384813019"
[1] "correlation of TNFa program to cluster3 : 0.532136383102456"
[1] "correlation of TNFa program to cluster4 : 0.534315691094927"
[1] "correlation of TNFa program to cluster5 : 0.697969009217511"
[1] "correlation of TNFa program to cluster6 : 0.326695267801052"
[1] "correlation of TNFa program to cluster7 : 0.0360547219515362"

26 clusters 4 and 5

chosen_clusters = c(4,5)
  chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster %in% chosen_clusters) %>% rownames() #take relevant genes
  
  hyp_obj <- hypeR(chosen_genes, genesets_env, test = "hypergeometric", fdr=1, plotting=F,background = rownames(gep_scores))

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

  hyp_dots(hyp_obj,size_by = "none")+
    FeaturePlot(object = xeno,features = score_name)

  cor_res = cor(xeno$TNFa,xeno[[score_name]])
print(paste("correlation of TNFa program to", score_name,":", cor_res))

[1] "correlation of TNFa program to cluster_4_5 : 0.754961916533887"

27 Title

library(biomaRt)
ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")

for (chosen_clusters in 1:num_of_clusters) {
  chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == chosen_clusters) %>% rownames() #take relevant genes
  genedesc <- getBM(attributes=c('external_gene_name','wikigene_description',"interpro_description"), filters = 'external_gene_name', values = chosen_genes, mart =ensembl)  
  print_tab(genedesc,title = chosen_clusters)
}

28 Programs regulation

metagenes_mean_compare(dataset = xeno,time.point_var = "treatment",prefix = "model",patient.ident_var = "orig.ident",pre_on = c("NT","OSI","res"),programs = c("TNFa","HALLMARK_TNFA_SIGNALING_VIA_NFKB","KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION","KEGG_ANTIGEN_PROCESSING_AND_PRESENTATIONintersected","cluster_4_5","cluster7"))

TNFa per patient

TNFa

HALLMARK_TNFA_SIGNALING_VIA_NFKB per patient

HALLMARK_TNFA_SIGNALING_VIA_NFKB

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION per patient

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATIONintersected per patient

KEGG_ANTIGEN_PROCESSING_AND_PRESENTATIONintersected

cluster_4_5 per patient

cluster_4_5

cluster7 per patient

cluster7

NA

29 HIF_targets- Hypoxia correlation

# plot correlation for every subset of hif targets
for (genes in list(hif_targets,xeno_cluster_3_genes,xeno_cluster_3_2_genes)) {
  hif_targets_by_tp = FetchData(object = xeno,vars = c(genes)) %>% rowSums() %>% as.data.frame() #mean expression
  hif_targets_by_tp[,2] = xeno$Hypoxia
  
  names(hif_targets_by_tp) = c("hif_targets","hypoxia_program")
  
  
  
  p1 = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
      geom_point()+
    geom_density_2d(aes(color = ..level..)) +
    geom_smooth(method=lm) +
    stat_cor(method = "pearson", label.x = 20, label.y = 1.1)+
    scale_color_viridis_c()
  
  p2 = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
    geom_bin2d() +
    theme_bw()+ scale_fill_gradientn(limits=c(0,1100), breaks=seq(0, 1100, by=200), colours=c("blue","yellow","red"))+ 
    stat_cor(method = "pearson", label.x = 20, label.y = 1.1)+
    geom_smooth(method=lm) 
  
  p = ggarrange(plotlist = list(p1,p2),nrow  = 2)  
  
  print_tab(plt = p,title = "geom_bin2d")
}

30 UMAPS

hif_targets_by_tp = FetchData(object = xeno,vars = c(hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression
hif_targets_by_tp[,2] = xeno$Hypoxia
names(hif_targets_by_tp) = c("hif_targets","hypoxia_program")

high_hif_low_hypoxia_cells = hif_targets_by_tp %>% filter(hif_targets>25 & hypoxia_program < 0.2) %>% rownames()
low_hif_high_hypoxia_cells = hif_targets_by_tp %>% filter(hif_targets<15 & hypoxia_program > 0.6) %>% rownames()

hif_targets_by_tp = FetchData(object = xeno,vars = c(hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression
xeno = AddMetaData(object = xeno, metadata = hif_targets_by_tp,col.name = "HIF_targets_score")
cells_to_highlight =  list(high_hif_low_hypoxia_cells = high_hif_low_hypoxia_cells, low_hif_high_hypoxia_cells = low_hif_high_hypoxia_cells)

DimPlot(object = xeno, cells.highlight = cells_to_highlight, cols.highlight = c("red","blue"), cols = "gray", order = TRUE)
FeaturePlot(object = xeno,features = c( "HIF_targets_score","Hypoxia","Cell_cycle" ))

hif_hypoxia_correlated_cells = colnames(xeno) [!colnames(xeno) %in% high_hif_low_hypoxia_cells & !colnames(xeno) %in% high_hif_low_hypoxia_cells]
hif_targets_score = FetchData(object = xeno,vars = c(hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression
hif_targets_score[,2] = xeno$Hypoxia
names(hif_targets_score) = c("hif_targets","hypoxia_program")
hif_targets_score  = hif_targets_score %>% mutate(type = case_when(hif_targets>25 & hypoxia_program < 0.2 ~ "high_hif_low_hypoxia",
                                                                   hif_targets<15 & hypoxia_program > 0.6 ~ "low_hif_high_hypoxia",
                                                                   hif_targets>25 & hypoxia_program > 0.6 ~ "high_hif_high_hypoxia",
                                                                   hif_targets<15 & hypoxia_program <0.2 ~ "high_hif_high_hypoxia",
                                                                   TRUE ~ "other"))
xeno = AddMetaData(object = xeno,metadata = hif_targets_score[,"type", drop = F],col.name = "score_correlation")
xeno = SetIdent(object = xeno,value = "score_correlation")
DimPlot(object = xeno,group.by = "score_correlation")
markers = FindMarkers(object = xeno, ident.1 = "high_hif_low_hypoxia",ident.2 = "high_hif_high_hypoxia",densify = T)

updeg = markers %>% filter(p_val_adj<0.05 & avg_log2FC>0) %>% rownames()

new_hif_targets = hif_targets[!hif_targets %in% updeg]
hif_targets_by_tp = FetchData(object = xeno,vars = c(new_hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression
  hif_targets_by_tp[,2] = xeno$Hypoxia
  
  names(hif_targets_by_tp) = c("hif_targets","hypoxia_program")
  
  
  
  p1 = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
      geom_point()+
    geom_density_2d(aes(color = ..level..)) +
    geom_smooth(method=lm) +
    stat_cor(method = "pearson", label.x = 20, label.y = 1.1)+
    scale_color_viridis_c()
  
  p2 = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
    geom_bin2d() +
    theme_bw()+ scale_fill_gradientn(limits=c(0,1100), breaks=seq(0, 1100, by=200), colours=c("blue","yellow","red"))+ 
    stat_cor(method = "pearson", label.x = 20, label.y = 1.1)+
    geom_smooth(method=lm) 
  
  p = ggarrange(plotlist = list(p1,p2),nrow  = 2)  
  
  print_tab(plt = p,title = "geom_bin2d")

upreg_hif_targets = hif_targets[hif_targets %in% updeg]
upreg_hif_targets_expr = FetchData(object = xeno,vars = c(upreg_hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression
xeno = AddMetaData(object = xeno, metadata = upreg_hif_targets_expr,col.name = "upreg_hif_targets_score")
FeaturePlot(object = xeno,features = "upreg_hif_targets_score")

31 Calculate usage without cc in sum to 1

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=3,sumTo1=False)

all_metagenes_noSumTo1 = py$usage_by_calc
tnf_and_hypoxia = all_metagenes_noSumTo1[,1:2]
tnf_and_hypoxia = apply(X = tnf_and_hypoxia, MARGIN = 1, sum_2_one) %>% t() %>%  as.data.frame()
tnf_and_hypoxia[is.na(tnf_and_hypoxia)] <- 0 #replace NAN's with 0.

# plot correlation for every subset of hif targets
for (genes in list(hif_targets,xeno_cluster_3_genes,xeno_cluster_3_2_genes)) {
  hif_targets_by_tp = FetchData(object = xeno,vars = c(genes)) %>% rowSums() %>% as.data.frame() #mean expression
  hif_targets_by_tp[,2] = tnf_and_hypoxia[,1]
  # hif_targets_by_tp[,2] = xeno$Hypoxia
  
  names(hif_targets_by_tp) = c("hif_targets","hypoxia_program")
  
  
  
  p1 = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
      geom_point()+
    geom_density_2d(aes(color = ..level..)) +
    geom_smooth(method=lm) +
    stat_cor(method = "pearson", label.x = 20, label.y = 1.1)+
    scale_color_viridis_c()
  
  p2 = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
    geom_bin2d() +
    theme_bw()+ scale_fill_gradientn(limits=c(0,1100), breaks=seq(0, 1100, by=200), colours=c("blue","yellow","red"))+ 
    stat_cor(method = "pearson", label.x = 20, label.y = 1.1)+
    geom_smooth(method=lm) 
  
  p = ggarrange(plotlist = list(p1,p2),nrow  = 2)  
  
  print_tab(plt = p,title = "geom_bin2d")
}

32 UMAPS

hif_targets_by_tp = FetchData(object = xeno,vars = c(hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression
hif_targets_by_tp[,2] = tnf_and_hypoxia[,1]
names(hif_targets_by_tp) = c("hif_targets","hypoxia_program")

high_hif_low_hypoxia_cells = hif_targets_by_tp %>% filter(hif_targets>25 & hypoxia_program < 0.2) %>% rownames()

hif_targets_by_tp = FetchData(object = xeno,vars = c(hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression
xeno = AddMetaData(object = xeno, metadata = hif_targets_by_tp,col.name = "HIF_targets_score")
xeno = AddMetaData(object = xeno, metadata = tnf_and_hypoxia[,1],col.name = "Hypoxia2")

DimPlot(object = xeno, cells.highlight = high_hif_low_hypoxia_cells, cols.highlight = "red", cols = "gray", order = TRUE)
FeaturePlot(object = xeno,features = c( "HIF_targets_score","Hypoxia2","Cell_cycle" ))
FeaturePlot(object = xeno,features = c("Hypoxia2"))
DimPlot(object = xeno,group.by = "orig.ident")

33 Hypoxia raw

xeno = AddMetaData(object = xeno,metadata = all_metagenes_noSumTo1[,1],col.name = "hypoxia_raw")
FeaturePlot(object = xeno,features = "hypoxia_raw") + scale_color_gradientn(colours = rainbow(5), limits = c(0, 3000))