harmony_all_analysis

1 Parameters

2 Functions

library(stringi)
library(reticulate)
source_from_github(repositoy = "DEG_functions",version = "0.2.24")
source_from_github(repositoy = "cNMF_functions",version = "0.3.84",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
}

# import python functions:
import types

get_norm_counts  = r.get_norm_counts
code_obj = compile(get_norm_counts, '<string>', 'exec')
get_norm_counts = types.FunctionType(code_obj.co_consts[0], globals())

get_usage_from_score  = r.get_usage_from_score
code_obj = compile(get_usage_from_score, '<string>', 'exec')
get_usage_from_score = types.FunctionType(code_obj.co_consts[0], globals())

compute_tpm  = r.compute_tpm
code_obj = compile(compute_tpm, '<string>', 'exec')
compute_tpm = types.FunctionType(code_obj.co_consts[0], globals())

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

4 Models 2K vargenes

suffix = r.suffix
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()

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

selected_k = 3
density_threshold = 0.1
# cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold)
usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

4.1 programs enrichment

gep_scores = py$gep_scores
gep_tpm = py$gep_tpm
usage_norm= py$usage_norm

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)

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

library(highcharter) 
options(highcharter.theme = hc_theme_smpl(tooltip = list(valueDecimals = 2)))

for (program  in names (gep_scores)) {
  print(
        hchart(
      density(gep_scores[,program]), 
      type = "area", name = program
      )
  )
}

for (program  in names (gep_scores)) {
  print(
    ggplot(gep_scores, aes(gep_scores[,program])) + stat_ecdf(geom = "step")

  )
}

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)

xeno_expression = FetchData(object = xeno_expression,vars = vargenes)

Error in UseMethod(generic = "FetchData", object = object) : 
  no applicable method for 'FetchData' applied to an object of class "data.frame"

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

6 Check if original cNMF score is like the calculated score

# cnmf_genes = py$cnmf_genes %>% as.data.frame()
# a = vargenes %>% as.data.frame()
usage_by_calc = py$usage_by_calc
usage_by_calc = usage_by_calc[,c(3,2,1)]
usage_norm = py$usage_norm
cor(usage_by_calc,usage_norm)

           1          2          3
3  0.9994225 -0.6424758 -0.3558304
2 -0.6333332  0.9995869 -0.4933343
1 -0.3536662 -0.4929504  0.9998665

7 calculate score for Xeno

def compute_tpm(input_counts):
    """
    Default TPM normalization
    """
    tpm = input_counts.copy()
    sc.pp.normalize_per_cell(tpm, counts_per_cell_after=1e6,copy=True)
    return(tpm)

all_metagenes = py$usage_by_calc

8 programs expression

all_metagenes = all_metagenes[,c(3,2,1)]

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

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

umap expression

NA

9 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

10 program assignment

larger_by = 1.5
xeno = program_assignment(dataset = xeno,larger_by = larger_by,program_names = colnames(all_metagenes))

print_tab(plt = 
            DimPlot(xeno,group.by = "program.assignment",cols = c(Hypoxia = "red",TNFa = "green",Cell_cycle = "blue","NA" = "grey"))
          ,title = "program.assignment",subtitle_num = 2)

program.assignment

print_tab(plt = 
              DimPlot(xeno,group.by = "orig.ident")
          ,title = "orig.ident",subtitle_num = 2)

orig.ident

print_tab(plt = 
            DimPlot(xeno,group.by = "treatment")
          ,title = "treatment",subtitle_num = 2)

treatment

p = cell_percentage(dataset = xeno,time.point_var = "treatment",by_program = T,x_order = c("NT","OSI","res"))
print_tab(plt = p,title = "by program",subtitle_num = 2)

by program

p = cell_percentage(dataset = xeno,time.point_var = "treatment",by_tp  = T,x_order =c("Hypoxia","TNFa","Cell_cycle","NA"))
print_tab(plt = p,title = "by time point",subtitle_num = 2)

by time point

NA

print_tab(plt = 
            DimPlot(xeno,group.by = "program.assignment",cols = c(Hypoxia = "red",TNFa = "green",Cell_cycle = "blue","NA" = "grey"))
          ,title = "program.assignment",subtitle_num = 3)
print_tab(plt = 
              DimPlot(xeno,group.by = "orig.ident")
          ,title = "orig.ident",subtitle_num = 3)
print_tab(plt = 
            DimPlot(xeno,group.by = "treatment")
          ,title = "treatment",subtitle_num = 3)

p = cell_percentage(dataset = xeno,time.point_var = "treatment",by_program = T,x_order = c("NT","OSI","res"))
print_tab(plt = p,title = "by program",subtitle_num = 3)

p = cell_percentage(dataset = xeno,time.point_var = "treatment",by_tp  = T,x_order =c("Hypoxia","TNFa","Cell_cycle","NA"))
print_tab(plt = p,title = "by time point",subtitle_num = 3)

