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.85",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
}

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)

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

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 = 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%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

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

7 Check if original cNMF score is like the calculated score

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

8 calculate score for Xeno

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.

all_metagenes = py$usage_by_calc

9 programs expression

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)
  # metage_metadata = scale(metage_metadata)
  xeno = AddMetaData(object = xeno,metadata = metage_metadata,col.name = names(all_metagenes)[i])
}

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

umap expression

NA

10 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

hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=xeno@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(xeno@assays$RNA@data[intersect(geneIds,var_features),],2,mean)
xeno=AddMetaData(xeno,score,"GO_MITOTIC_CC")

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

GO_MITOTIC_CC per patient

GO_MITOTIC_CC

NA

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

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

Warning: Use of `cc_scores$GO_MITOTIC_CC` is discouraged. Use `GO_MITOTIC_CC` instead.
Warning: Use of `cc_scores$GO_MITOTIC_CC` is discouraged. Use `GO_MITOTIC_CC` instead.
Warning: Use of `cc_scores$GO_MITOTIC_CC` is discouraged. Use `GO_MITOTIC_CC` instead.
Warning: Use of `cc_scores$GO_MITOTIC_CC` is discouraged. Use `GO_MITOTIC_CC` instead.

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

 df  = FetchData(object = xeno,vars = c("is_cycling","treatment")) %>% 
    filter (treatment %in% c("NT","OSI")) %>% 
    droplevels() 
  test = fisher.test(table(df))
    
  library(ggstatsplot)
print(
    ggbarstats(
    df, is_cycling, treatment,
    results.subtitle = FALSE,
    subtitle = paste0(
      "Fisher's exact test", ", p-value = ",
       round(test$p.value,10))
    )
  )

11 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)
print_tab(plt = 
              DimPlot(xeno,group.by = "orig.ident")
          ,title = "orig.ident",subtitle_num = 2)
print_tab(plt = 
            DimPlot(xeno,group.by = "treatment")
          ,title = "treatment",subtitle_num = 2)

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)

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)

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)

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

12 CC signature

# get genes and plot cor heatmap
hallmark_name = "HALLMARK_G2M_CHECKPOINT"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
geneIds= genesets[[hallmark_name]]@geneIds
hallmars_exp = FetchData(object = xeno,vars = c(geneIds))

Warning in FetchData.Seurat(object = xeno, vars = c(geneIds)) : The following requested variables were not found: AC027237.1, AC091021.1, PTTG3P

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

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

genes expression heatmap

NA

#choose clusters
chosen_clusters = c(1)
print ("chosen_clusters= ", chosen_clusters)

[1] “chosen_clusters=”

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

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

Expression

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

dist

Warning: Use of cc_scores$HALLMARK_G2M_CHECKPOINT is discouraged. Use HALLMARK_G2M_CHECKPOINT instead. Warning: Use of cc_scores$HALLMARK_G2M_CHECKPOINT is discouraged. Use HALLMARK_G2M_CHECKPOINT instead. Warning: Use of cc_scores$HALLMARK_G2M_CHECKPOINT is discouraged. Use HALLMARK_G2M_CHECKPOINT instead. Warning: Use of cc_scores$HALLMARK_G2M_CHECKPOINT is discouraged. Use HALLMARK_G2M_CHECKPOINT instead.

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

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

assignment 2 sd

NA

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

fisher

NA

12.1 correlation to nmf program

cor_res = cor(xeno$Cell_cycle,xeno$HALLMARK_G2M_CHECKPOINT)
print(paste("correlation of cc prgoram to HALLMARK_G2M_CHECKPOINT:", cor_res))

[1] "correlation of cc prgoram to HALLMARK_G2M_CHECKPOINT: 0.877586094823362"

13 TNFa signature

hallmark_name = "HALLMARK_TNFA_SIGNALING_VIA_NFKB"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
geneIds= genesets[[hallmark_name]]@geneIds
hallmars_exp = FetchData(object = xeno,vars = c(geneIds))

Warning in FetchData.Seurat(object = xeno, vars = c(geneIds)) : The following requested variables were not found: CCN1

