SiPSiC

1 Functions

library(stringi)
library(reticulate)
source_from_github(repositoy = "DEG_functions",version = "0.2.47")

ℹ SHA-1 hash of file is f5bb1cd741d13bded83fe3b6fd43169e29731216

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.28",script_name = "functions.R")

ℹ SHA-1 hash of file is 737683e4e0a82ffa22769bbd4a0842a271fe63fc

2 Data

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)
xeno = readRDS("./Data/10x_xeno_1000.Rds")


genesets <- msigdb_download("Homo sapiens",category="H") 
hif_targets = c("LDHA", "PPFIA4", "PFKFB4", "AK4", "ERO1A", "BNIP3L", "BNIP3",  #https://pubmed.ncbi.nlm.nih.gov/36384128/
"PFKL", "ENO1", "MIR210HG", "TPI1", "HILPDA", "PGK1", "INSIG2", 
"AK4P1", "AC114803.1", "ENO2", "FUT11", "EGLN3", "NDRG1", "BNIP3P1", 
"WDR54", "OVOL1-AS1", "STC2", "DDX41", "C4orf47", "GYS1", "ANKRD37", 
"BICDL2", "P4HA1", "PDK1", "OSMR", "PKM", "LDHAP5", "AL158201.1", 
"PFKP", "MIR210", "NLRP3P1", "GSX1", "AL109946.1", "PGAM1", "SLC16A3", 
"ISM2", "TCAF2", "ARID3A", "KDM3A", "JMJD6", "C4orf3")
genesets[["HIF_targets"]] = hif_targets

3 Patients log2(TPM)

scoresAndIndices_patients_tpm = list()
i=0
fun <- function(x,data) {
    i<<-i+1
    message(paste("pathway "),i)
    all_patients = c()
    for (patient in unique(lung$patient.ident)) {
        patient_data = subset(lung,subset = patient.ident == patient)
        res = getPathwayScores((2^patient_data@assays$RNA@data)-1, unlist(x))[["pathwayScores"]]
        all_patients = c(all_patients,res)
    }
    scoresAndIndices_patients_tpm = append(scoresAndIndices_patients_tpm,all_patients)


}


scoresAndIndices_patients_tpm = lapply(genesets,data=lung, FUN = fun) 

for (pathway_name in names(scoresAndIndices_patients_tpm)) {
  lung %<>% AddMetaData(metadata = unlist(scoresAndIndices_patients_tpm[[pathway_name]]),col.name = pathway_name)

}

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,p.adjust.method = "fdr")%>% 
              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)
}

final_df = signf_plot_pre_vs_on(dataset = lung,time.point_var = "time.point",prefix = "patient",patient.ident_var = "patient.ident",pre_on = c("pre-treatment","on-treatment"),test = "wilcox.test",programs = names(scoresAndIndices) )
final_df = reshape2::dcast(final_df, patient.ident  ~.y.,value.var = "p.adj") %>% column_to_rownames("patient.ident")

sig_heatmap(all_patients_result = t(final_df) %>% as.data.frame(),title = "ad")

all_pathways = list()
for (hallmark in names(scoresAndIndices_patients_tpm)) {
  data = FetchData(object = lung,vars = c(hallmark, "time.point", "patient.ident"))
  all_patients = list()
  for (patient in unique(lung$patient.ident)) {
    mean1 = data %>% filter(patient.ident == patient, time.point == "pre-treatment") %>% pull(1) %>% mean()
    mean2 = data %>% filter(patient.ident == patient, time.point == "on-treatment") %>% pull(1) %>% mean()
    fc = mean1 / mean2
    all_patients[[patient]] = fc
  }
  all_pathways[[hallmark]] = all_patients
}

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

pheatmap(a,color = colors_for_plot,breaks = breaks,display_numbers = T)

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

4 Patients TPM

scoresAndIndices_patients_logtpm = list()
i=0
fun <- function(x,data) {
    i<<-i+1
    message(paste("pathway "),i)
    all_patients = c()
    for (patient in unique(lung$patient.ident)) {
        patient_data = subset(lung,subset = patient.ident == patient)
        res = getPathwayScores((2^patient_data@assays$RNA@data)-1, unlist(x))[["pathwayScores"]]
        all_patients = c(all_patients,res)
    }
    scoresAndIndices_patients_logtpm = append(scoresAndIndices_patients_logtpm,all_patients)


}


