analysis

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

source_from_github(repositoy = "DEG_functions",version = "0.2.47")

ℹ SHA-1 hash of file is f5bb1cd741d13bded83fe3b6fd43169e29731216

3 ACC all cells UMAP

DimPlot(acc_all)

4 Neuronal signatures

for (neural_name in colnames(neuronal_signatures)) {
  genes = neuronal_signatures[,neural_name,drop=T]
  # 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)
}

Warning in FetchData.Seurat(object = acc_all, vars = genes, slot = “data”) : The following requested variables were not found: SNHG29, NOP53 Warning in FetchData.Seurat(object = acc_all, vars = genes, slot = “data”) : The following requested variables were not found: NOP53 Warning in FetchData.Seurat(object = acc_all, vars = genes, slot = “data”) : The following requested variables were not found: NOP53

plt = VlnPlot(object = acc_all,features = colnames(neuronal_signatures))
plt[[1]] = plt[[1]]+ ylab("TPM")
plt[[4]] = plt[[4]]+ ylab("TPM")
plt[[7]] = plt[[7]]+ ylab("TPM")
plt

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.1 Neuronal signatures expression

for (neural_name in colnames(neuronal_signatures)) {
  genes = neuronal_signatures[,neural_name,drop=T]
  # 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 = 3)
}

Sympathetic_cholinergic_neuron

Sympathetic_noradrenergic_neuron

Peripheral_nervous_system_neuron_

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

Peripheral_neuron

Sensory_neuron

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

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

Parasympathetic_neuron

NA

4.2 WBC

wbc_cells = subset(acc_all, subset = cell.type == "WBC")
wbc_cells <- FindNeighbors(wbc_cells, dims = 1:10,verbose = F) %>%  FindClusters(resolution = 0.5) %>%  RunUMAP(dims = 1:10,verbose = F)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 33
Number of edges: 528

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.5572
Number of communities: 2
Elapsed time: 0 seconds

DimPlot(wbc_cells)

WBC devided into 2 clusters:

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

4.3 WBC Neuronal signatures

4.4 CD79A neuronal_signatures correlation

library(ggpubr)
for (neural_name in colnames(neuronal_signatures)) {
  wbc_cells = subset(acc_all, subset = cell.type == "WBC")
  data = FetchData(object = wbc_cells,vars = c(neural_name,"CD79A"))
  p = ggplot(data, aes_string(x=neural_name, y="CD79A")) + 
           geom_point()+
           geom_smooth(method=lm) +
           stat_cor()
  print_tab(p,title  = neural_name,subtitle_num = 3)
}

Sympathetic_cholinergic_neuron

geom_smooth() using formula ‘y ~ x’

Sympathetic_noradrenergic_neuron

geom_smooth() using formula ‘y ~ x’

Peripheral_nervous_system_neuron_

geom_smooth() using formula ‘y ~ x’

Autonomic_neuron

geom_smooth() using formula ‘y ~ x’

Peripheral_neuron

geom_smooth() using formula ‘y ~ x’

Sensory_neuron

geom_smooth() using formula ‘y ~ x’

Adrenergic_neuron_tabula_sapiens

geom_smooth() using formula ‘y ~ x’

Peripheral_sensory_neuron

geom_smooth() using formula ‘y ~ x’

Parasympathetic_neuron

geom_smooth() using formula ‘y ~ x’

NA

4.5 PTPRC neuronal_signatures correlation

library(ggpubr)
for (neural_name in colnames(neuronal_signatures)) {
  wbc_cells = subset(acc_all, subset = cell.type == "WBC")
  data = FetchData(object = wbc_cells,vars = c(neural_name,"PTPRC"))
  p = ggplot(data, aes_string(x=neural_name, y="PTPRC")) + 
           geom_point()+
           geom_smooth(method=lm) +
           stat_cor()
  print_tab(p,title  = neural_name,subtitle_num = 3)
}

Sympathetic_cholinergic_neuron