11 Patients programs expression

quit
lung = FindVariableFeatures(object = lung,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% [—-|—-|—-|—-|—-|—-|—-|—-|—-|—-| **************************************************|

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

all_metagenes = py$usage_by_calc
all_metagenes = all_metagenes[,c(3,2,1)]

names (all_metagenes) = c("Hypoxia","TNFa","Cell_cycle")
#add each metagene to metadata
for (i  in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% dplyr::select(i)
  lung = AddMetaData(object = lung,metadata = metage_metadata)
}

FeaturePlot(object = lung,features = colnames(all_metagenes))

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

metagenes_mean_compare(dataset = lung,time.point_var = "time.point",prefix = "patient",patient.ident_var = "patient.ident",pre_on = c("pre-treatment","on-treatment"))

11.1 lung program assignment

larger_by = 1.25
lung = program_assignment(dataset = lung,larger_by = larger_by,program_names = colnames(all_metagenes))

print_tab(plt = 
            DimPlot(lung,group.by = "program.assignment",cols = c(Hypoxia = "red",TNFa = "green",Cell_cycle = "blue","NA" = "grey"))
          ,title = "program.assignment",subtitle_num = 3)
print_tab(plt = 
              DimPlot(lung,group.by = "patient.ident")
          ,title = "patient.ident",subtitle_num = 3)
print_tab(plt = 
            DimPlot(lung,group.by = "time.point")
          ,title = "time.point",subtitle_num = 3)

p = cell_percentage(dataset = lung,time.point_var = "time.point",by_program = T,x_order = c("pre-treatment","on-treatment","resistant"))
print_tab(plt = p,title = "by program",subtitle_num = 3)

p = cell_percentage(dataset = lung,time.point_var = "time.point",by_tp  = T,x_order =c("Hypoxia","TNFa","Cell_cycle","NA"))
print_tab(plt = p,title = "by time point",subtitle_num = 3)

  top_genes = gep_scores  %>%  arrange(desc(gep_scores["Hypoxia"])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  expr = xeno_expression[,colnames(xeno_expression) %in% top]
  expr_cor = cor(expr)

  pht1 = pheatmap(expr_cor,show_colnames = F,show_rownames = F, silent = T)
      
  
  num_of_clusters = 4
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)
  
  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 = expr_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 = "Hypoxia")
  

  top_genes = gep_scores  %>%  arrange(desc(gep_scores["TNFa"])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  expr = xeno_expression[,colnames(xeno_expression) %in% top]
  expr_cor = cor(expr)

  pht1 = pheatmap(expr_cor,show_colnames = F,show_rownames = F, silent = T)
      
  
  num_of_clusters = 4
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)
  
  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 = expr_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 = "TNFa")
  

  top_genes = gep_scores  %>%  arrange(desc(gep_scores["Cell_cycle"])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  expr = xeno_expression[,colnames(xeno_expression) %in% top]
  expr_cor = cor(expr)

  pht1 = pheatmap(expr_cor,show_colnames = F,show_rownames = F, silent = T)
      
  
  num_of_clusters = 4
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)
  
  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 = expr_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 = "Cell_cycle")
  

top_genes = gep_scores  %>%  arrange(desc(gep_scores["Hypoxia"])) #sort by score a
hypoxia_genes = head(rownames(top_genes),20) #take top top_genes_num
intersect(cluster_3_genes,hypoxia_genes)

[1] "ERO1A" "PGK1" 

library(ggvenn)
all = list(hypoxia_genes = hypoxia_genes, hif_targets = cluster_3_genes)
ggvenn(
  all
)

12 HIF_targets- Hypoxia correlation

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

p = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
  geom_point()+
  geom_smooth(method=lm)
print_tab(plt = p,title = "points")

points

geom_smooth() using formula ‘y ~ x’

p = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
  # geom_point()+
  geom_smooth(method=lm) +
  stat_density_2d(geom = "point", aes(size = after_stat(density)), n = 20, contour = FALSE)+ stat_cor(method = "pearson", label.x = 20, label.y = 1.1)
print_tab(plt = p,title = "points density")

points density

geom_smooth() using formula ‘y ~ x’

p = ggplot(hif_targets_by_tp, aes(x=hif_targets, y=hypoxia_program)) + 
    geom_smooth(method=lm) +
    geom_point()+
  geom_density_2d(aes(color = ..level..)) +
  scale_color_viridis_c()
print_tab(plt = p,title = "polygon")

polygon

geom_smooth() using formula ‘y ~ x’

NA

hif_targets_score = FetchData(object = lung,vars = c(cluster_3_genes)) %>% rowSums() %>% as.data.frame() #mean expression
lung = AddMetaData(object = lung,metadata = hif_targets_score,col.name = "HIF_targets_score")
FeaturePlot(lung,features = "HIF_targets_score")