paper_figures

1 Functions

library(stringi)
source_from_github(repositoy = "DEG_functions",version = "0.2.47")
source_from_github(repositoy = "cNMF_functions",version = "0.4.0",script_name = "cnmf_functions_V3.R")
source_from_github(repositoy = "sc_general_functions",version = "0.1.28",script_name = "functions.R")

2 Data

library("readxl")
acc_all = readRDS(file = "./Data/acc_tpm_nCount_mito_no146_15k_alldata.rds")
acc_cancer_pri = readRDS(file = "./Data/acc_cancer_no146_primaryonly15k_cancercells.rds")
acc_cancer = readRDS(file = "./Data/acc_tpm_nCount_mito_no146_15k_cancercells.rds")


neuronal_signatures <- read_excel("./Data/Neuronal Signatures.xlsx",col_names =F)
neuronal_signatures  %<>%  t() %>% as.data.frame() %>% janitor::row_to_names(1) %>%  filter(!row_number() == 1)
rownames(neuronal_signatures) <- NULL
colnames(neuronal_signatures)   %<>%  gsub(replacement = "", pattern = "_\\d.*") #remove any _numbers
colnames(neuronal_signatures)   %<>%  gsub(replacement = "", pattern = "\\(.*") #rename "(" to the end
colnames(neuronal_signatures)   %<>%  gsub(replacement = "_", pattern = " ") #rename all spaces

neuronal_pathways <- read_excel("./Data/Pathway analysis Gene sets.xlsx",col_names =F)
neuronal_pathways  %<>%  t() %>% as.data.frame() %>% janitor::row_to_names(1) %>%  filter(!row_number() == 1)  %>%  as.list() 
neuronal_pathways = lapply(neuronal_pathways, na.omit)
neuronal_pathways = lapply(neuronal_pathways, as.character)
b2r_genes = c("RAF1", "PAMPK1", "PRKACA", "PIK3CA", "AKT1", "AKT2", "AKT3", "CREB1", "CREBBP", "BCL2", "BAD", "JUN", "FOS",  "MYC")
neuronal_pathways[["b2r_signature"]] = b2r_genes

neuronal_signatures_no_RP = lapply(neuronal_signatures %>% as.list(),  function(x) {x[!startsWith(x = x,prefix = "RP")]})
names(neuronal_signatures_no_RP) =  paste0(names(neuronal_signatures_no_RP),"_noRP_TPM")
for (neural_name in names(neuronal_signatures_no_RP)) {
  genes = neuronal_signatures_no_RP[[neural_name]]
  # Assuming df is your data frame

  pathways_scores = FetchData(object = acc_all,vars = genes,slot = "data") %>% 
        mutate_all(~ 2^.)%>% mutate_all(~ .-1) %>%  #covert log(TPM+1) to TPM
    rowwise() %>% mutate(score = mean(c_across(everything()))) #mean
  acc_all  %<>% AddMetaData(metadata = pathways_scores$score,col.name = neural_name)
}

FeaturePlot(acc_all,features = "Sympathetic_cholinergic_neuron_noRP_TPM")

df = acc_all$Sympathetic_cholinergic_neuron_noRP_TPM %>% as.data.frame()
names(df) = "epi"
ggplot(df, aes(x=epi)) +
    geom_density()

VlnPlot(object = acc_all,features = names(neuronal_signatures_no_RP))

3 Genes expression

count.data <- GetAssayData(object = acc_all[["RNA"]], slot = "data")
count.data <- as.matrix(x = (2**count.data) - 1)
acc_tpm <- SetAssayData(
    object = acc_all,
    slot = "counts",
    new.data = count.data,
    assay = "RNA"
)

for (neural_name in names(neuronal_signatures_no_RP)) {
  genes = neuronal_signatures_no_RP[[neural_name]]
  # Assuming df is your data frame
  print_tab(
    DoHeatmap(object = acc_tpm,features = genes,slot = "counts",combine = T,disp.max = 5000)+ labs(fill = "TPM") 
  ,title = neural_name,subtitle_num = 2)
}

Sympathetic_cholinergic_neuron_noRP_TPM

Sympathetic_noradrenergic_neuron_noRP_TPM

3.1 Peripheral_nervous_system_neuron__noRP_TPM {.unnumbered }