geom_smooth() using formula ‘y ~ x’

Sympathetic_noradrenergic_neuron

geom_smooth() using formula ‘y ~ x’

Peripheral_nervous_system_neuron_

geom_smooth() using formula ‘y ~ x’

Autonomic_neuron

geom_smooth() using formula ‘y ~ x’

Peripheral_neuron

geom_smooth() using formula ‘y ~ x’

Sensory_neuron

geom_smooth() using formula ‘y ~ x’

Adrenergic_neuron_tabula_sapiens

geom_smooth() using formula ‘y ~ x’

Peripheral_sensory_neuron

geom_smooth() using formula ‘y ~ x’

Parasympathetic_neuron

geom_smooth() using formula ‘y ~ x’

NA

5 CSF3

Single-cell RNA sequencing reveals intratumoral heterogeneity and potential mechanisms of malignant progression in prostate cancer with perineural invasion

VlnPlot(object = acc_all,features = "CSF3",group.by = "patient.ident")+ ylab("log2 (TPM+1)")

6 Neuronal pathways in ACC cancer cells

6.1 dim reduction

acc_cancer <- RunPCA(acc_cancer, features = VariableFeatures(object = acc_cancer),verbose = F)
ElbowPlot(acc_cancer)
acc_cancer <- FindNeighbors(acc_cancer, dims = 1:10,verbose = F) %>%  FindClusters(resolution = 0.5) %>%  RunUMAP(dims = 1:10,verbose = F)

DimPlot(acc_cancer,group.by = "patient.ident")