hallmars_exp = hallmars_exp[,colSums(hallmars_exp[])>0] #remove no expression genes
hallmark_cor = cor(hallmars_exp)
pht1 = pheatmap(mat = hallmark_cor,silent = T)

num_of_clusters = 6
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 = 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")

genes expression heatmap

NA

chosen_genes = annotation[["myannotation"]] %>% dplyr::filter(cluster == 3 | cluster == 1 | cluster == 5) %>% rownames() #take relevant genes
score <- apply(xeno@assays$RNA@data[chosen_genes,],2,mean)
xeno=AddMetaData(xeno,score,hallmark_name)

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

Expression

NA

13.1 correlation to nmf program

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

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

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

dist

NA

14 Signature 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("HALLMARK_G2M_CHECKPOINT","HALLMARK_TNFA_SIGNALING_VIA_NFKB"))

HALLMARK_G2M_CHECKPOINT per patient

HALLMARK_G2M_CHECKPOINT

HALLMARK_TNFA_SIGNALING_VIA_NFKB per patient

HALLMARK_TNFA_SIGNALING_VIA_NFKB

NA

15 HIF_targets- Hypoxia correlation

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_hypoxia2[,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")
}

Warning in FetchData.Seurat(object = xeno, vars = c(genes)) : The following requested variables were not found: AK4P1, BNIP3P1, LDHAP5, AL158201.1, MIR210, NLRP3P1, AL109946.1 geom_smooth() using formula ‘y ~ x’ geom_smooth() using formula ‘y ~ x’ ## geom_bin2d {.unnumbered }

geom_smooth() using formula ‘y ~ x’ geom_smooth() using formula ‘y ~ x’ ## geom_bin2d {.unnumbered }

geom_smooth() using formula ‘y ~ x’ geom_smooth() using formula ‘y ~ x’ ## geom_bin2d {.unnumbered }

NA

16 UMAPS

hif_targets_by_tp = FetchData(object = xeno,vars = c(hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression

Warning in FetchData.Seurat(object = xeno, vars = c(hif_targets)) :
  The following requested variables were not found: AK4P1, BNIP3P1, LDHAP5, AL158201.1, MIR210, NLRP3P1, AL109946.1

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

Warning in FetchData.Seurat(object = xeno, vars = c(hif_targets)) :
  The following requested variables were not found: AK4P1, BNIP3P1, LDHAP5, AL158201.1, MIR210, NLRP3P1, AL109946.1

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

Warning in FetchData.Seurat(object = xeno, vars = c(new_hif_targets)) : The following requested variables were not found: AK4P1, BNIP3P1, LDHAP5, AL158201.1, MIR210, NLRP3P1, AL109946.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)  

geom_smooth() using formula ‘y ~ x’ geom_smooth() using formula ‘y ~ x’

  
  print_tab(plt = p,title = "geom_bin2d")

geom_bin2d

NA

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

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

/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_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")
}

Warning in FetchData.Seurat(object = xeno, vars = c(genes)) : The following requested variables were not found: AK4P1, BNIP3P1, LDHAP5, AL158201.1, MIR210, NLRP3P1, AL109946.1 geom_smooth() using formula ‘y ~ x’ geom_smooth() using formula ‘y ~ x’ ## geom_bin2d {.unnumbered }

geom_smooth() using formula ‘y ~ x’ geom_smooth() using formula ‘y ~ x’ ## geom_bin2d {.unnumbered }

geom_smooth() using formula ‘y ~ x’ geom_smooth() using formula ‘y ~ x’ ## geom_bin2d {.unnumbered }

NA

18 UMAPS

hif_targets_by_tp = FetchData(object = xeno,vars = c(hif_targets)) %>% rowSums() %>% as.data.frame() #mean expression

Warning in FetchData.Seurat(object = xeno, vars = c(hif_targets)) :
  The following requested variables were not found: AK4P1, BNIP3P1, LDHAP5, AL158201.1, MIR210, NLRP3P1, AL109946.1

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

Warning in FetchData.Seurat(object = xeno, vars = c(hif_targets)) :
  The following requested variables were not found: AK4P1, BNIP3P1, LDHAP5, AL158201.1, MIR210, NLRP3P1, AL109946.1

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

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

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