scoresAndIndices_patients_logtpm = lapply(genesets,data=lung, FUN = fun) 

for (pathway_name in names(scoresAndIndices_patients_logtpm)) {
  lung %<>% AddMetaData(metadata = unlist(scoresAndIndices_patients_logtpm[[pathway_name]]),col.name = pathway_name)

}

all_pathways = list()
for (hallmark in names(scoresAndIndices_patients_tpm)) {
  data = FetchData(object = lung,vars = c(hallmark, "time.point", "patient.ident"))
  all_patients = list()
  for (patient in unique(lung$patient.ident)) {
    mean1 = data %>% filter(patient.ident == patient, time.point == "pre-treatment") %>% pull(1) %>% mean()
    mean2 = data %>% filter(patient.ident == patient, time.point == "on-treatment") %>% pull(1) %>% mean()
    fc = mean1 / mean2
    all_patients[[patient]] = fc
  }
  all_pathways[[hallmark]] = all_patients
}

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

pheatmap(a,color = colors_for_plot,breaks = breaks,display_numbers = T)

5 Xeno TPM

scoresAndIndices = list()
i=0
fun <- function(x,dataset) {
    i<<-i+1
    message(paste("pathway "),i)
    all_patients = c()
    for (patient in unique(dataset$orig.ident)) {
        patient_data = subset(dataset,subset = orig.ident == patient)
        res = getPathwayScores((2^patient_data@assays$RNA@data)-1, unlist(x))[["pathwayScores"]]
        all_patients = c(all_patients,res)
    }
    scoresAndIndices = append(scoresAndIndices,all_patients)


}


scoresAndIndices = lapply(genesets,dataset=xeno, FUN = fun) 

for (pathway_name in names(scoresAndIndices)) {
  xeno %<>% AddMetaData(metadata = unlist(scoresAndIndices[[pathway_name]]),col.name = pathway_name)
}

all_pathways = list()
for (hallmark in names(scoresAndIndices)) {
  data = FetchData(object = lung,vars = c(hallmark, "time.point", "patient.ident"))
  all_patients = list()
  for (patient in unique(lung$patient.ident)) {
    mean1 = data %>% filter(patient.ident == patient, time.point == "pre-treatment") %>% pull(1) %>% mean()
    mean2 = data %>% filter(patient.ident == patient, time.point == "on-treatment") %>% pull(1) %>% mean()
    fc = mean1 / mean2
    all_patients[[patient]] = fc
  }
  all_pathways[[hallmark]] = all_patients
}

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

pheatmap(a,color = colors_for_plot,breaks = breaks,display_numbers = T)

6 Xeno log2(TPM)

scoresAndIndices = list()
i=0
fun <- function(x,dataset) {
    i<<-i+1
    message(paste("pathway "),i)
    all_patients = c()
    for (patient in unique(dataset$orig.ident)) {
        patient_data = subset(dataset,subset = orig.ident == patient)
        res = getPathwayScores((2^patient_data@assays$RNA@data)-1, unlist(x))[["pathwayScores"]]
        all_patients = c(all_patients,res)
    }
    scoresAndIndices = append(scoresAndIndices,all_patients)


}


scoresAndIndices = lapply(genesets[1:2],dataset=xeno, FUN = fun) 

for (pathway_name in names(scoresAndIndices)) {
  xeno %<>% AddMetaData(metadata = unlist(scoresAndIndices[[pathway_name]]),col.name = pathway_name)
}

all_pathways = list()
for (hallmark in names(scoresAndIndices)) {
  data = FetchData(object = lung,vars = c(hallmark, "time.point", "patient.ident"))
  all_patients = list()
  for (patient in unique(lung$patient.ident)) {
    mean1 = data %>% filter(patient.ident == patient, time.point == "pre-treatment") %>% pull(1) %>% mean()
    mean2 = data %>% filter(patient.ident == patient, time.point == "on-treatment") %>% pull(1) %>% mean()
    fc = mean1 / mean2
    all_patients[[patient]] = fc
  }
  all_pathways[[hallmark]] = all_patients
}

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

pheatmap(a,color = colors_for_plot,breaks = breaks,display_numbers = T)