for (pathway_name in names(neuronal_pathways)) {
  genes = neuronal_pathways[[pathway_name]]
  pathways_scores = FetchData(object = acc_cancer,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_cancer  %<>% AddMetaData(metadata = pathways_scores$score,col.name = pathway_name)
}

6.2 Scores UMAP

plt = FeaturePlot(object = acc_cancer,features = names(neuronal_pathways))
for (i in 1:(length(plt$patches$plots)+1)) {
  plt[[i]] = plt[[i]] + theme(plot.title = element_text(size=10.5))+ labs(color='TPM')
}
print_tab(plt,title = "UMAP" ,subtitle_num = 3)

UMAP

plt = FeaturePlot(object = acc_cancer,features = names(neuronal_pathways),max.cutoff = 50)
for (i in 1:(length(plt$patches$plots)+1)) {
  plt[[i]] = plt[[i]] + theme(plot.title = element_text(size=10.5))+ labs(color='TPM')
}
print_tab(plt,title = "UMAP max=50" ,subtitle_num = 3)

UMAP max=50

plt = FeaturePlot(object = acc_cancer,features = names(neuronal_pathways),max.cutoff = 250)
for (i in 1:(length(plt$patches$plots)+1)) {
  plt[[i]] = plt[[i]] + theme(plot.title = element_text(size=10.5))+ labs(color='TPM')
}
print_tab(plt,title = "UMAP max=250" ,subtitle_num = 3)

UMAP max=250

plt = FeaturePlot(object = acc_cancer,features = names(neuronal_pathways))
for (i in 1:(length(plt$patches$plots)+1)) {
  plt[[i]] = plt[[i]] + theme(plot.title = element_text(size=10.5))+ labs(color='TPM') + scale_color_gradientn(colours = rainbow(5))
}

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

print_tab(plt,title = "UMAP rainbow" ,subtitle_num = 3)

UMAP rainbow

NA

6.3 Scores violin

 plt = VlnPlot(object = acc_cancer,features = names(neuronal_pathways),group.by = "patient.ident")+theme(plot.title = element_text(size = 3))
 
for (i in 1:(length(plt$patches$plots)+1)) {
  plt[[i]] = plt[[i]] + theme(plot.title = element_text(size=9.5))
  if (i %in% c(1,4,7)) {
    plt[[i]] = plt[[i]]+ylab("TPM")
  }
}
print(plt)

7 ACC2 DEG

acc_cancer_pri = SetIdent(object = acc_cancer_pri,value = "patient.ident")
markers = FindMarkers(object = acc_cancer_pri,ident.1 = "ACC2",features = VariableFeatures(acc_cancer_pri),densify = T)
volcano_plot(de_genes = markers,max_names = 10,title = "",ident1 = "ACC2",ident2 = "PNI patients",show_graph = F,log2fc_cutoff = 1)

7.1 GSEA- canonical pathways

up = up in ACC2

library(hypeR)
genesets <- msigdb_download("Homo sapiens",category="H") %>% append( msigdb_download("Homo sapiens",category="C2",subcategory = "CP"))
all_genes = markers  %>%  arrange(desc(avg_log2FC)) %>% select("avg_log2FC") 
ranked_list   <- setNames(all_genes$avg_log2FC, rownames(all_genes))
hyp_obj <- hypeR_fgsea(signature = ranked_list,genesets =  genesets,up_only = F)
hyp_dots(hyp_obj)

7.2 GSEA- GO BP

genesets <- msigdb_download("Homo sapiens",category="H") %>% append( msigdb_download("Homo sapiens",category="C5",subcategory = "GO:BP"))
all_genes = markers  %>%  arrange(desc(avg_log2FC)) %>% select("avg_log2FC") 
ranked_list   <- setNames(all_genes$avg_log2FC, rownames(all_genes))
hyp_obj <- hypeR_fgsea(signature = ranked_list,genesets =  genesets,up_only = F)
hyp_dots(hyp_obj)

8 Neuronal pathways in ACC primary cancer cells

lum_score = FetchData(acc_cancer_pri,"luminal_over_myo")
lum_score  %<>% mutate (lum_or_myo = case_when(
         luminal_over_myo > 1 ~ "luminal",
         luminal_over_myo < (-1) ~ "myo",
         TRUE ~ "unknown"))
acc_cancer_pri  %<>% AddMetaData(metadata = lum_score$lum_or_myo,col.name = "lum_or_myo")

for (pathway_name in names(neuronal_pathways)) {
  genes = neuronal_pathways[[pathway_name]]
 pathways_scores = FetchData(object = acc_cancer_pri,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_cancer_pri  %<>% AddMetaData(metadata = pathways_scores$score,col.name = pathway_name)
}

plt = FeaturePlot(object = acc_cancer_pri,features =  names(neuronal_pathways))

for (i in 1:(length(plt$patches$plots)+1)) {
  plt[[i]] = plt[[i]] + theme(plot.title = element_text(size=10.5))+ labs(color='TPM') 
}
print(plt)

9 Lum vs Myo

library(ggpubr)
for (pathway in names(neuronal_pathways)) {
  data = FetchData(object = acc_cancer_pri,vars = c("lum_or_myo",pathway)) %>% filter(lum_or_myo != "unknown")
  p = ggboxplot(data, x = "lum_or_myo", y =pathway,
                palette = "jco",
                add = "jitter")+ 
    stat_compare_means(method = "wilcox.test",comparisons = list(c("luminal","myo")))+ stat_summary(fun.data = function(x) data.frame(y=max(x)*1.2, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("TPM")+ggtitle(pathway)
  print_tab(p,title = pathway)
}

GOBP_NOREPINEPHRINE_SECRETION

GOBP_NOREPINEPHRINE_TRANSPORT

GOBP_NORADRENERGIC_NEURON_DEVELOPMENT

GOBP_NORADRENERGIC_NEURON_DIFFERENTIATION

GOBP_AUTONOMIC_NERVOUS_SYSTEM_DEVELOPMENT

GOBP_PARASYMPATHETIC_NERVOUS_SYSTEM_DEVELOPMENT

GOBP_ADRENERGIC_RECEPTOR_SIGNALING_PATHWAY

GOMF_BETA_2_ADRENERGIC_RECEPTOR_BINDING

NA