Warning in DoHeatmap(object = acc_tpm, features = genes, slot = “counts”, : The following features were omitted as they were not found in the counts slot for the RNA assay: NOP53, SNHG29

Autonomic_neuron_noRP_TPM

Peripheral_neuron_noRP_TPM

Sensory_neuron_noRP_TPM

Warning in DoHeatmap(object = acc_tpm, features = genes, slot = “counts”, : The following features were omitted as they were not found in the counts slot for the RNA assay: NOP53

Adrenergic_neuron_tabula_sapiens_noRP_TPM

Warning in DoHeatmap(object = acc_tpm, features = genes, slot = “counts”, : The following features were omitted as they were not found in the counts slot for the RNA assay: NOP53

Peripheral_sensory_neuron_noRP_TPM

Parasympathetic_neuron_noRP_TPM

NA

4 Violin plot

VlnPlot(object = wbc_cells,features = c("CD79A","PTPRC","CXCR4"))

high_cholinergic_neuron = list( high_cholinergic_neuron = colnames(acc_all)[acc_all$Sympathetic_cholinergic_neuron_noRP_TPM>1000])
DimPlot(object = acc_all, cells.highlight = high_cholinergic_neuron, cols.highlight = "red", cols = "gray", order = TRUE)

FeaturePlot(object = acc_all,features = names(neuronal_signatures_no_RP))& theme(plot.title = element_text(size=13))

neuronal_signatures_no_RP = lapply(neuronal_signatures %>% as.list(),  function(x) {x[!startsWith(x = x,prefix = "RP")]})
names(neuronal_signatures_no_RP) =  paste0(names(neuronal_signatures_no_RP),"_noRP_logTPM")
for (neural_name in names(neuronal_signatures_no_RP)) {
  genes = neuronal_signatures_no_RP[[neural_name]]
  # Assuming df is your data frame

  pathways_scores = FetchData(object = acc_all,vars = genes,slot = "data") %>% 
        # mutate_all(~ 2^.)%>% mutate_all(~ .-1) %>%  #covert log(TPM+1) to TPM
    rowwise() %>% mutate(score = mean(c_across(everything()))) #mean
  acc_all  %<>% AddMetaData(metadata = pathways_scores$score,col.name = neural_name)
}

FeaturePlot(acc_all,features = "Sympathetic_cholinergic_neuron_noRP_logTPM")
df = acc_all$Sympathetic_cholinergic_neuron_noRP_logTPM %>% as.data.frame()
names(df) = "epi"
ggplot(df, aes(x=epi)) +
    geom_density()

for (genes in (neuronal_signatures)) {
  p = DoHeatmap(object = acc_all,
            features = genes,
            group.by = "cell.type")
  print(p)
}

5 Pathways

for (pathway_name in names(neuronal_pathways)) {
  genes = neuronal_pathways[[pathway_name]]
  pathways_scores = FetchData(object = acc_cancer_pri,vars = genes,slot = "data") %>% 
    rowwise() %>% mutate(score = mean(c_across(everything()))) #mean
  acc_cancer_pri  %<>% AddMetaData(metadata = pathways_scores$score,col.name = pathway_name)
}

FeaturePlot(object = acc_cancer_pri,features = names(neuronal_pathways)[c(2,9)],pt.size = 0.5,slot = "data") & labs(color='log TPM')&
  theme(plot.title = element_text(size=10.5))

FeaturePlot(object = acc_cancer_pri,features = names(neuronal_pathways)[c(1,2,9)],pt.size = 0.5,slot = "data") & labs(color='log TPM')&
  theme(plot.title = element_text(size=10.5))

library(ggrepel)
data = FetchData(object = acc_cancer_pri,vars =  c("b2r_signature","GOBP_NOREPINEPHRINE_TRANSPORT","patient.ident"))
average_data1  = data %>% group_by(patient.ident) %>%
    dplyr::summarize(b2r_mean = mean(b2r_signature, na.rm=TRUE))

average_data2  = data %>% group_by(patient.ident) %>%
    dplyr::summarize(epi_trans_mean = mean(GOBP_NOREPINEPHRINE_TRANSPORT, na.rm=TRUE))
average_all = cbind(average_data1,average_data2)
average_all = average_all[,c(-1)]
ggplot(average_all,
           aes(x = b2r_mean, y = epi_trans_mean, label=patient.ident)) +  geom_smooth(method = lm) +
  geom_point() + stat_cor(method = "pearson")+geom_text_repel()

`geom_smooth()` using formula 'y ~ x'

library(ggplot2)
library(cdata)
library('WVPlots')

colnames(df) = c("patient.ident","epi_sec","epi_trans","noradrenergic_dev", "noradrenergic_dif","auto_nerv_dev","parasym_nerv","ADREN_sig","b2ADRENE","b2r_signature")
PairPlot(df, 
         colnames(df)[2:length(df)], 
         "Anderson's Iris Data -- 3 species", 
         group_var = "patient.ident") &  geom_smooth(method = lm) &stat_cor(method = "pearson",size=3)