bivona_data

1 Functions

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

ℹ SHA-1 hash of file is 29c04ac248b96d9a0ee5b6f3b8745d3df770eb54
Loading required package: facefuns

2 Data

load("./Data/Bivona_scRNAseq/NI04_tumor_seurat_object.RData")
tiss_subset_tumor2 = UpdateSeuratObject(tiss_subset_tumor2)

Validating object structure
Updating object slots

wer

Error: object 'wer' not found

tiss_subset_tumor2

calculate nmf metagenes

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

k = 5
density_threshold = 0.1 
_, gep_scores, _, _ = cnmf_obj.load_results(K=k, density_threshold=density_threshold)

3 Calculate usage by counts before Harmony

# get expression with genes in cnmf input
genes = rownames(py$gep_scores)
genes = genes [genes %in% rownames(tiss_subset_tumor2)]
bivona_expression = t(as.matrix(GetAssayData(tiss_subset_tumor2,slot='data'))) 
bivona_expression = 2**bivona_expression #convert from log2(tpm+1) to tpm
bivona_expression = bivona_expression-1
bivona_expression = bivona_expression[,genes] %>% as.data.frame()

all_0_genes = colnames(bivona_expression)[colSums(bivona_expression==0, na.rm=TRUE)==nrow(bivona_expression)] #delete rows that have all 0
genes = genes[!genes %in% all_0_genes]
bivona_expression = bivona_expression[,!colnames(bivona_expression) %in% all_0_genes]
gc(verbose = F)

reticulate::repl_python()
usage_by_calc = get_usage_from_score(counts=bivona_expression,tpm=bivona_expression,genes=genes,cnmf_obj=cnmf_obj,k=5,sumTo1=True)

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

bivona_5_metagenes = py$usage_by_calc

all_metagenes = bivona_5_metagenes
colnames(all_metagenes) = c("IFNa_nmf","immune_response_nmf", "hypoxia_nmf","cell_cycle_nmf","unknown")

4 programs expression

metagenes_violin_compare.2 = function(dataset,prefix = "",pre_on = c("OSI","NT"),axis.text.x = 11,test = "t.test", programs = c("Hypoxia","TNFa","Cell_cycle"),return_list = F,combine_patients = F){
  plt.lst = list()
  if(combine_patients){
    genes_by_tp = FetchData(object = dataset,vars =  c("analysis",programs)) %>% filter(analysis %in% pre_on)  %>% as.data.frame() #mean expression
    formula <- as.formula( paste("c(", paste(programs, collapse = ","), ")~ analysis ") )
    
    #plot and split by patient:   
    stat.test = compare_means(formula = formula ,data = genes_by_tp,method = test,p.adjust.method = "fdr")%>% # Add pairwise comparisons p-value
      dplyr::filter(group1 == pre_on[1] & group2 == pre_on[2])  #filter for pre vs on analysis only
    
    stat.test$p.format =stat.test$p.adj #modift 0 pvalue to be lowest possible float
    stat.test$p.format[!stat.test$p.format == 0 ] <- paste("=",stat.test$p.format[!stat.test$p.format == 0 ])
    stat.test$p.format[stat.test$p.format == 0 ] <- paste("<",.Machine$double.xmin %>% signif(digits = 3))
    
    
    genes_by_tp = reshape2::melt(genes_by_tp, id.vars = c("analysis"),value.name = "score")
    plt = ggplot(genes_by_tp, aes(x = variable, y = score,fill = analysis)) + geom_split_violin(scale = 'width')+ 
      geom_boxplot(width = 0.25, notch = FALSE, notchwidth = .4, outlier.shape = NA, coef=0)+
      ylim(min(genes_by_tp$score),max(genes_by_tp$score)*1.25)
    plt = plt +stat_pvalue_manual(stat.test, label = "p {p.format}",  #add p value
                                  y.position = max(genes_by_tp$score)*1.08,inherit.aes = F,size = 3.3,x = ".y.") # set position at the top value
    return(plt)
  }
  
  for (metegene in programs) {
    #create data:
    genes_by_tp = FetchData(object = dataset,vars =  c("orig.ident","treatment",metegene)) %>% filter(treatment %in% pre_on)  %>% as.data.frame() #mean expression
    names(genes_by_tp)[3] = "Metagene_mean"
    
    fm <- as.formula(paste("Metagene_mean", "~", "treatment")) #make formula to plot
    
    #plot and split by patient:   
    stat.test = compare_means(formula = fm ,data = genes_by_tp,method = test,group.by = "orig.ident",p.adjust.method = "fdr")%>% # Add pairwise comparisons p-value
      dplyr::filter(group1 == pre_on[1] & group2 == pre_on[2])  #filter for pre vs on treatment only
    
    stat.test$p.format =stat.test$p.adj #modift 0 pvalue to be lowest possible float
    stat.test$p.format[!stat.test$p.format == 0 ] <- paste("=",stat.test$p.format[!stat.test$p.format == 0 ])
    stat.test$p.format[stat.test$p.format == 0 ] <- paste("<",.Machine$double.xmin %>% signif(digits = 3))
    
    plt = ggplot(genes_by_tp, aes(x = orig.ident, y = Metagene_mean,fill = treatment)) + geom_split_violin(scale = 'width')+ylab(metegene)+ 
      geom_boxplot(width = 0.25, notch = FALSE, notchwidth = .4, outlier.shape = NA, coef=0)+
      ylim(min(genes_by_tp$Metagene_mean),max(genes_by_tp$Metagene_mean)*1.25)
    plt = plt +stat_pvalue_manual(stat.test, label = "p {p.format}",  #add p value
                                  y.position = max(genes_by_tp$Metagene_mean)*1.08,x = "orig.ident",inherit.aes = F,size = 3.3) # set position at the top value
    
    plt.lst[[metegene]] = plt
    if (!return_list) {
      print(plt)
    }
  }
  
  
  
  
  if (return_list) {
    return(plt.lst)
  }
}

5 NMF programs

metagenes_violin_compare.2(dataset = tiss_subset_tumor2,prefix = "patient",pre_on = c("naive","grouped_pr"),test = "wilcox.test",programs = colnames(all_metagenes)[1:3], return_list = F,combine_patients = T)

genesets <- msigdb_download("Homo sapiens",category="H") 
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)

for (i in seq_along(gene_list)) {
  genes = gene_list[[i]]
  genes = genes[genes %in% rownames(tiss_subset_tumor2)]
  name = names(gene_list)[i]
  scores = FetchData(object = tiss_subset_tumor2,vars = c(genes))
  scores = scores %>% rowMeans() %>% as.data.frame()
  tiss_subset_tumor2 %<>% AddMetaData(metadata = scores,col.name = name)
  
}

6 Signatures

metagenes_violin_compare.2(dataset = tiss_subset_tumor2,prefix = "patient",pre_on = c("naive","grouped_pr"),test = "wilcox.test",programs = names(gene_list), return_list = F,combine_patients = T)

7 metadata

tiss_subset_tumor2@meta.data

23

[1] 23

123

[1] 123

23

[1] 23