analysis

1 Functions

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

ℹ SHA-1 hash of file is f5bb1cd741d13bded83fe3b6fd43169e29731216
Welcome to enrichR
Checking connection ... 
Enrichr ... Connection is Live!
FlyEnrichr ... Connection is available!
WormEnrichr ... Connection is available!
YeastEnrichr ... Connection is available!
FishEnrichr ... Connection is available!

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

ℹ SHA-1 hash of file is 7b39dfd215cf7e1d29ee976ecc6bfa0e3f1532f6
Loading required package: reticulate

source_from_github(repositoy = "sc_general_functions",version = "0.1.28",script_name = "functions.R")

ℹ SHA-1 hash of file is 737683e4e0a82ffa22769bbd4a0842a271fe63fc
Loading required package: RColorBrewer

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)

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

5 Neuronal signatures NO RP 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")
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)
}

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 = names(neuronal_signatures_no_RP))
plt[[1]] = plt[[1]]+ ylab("TPM")
plt[[4]] = plt[[4]]+ ylab("TPM")
plt[[7]] = plt[[7]]+ ylab("TPM")
plt

5.1 Neuronal signatures 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 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

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

5.3 WBC Neuronal signatures

acc_all$cluster_idents = acc_all@active.ident
data = FetchData(object = acc_all,vars = c("cluster_idents")) %>% mutate(cluster_idents = as.character(cluster_idents))
cluster_0_cells = subset(wbc_cells, subset = seurat_clusters == 0) %>% colnames()
cluster_1_cells = subset(wbc_cells, subset = seurat_clusters == 1) %>% colnames()

data[rownames(data) %in% cluster_0_cells,] ="WBC_cluster_0"
data[rownames(data) %in% cluster_1_cells,] ="WBC_cluster_1"
acc_all$cluster_idents = data

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

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

5.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.6 PTPRC neuronal_signatures correlation no RP

library(ggpubr)
for (neural_name in names(neuronal_signatures_no_RP)) {
  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_noRP

geom_smooth() using formula ‘y ~ x’

Sympathetic_noradrenergic_neuron_noRP

geom_smooth() using formula ‘y ~ x’

5.6.1 Peripheral_nervous_system_neuron__noRP {.unnumbered }

geom_smooth() using formula ‘y ~ x’

Autonomic_neuron_noRP

geom_smooth() using formula ‘y ~ x’

Peripheral_neuron_noRP

geom_smooth() using formula ‘y ~ x’

Sensory_neuron_noRP

geom_smooth() using formula ‘y ~ x’

Adrenergic_neuron_tabula_sapiens_noRP

geom_smooth() using formula ‘y ~ x’

Peripheral_sensory_neuron_noRP

geom_smooth() using formula ‘y ~ x’

Parasympathetic_neuron_noRP

geom_smooth() using formula ‘y ~ x’

NA

FeaturePlot(object = wbc_cells,features = c("PTPRC", colnames(neuronal_signatures)[1:5]))

FeaturePlot(object = wbc_cells,features = c("PTPRC", colnames(neuronal_signatures)[6:9]))

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

7 Neuronal pathways in ACC cancer cells

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

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

Number of nodes: 951
Number of edges: 28888

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.8456
Number of communities: 7
Elapsed time: 0 seconds

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

Warning in FetchData.Seurat(object = acc_cancer, vars = genes, slot = "data") :
  The following requested variables were not found: ADORA3
Warning in FetchData.Seurat(object = acc_cancer, vars = genes, slot = "data") :
  The following requested variables were not found: ADORA3
Warning in FetchData.Seurat(object = acc_cancer, vars = genes, slot = "data") :
  The following requested variables were not found: ASCL1, PHOX2B
Warning in FetchData.Seurat(object = acc_cancer, vars = genes, slot = "data") :
  The following requested variables were not found: ASCL1, PHOX2B
Warning in FetchData.Seurat(object = acc_cancer, vars = genes, slot = "data") :
  The following requested variables were not found: ARX, HOXB1, ASCL1, PRLH, PHOX2B, ISX
Warning in FetchData.Seurat(object = acc_cancer, vars = genes, slot = "data") :
  The following requested variables were not found: HOXB1, PHOX2B
Warning in FetchData.Seurat(object = acc_cancer, vars = genes, slot = "data") :
  The following requested variables were not found: DRD3, MIR1-1, MIR133A1

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

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

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

8.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)

$up

$dn

8.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)

Warning in fgsea::fgseaMultilevel(stats = signature, pathways = gsets.obj$genesets,  :
  For some of the pathways the P-values were likely overestimated. For such pathways log2err is set to NA.

hyp_dots(hyp_obj)

$up

$dn

9 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)
}

Warning in FetchData.Seurat(object = acc_cancer_pri, vars = genes, slot = "data") :
  The following requested variables were not found: ADORA3
Warning in FetchData.Seurat(object = acc_cancer_pri, vars = genes, slot = "data") :
  The following requested variables were not found: ADORA3
Warning in FetchData.Seurat(object = acc_cancer_pri, vars = genes, slot = "data") :
  The following requested variables were not found: ASCL1, PHOX2B
Warning in FetchData.Seurat(object = acc_cancer_pri, vars = genes, slot = "data") :
  The following requested variables were not found: ASCL1, PHOX2B
Warning in FetchData.Seurat(object = acc_cancer_pri, vars = genes, slot = "data") :
  The following requested variables were not found: ARX, HOXB1, ASCL1, PRLH, PHOX2B, ISX
Warning in FetchData.Seurat(object = acc_cancer_pri, vars = genes, slot = "data") :
  The following requested variables were not found: HOXB1, PHOX2B
Warning in FetchData.Seurat(object = acc_cancer_pri, vars = genes, slot = "data") :
  The following requested variables were not found: DRD3, MIR1-1, MIR133A1

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)

10 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