A human lysosomal storage disorder toolkit for decoding proteome landscapes in cortical-like and dopaminergic-like induced neurons.

Code ↔ Paper

7 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 7 matches

[1] § Results › A Stem Cell Toolkit for Systematic Analysis of LSD Genes. ↔ proteome/diff132_d50_nDIA.Rmd, lines 243–306 · score 0.81 · ATP13A2, CLN1, CLN2, CLN7, CTSF, DNAJC5
[2] § Results › LSD Mutant iNeuron Organelle Landscapes. ↔ proteome/diff132_d50_nDIA.Rmd, lines 510–602 · score 0.76 · technical noise, biological signals, genotype variance, organelle annotations, log2fc, ratio
[3] § Results › Proteomic Fingerprint Resource for LSD Mutant iN and iDA Neurons. ↔ proteome/diff132_d50_nDIA.Rmd, lines 2310–2349 · score 0.66 · Integral membrane protein, SynapseSV, recycling endosome, Disease class, Ternary, log2FC
[4] § Results › LSD Mutant iNeuron Organelle Landscapes. ↔ proteome/diff132_d50_nDIA.Rmd, lines 310–349 · score 0.62 · neuronal ceroid lipofuscinoses, integral membrane protein, disease class, IDs, disorders, GRN
[5] § Results › LSD Mutant iNeuron Organelle Landscapes. ↔ proteome/diff132_d50_nDIA.Rmd, lines 12–114 · score 0.61 · subcellular compartments, synaptic compartments, impact score, disease classes, knockout, variance
[6] § Results › Visualization of ASAH1−/− Endolysosomal Compartments by Cryo-ET. ↔ lipidome/Lipidomics_HeLa_iN-diff133_ASAH1-WC-OrganellIeIP.Rmd, lines 85–198 · score 0.53 · LysoIP, HeLa, GM1, lipidomics, GM2, GM3
[7] § Results › Proteomic Fingerprint Resource for LSD Mutant iN and iDA Neurons. ↔ proteome/diff132_d50_nDIA.Rmd, lines 243–306 · score 0.51 · hrMS2, nDIA, Astral, GRN, triplicate, SMPD1

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 4,896 lines · 195 KB · MIT · 6 matches

---
title: "diff132_d50_nDIA"
output: html_document
date: "`r Sys.Date()`"
chunk_output_type: console
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# =================================================
# Module 0: Setup & Configuration & Experimental Background
# =================================================
# Module overview & aims
Module 1: Imports replicate-level proteomics data, performs quality control through replicate correlation analysis, F-ratio variance testing (biological vs. technical noise), and PCA visualization across genotypes, neuron types, and differentiation timepoints.
Module 2: Processes summarized fold-change data with subcellular annotations, generates global heatmaps, validates LSD knockout efficiency, and visualizes organelle-level proteome changes using ternary plots, half-circos plots, and annotation-specific PCA for synaptic compartments.
Module 3: Validates neuronal identity through pre/postsynaptic and dopaminergic marker expression, then profiles synaptic machinery (v-ATPase, SNARE, SV docking complexes) across genotypes with ranked log2FC plots and ChimeraX-compatible color tables for structure mapping.
Module 4: Compares day 30 versus day 50 proteomes through neuronal maturation marker trajectories, pre/postsynaptic protein dynamics, and annotation-level log2FC violin plots with Wilcoxon statistics for ASAH1 and GBA1 knockouts across differentiation timepoints.
Module 5: Integrates day 70 ASAH1 knockout data with day 30/50 to build a three-timepoint longitudinal analysis, including PCA trajectories, annotation-level violin plots across differentiation, pairwise annotation correlations, and ChimeraX color-mapping commands for mitochondrial Complex I structure visualization.
Module 6: Provides reusable functions for generating annotation-stratified violin plots of log2FC across genotypes and volcano plots with customizable annotation highlights, plus iN-vs-iDA correlation scatter plots preserving top hits from both neuron types.
Module 7: Computes genotype-by-genotype and annotation-by-annotation Spearman correlations of log2FC values, generates impact scoring heatmaps ranking genotypes by their deviation from the group, and produces bubble plots comparing iN vs iDA correlation patterns across disease classes and subcellular compartments.
Module 8: Performs disease-class-focused analyses including sphingolipidosis subset heatmaps and correlations, cross-experiment ASAH1 timecourse comparisons (diff118 vs diff132), Euler diagrams for genotype overlap within annotations, limma-based differential analysis between genotypes, and circos plots visualizing annotation-level correlations stratified by disease class.
Module 9:Linear Regression
Module 10: Performs a comprehensive cross-omics correlation analysis between neuron (iN/iDA) nDIA proteomes and HeLa nMOST proteomes across LSD genotypes, integrating GO and organelle annotations to visualize genotype-, cluster-, and annotation-specific correlation structure using boxplots, heatmaps, and multi-layer circos plots.
Module 11: Builds an integrated PPI vulnerability and rewiring pipeline for neuronal proteomes (iN and iDA) across lysosomal gene knockouts by merging experimental proteomics with predicted and validated interaction networks. It quantifies network destabilization, identifies organelle- and neuron-type–specific vulnerabilities, and visualizes PPI loss, rewiring, and functional impacts using network, enrichment, and comparative analyses.
Module 12: Extends the PPI vulnerability framework to HeLa whole-cell proteome data from nMOST, computes control-normalized baseline stability, identifies lost/retained edges and broken trimers/complexes per KO, then performs cross-cell-type comparisons (iN, iDA, HeLa) using streamgraphs and matching pattern analysis to identify cell-type-specific versus shared PPI vulnerabilities.
Module 13: Validates predicted PPI networks against XL-MS cross-linking data (DSSO-CLASP from Zhu et al.) by mapping gene names to UniProt IDs, computing edge-key overlaps between CLASP-detected interactions and neuronal/HeLa baseline networks, and quantifying both forward (CLASP in baseline) and reverse (mito baseline in CLASP) validation rates with mitochondrial annotation breakdowns.
# Load packackes & set project directory
```{r setup, include=FALSE}
# install.packages(c(
# "plyr", "ggrepel", "ggpubr", "ggpmisc", "patchwork", "gghighlight",
# "ggVennDiagram", "UpSetR", "ComplexUpset", "pheatmap", "circlize",
# "magick", "viridis", "NatParksPalettes", "limma", "bigstatsr",
# "irlba", "ggbiplot", "plotly", "Cairo", "devtools", "lintr"
# ))
# Core Data Manipulation & Tidyverse
library(tidyverse) # Includes ggplot2, dplyr, tidyr, readr, tibble, purrr, stringr, forcats
library(plyr)
library(dplyr)
library(data.table)
library(broom)
library(rlang)
library(tidyr)
# Visualization: ggplot2 Extensions & Plots
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(ggpmisc)
library(cowplot)
library(patchwork)
library(scales)
library(gghighlight)
library(superheat)
library(ggVennDiagram)
library(UpSetR)
library(ComplexUpset)
install.packages("BiocManager")
BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)
#install.packages("eulerr")
library(eulerr)
#install.packages("Ternary")
library(Ternary)
#devtools::install_github("davidsjoberg/ggsankey")
library(ggsankey)
# Heatmaps & Clustering
library(pheatmap)
library(RColorBrewer)
library(circlize)
library(magick)
# Color Palettes
library(viridis)
library(NatParksPalettes)
# Dimension Reduction / Statistics
library(limma)
library(bigstatsr)
library(irlba)
library(ggbiplot)
# Interactive & Advanced Plotting
library(plotly)
library(circlize)
# Utilities / Dev / Misc
library(grid)
library(gridExtra)
#install.packages("Cairo")
library(Cairo)
#library(png)
library(devtools)
library(lintr)
```
# Output dirs overview
```{r}
out_dir_d50_QC <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/QC"
dir.create(out_dir_d50_QC, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_timecourse <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/timecourse"
dir.create(out_dir_d50_timecourse, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_euler <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/QC/euler_annotations"
dir.create(out_dir_d50_euler, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_ascd_avg <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/ascendingAvg"
dir.create(out_dir_d50_ascd_avg, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_ascd_avg_SV <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/ascendingAvg/synapsevar"
dir.create(out_dir_d50_ascd_avg_SV, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_heatmaps <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/heatmaps"
dir.create(out_dir_d50_heatmaps, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_dataframes <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/dataframes"
dir.create(out_dir_d50_dataframes, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_violins <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/Organelle_Violin"
dir.create(out_dir_d50_violins, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_volcano <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/Volcano"
dir.create(out_dir_d50_volcano, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_SphMut <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/SphMut"
dir.create(out_dir_d50_SphMut, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_SphMuteuler <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/SphMut/euler"
dir.create(out_dir_d50_SphMuteuler, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_DisClass <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/DisClass"
dir.create(out_dir_d50_DisClass, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_DisClass_circos <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/DisClass/Cirosdiagram"
dir.create(out_dir_d50_DisClass_circos, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_DisClass_splitcorr_dir <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/DisClass/Cirosdiagram/splitcorr"
dir.create(out_dir_d50_DisClass_splitcorr_dir, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_GRN <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/GRN"
dir.create(out_dir_d50_GRN, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_CrossCorr <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/Organelle_CrossCorr"
dir.create(out_dir_d50_CrossCorr, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_CrossCorr_annotavg <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/Organelle_CrossCorr/annot_avg"
dir.create(out_dir_d50_CrossCorr_annotavg, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_CrossCorr_nMOSTGOavgannot <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/Organelle_CrossCorr/nMOSTGOavgannnot"
dir.create(out_dir_d50_CrossCorr_nMOSTGOavgannot, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_CrossCorr_genoavg <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/Organelle_CrossCorr/geno_avg"
dir.create(out_dir_d50_CrossCorr_genoavg, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_CrossCorr_bubble <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/Organelle_CrossCorr/bubble"
dir.create(out_dir_d50_CrossCorr_bubble, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_CrossCorr_bubble_ASAH1 <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/Organelle_CrossCorr/ASAH1_c23d30d50"
dir.create(out_dir_d50_CrossCorr_bubble_ASAH1, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_corrProtein <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/proteincorr"
dir.create(out_dir_d50_corrProtein, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_SynPRM <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/SynPRM"
dir.create(out_dir_d50_SynPRM, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_LinReg <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/LinReg"
dir.create(out_dir_d50_LinReg, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_LinReg_limma <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/LinReg/limma_all_coef"
dir.create(out_dir_d50_LinReg_limma, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_LinReg_limma_select <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/LinReg/limma_select"
dir.create(out_dir_d50_LinReg_limma_select, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_nMOST_nDIAcorrelation <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/nMOST-nDIAcorrelation"
dir.create(out_dir_d50_nMOST_nDIAcorrelation, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_nMOST_nDIAcorrelationS_single <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/nMOST-nDIAcorrelation/SingleGeno"
dir.create(out_dir_d50_nMOST_nDIAcorrelationS_single, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_nMOST_nDIAcorrelationS_OrganelleComp <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/nMOST-nDIAcorrelation/HeLa_Neuron_SelectAnnotation_comparisons"
dir.create(out_dir_d50_nMOST_nDIAcorrelationS_OrganelleComp, recursive = TRUE, showWarnings = FALSE)
out_dir_d50_nMOST_nDIAcorrelationS_Circos <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/nMOST-nDIAcorrelation/HeLa_Neuron_SelectAnnotation_comparisons"
dir.create(out_dir_d50_nMOST_nDIAcorrelationS_Circos, recursive = TRUE, showWarnings = FALSE)
out_dir_d70_dataframes <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day70/dataframes"
dir.create(out_dir_d70_dataframes, recursive = TRUE, showWarnings = FALSE)
out_dir_d70_heatmap <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day70/heatmap"
dir.create(out_dir_d70_heatmap, recursive = TRUE, showWarnings = FALSE)
out_dir_d70_timecourse <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day70/timecourse"
dir.create(out_dir_d70_timecourse, recursive = TRUE, showWarnings = FALSE)
out_dir_d70_pca <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day70/pca"
dir.create(out_dir_d70_pca, recursive = TRUE, showWarnings = FALSE)
out_dir_d70_violins <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day70/violins"
dir.create(out_dir_d70_violins, recursive = TRUE, showWarnings = FALSE)
out_dir_d70_corr <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day70/corr"
dir.create(out_dir_d70_corr, recursive = TRUE, showWarnings = FALSE)
out_dir_PPI <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/PPI"
dir.create(out_dir_PPI, recursive = TRUE, showWarnings = FALSE)
```
# Experimental Background: nDIA proteomics iN & iDA d30+d50
- Samples are in triplicates whole-cell samples from 12-well dishes.
- day 30 of in vitro differentiation of iN and iDA:
Genotypes: Ctrl, ASAH1, CLN11.GRN, GBA1, SMPD1
- day 50 of in vitro differentiation of iN and iDA:
Genotypes: Ctrl, ASAH1, CLN1.PPT1, CLN3, CLN4.DNAJC5, CLN5, CLN6, CLN7.MFSD8, CLN8, CLN10.CTSD, CLN11.GRN, CLN12.ATP13A2, CLN13.CTSF, CTNS, GAA, GBA1, GLB1, HEXA, HEXB, LIPA, MCOLN1.TRPML1, NPC1, NPC2, PSAP, SMPD1
run on Thermo Orbitrap Astral, nDIA, hrMS2
# Circos plot for neuroLSD
```{r}
library(circlize)
#read the csv file
#read in first column (==gene) as row title
LSD <- read.csv('/Users/felix/Documents/PostDoc/Harvard/03_LSD/Proteomics_Lipidomics/Data_forplots/Circos/LSD_circos.csv',
header=T, row.names=1)
#convert the table to a martix
data <- as.matrix(LSD)
#create a Ciros diagram
chordDiagram(data)
############## CIRCOS Plot with confirmed KOs #####################
## plot only confirmed KOs in color, rest of genes in grey/ per disease group
## save plot as pdf
# set colors for disease groups
col = c(Sphingolipidoses="#D53E4F", Mucopolysaccharidoses="#F46D43", Glycoproteinoses="#FDAE61", Neuronal.ceroid.lipofuscinoses="#FEE08B",
Integral.membrane.protein.disorders="#E6F598", PTM.defects="#ABDDA4", Lipid.storage.diseases= "#66C2A5",
AGA="grey76", ARSA="grey76", ARSB="grey76", ASAH1="#D53E4F", CLN1.PPT1="#FEE08B", CLN10.CTSD="#FEE08B", CLN11.GRN="#FEE08B", CLN12.ATP13A2="#FEE08B",CLN12.ATP13A2="#E6F598", CLN13.CTSF="#FEE08B", CLN14.KCTD7="grey76", CLN2.TPP1="#FEE08B", CLN3="#FEE08B", CLN4.DNAJC5="#FEE08B", CLN5="#FEE08B", CLN6="#FEE08B", CLN7.MFSD8="#FEE08B", CLN8="#FEE08B", CTNS="grey76", CTSA="grey76", FUCA="grey76", GAA="#66C2A5", GALC="grey76", GALNS="grey76", GBA1="#D53E4F", GLA="grey76", GLB1="grey76", GM2A="grey76", GNPTAB="grey76", GNPTG="grey76", GNS="grey76", GUSB="grey76", HEXA="#D53E4F", HEXB="#D53E4F", HGSNAT="grey76", HYAL1="grey76", IDS="grey76", IDUA="grey76", LAMP2="grey76", LIPA="#66C2A5", MAN2B1="grey76", MANBA="grey76", MCOLN1="#E6F598", NAGA="grey76", NAGLU="grey76", NEU1="grey76", NPC1="#E6F598", NPC2="#E6F598", PSAP="#D53E4F", SCARB2="grey76", SGSH="grey76", SLC17A5="grey76", SMPD1="#D53E4F", SUMF1="grey76")
# Set output path and size
pdf("/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/LSD_circos_plot_confirmed_KOs.pdf", width = 10, height = 10) # adjust size as needed
# Plot: chord diagram
chordDiagram(data,
grid.col = col,
annotationTrack = "grid",
preAllocateTracks = 1,
transparency = 0.1,
link.lwd = 1,
link.lty = 1,
link.border = 1)
# Add labels
circos.trackPlotRegion(track.index = 2, panel.fun = function(x, y) {
xlim = get.cell.meta.data("xlim")
ylim = get.cell.meta.data("ylim")
sector.name = get.cell.meta.data("sector.index")
circos.text(mean(xlim), ylim[1] + 2.5, sector.name,
facing = "clockwise", niceFacing = TRUE,
adj = c(0, 0.5), cex = 0.6)
}, bg.border = "black")
# Finish writing the PDF
dev.off()
```
# =================================================
# Module 1: Data Import, QC (Replicate-level)
# =================================================
# top level disese-class annotation
```{r}
disease_classes <- list(
Sphingolipidoses = c("ARSA", "ASAH1", "GALC", "GBA1", "GLA", "GLB1", "GM2A", "HEXA", "HEXB", "PSAP", "SMPD1"),
Mucopolysaccharidoses = c("GLB1", "ARSB", "GALNS", "GNS", "GUSB", "HGSNAT", "HYAL1", "IDS", "IDUA", "NAGLU", "SGSH"),
Glycoproteinoses = c("AGA", "CTSA", "FUCA", "MAN2B1", "MANBA", "NAGA", "NEU1"),
Neuronal.ceroid.lipofuscinoses = c("PPT1", "CTSD", "GRN", "ATP13A2", "CTSF", "KCTD7", "TPP1", "CLN3", "DNAJC5", "CLN5", "CLN6", "MFSD8", "CLN8"),
Integral.membrane.protein.disorders = c("ATP13A2", "CLN3", "CTNS", "LAMP2", "MCOLN1", "NPC1", "NPC2", "SCARB2", "SLC17A5"),
PTM.defects = c("GNPTAB", "GNPTG", "SUMF1"),
Lipid.storage.diseases = c("GAA", "LIPA")
)
# convert into df
disease_classes_dataframe <- stack(disease_classes)
colnames(disease_classes_dataframe) <- c("sample", "DiseaseClass")
# give it a color
disease_class_palette <- c(
Sphingolipidoses = "#D53E4F",
Mucopolysaccharidoses = "#F46D43",
Glycoproteinoses = "#FDAE61",
Neuronal.ceroid.lipofuscinoses = "#FEE08B",
Integral.membrane.protein.disorders = "#E6F598",
PTM.defects = "#ABDDA4",
Lipid.storage.diseases = "#66C2A5"
)
# create disease_class_df
disease_class_df <- data.frame(
sample = c("ARSA", "ASAH1", "GALC", "GBA1", "GLA", "GLB1", "GM2A", "HEXA", "HEXB", "PSAP", "SMPD1", "GLB1", "ARSB", "GALNS", "GNS", "GUSB", "HGSNAT", "HYAL1", "IDS", "IDUA", "NAGLU", "SGSH", "AGA", "CTSA", "FUCA", "MAN2B1", "MANBA", "NAGA", "NEU1", "PPT1", "CTSD", "GRN", "ATP13A2", "CTSF", "KCTD7", "TPP1", "CLN3", "DNAJC5", "CLN5", "CLN6", "MFSD8", "CLN8", "ATP13A2", "CLN3", "CTNS", "LAMP2", "MCOLN1", "NPC1", "NPC2", "SCARB2", "SLC17A5", "GNPTAB", "GNPTG", "SUMF1", "GAA", "LIPA"),
DiseaseClass = c(rep("Sphingolipidoses", 11), rep("Mucopolysaccharidoses", 11), rep("Glycoproteinoses", 7), rep("Neuronal.ceroid.lipofuscinoses", 13), rep("Integral.membrane.protein.disorders", 9), rep("PTM.defects", 3), rep("Lipid.storage.diseases", 2))
)
```
# Data cleaning of replicate data
plot PCA plots to check data landscape on top level
Have following genotypes and cells in there
iN, DA, pool (== QC mmaker)
3 replicates for genotypes
```{r}
# --------------------------------------------------- #
# Data cleaning and preparation
# --------------------------------------------------- #
# --- Step 1: Read in the data from the CSV file
# The CSV file has been pre-processed in Excel to replace all NA values with 0.
diff132_d50_PCA_df_full <- read.csv('/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/datasets/lsd_pro_long.csv', header = TRUE, sep = ',', stringsAsFactors = FALSE)
# --- Step 1.5: Fill NA or empty values in 'neuron' and 'replicate' with QC info
diff132_d50_PCA_df_full <- diff132_d50_PCA_df_full %>%
mutate(
is_qc = is.na(neuron) | neuron == "",
neuron = ifelse(is_qc, "QC", neuron),
replicate = ifelse(is_qc, as.character(seq_len(sum(is_qc))), replicate)
)
# --- Step 2: Inspect the structure of the data
# This helps to understand the data types and identify which columns to retain or remove.
str(diff132_d50_PCA_df_full)
# --- Step 3: Remove unnecessary columns
# Assuming columns 1, 2, 4, 5, and 10 are not needed for PCA, they are removed.
diff132_d50_PCA_df <- diff132_d50_PCA_df_full[, -c(1,2,4,6,7,13,14)]
# --- Step 4: Identify and remove rows with missing or empty 'Genes' entries
# It's crucial to ensure that each gene has a valid identifier for accurate analysis.
# Identify problematic rows where 'Genes' is NA or an empty string
problematic_rows <- diff132_d50_PCA_df %>%
filter(is.na(Genes) | Genes == "")
# Remove these problematic rows from the dataset
cleaned_df <- diff132_d50_PCA_df %>%
filter(!is.na(Genes) & Genes != "") %>% # remove empty gene names
filter(!is.na(quan)) # remove missing values in quan
# --- Step 4a: Load subcellular annotations
cols_to_read <- c(1,2,3,4,5,6,7,9,10,12,15,16,17,18,19,20,21,22,23,25,26,27,28,29,30,31,32,33,34)
subcell_df <- read.csv("/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/tmpneurogit/proteome/SubCellAnnotation.csv", stringsAsFactors = FALSE)[, cols_to_read]
subcell_df[] <- lapply(subcell_df, as.character)
long_subcell_df <- subcell_df %>%
pivot_longer(cols = everything(), names_to = "Localization", values_to = "Genes") %>%
filter(!is.na(Genes) & Genes != "") %>%
mutate(Genes = trimws(Genes)) %>%
distinct()
```
# Replicate correlation analysis
```{r}
# --- Step 4b: Replicate correlation summary
rep_corr <- cleaned_df %>%
dplyr::filter(neuron != "QC", !is.na(genotype)) %>%
dplyr::group_by(Genes, day, neuron, genotype, replicate) %>%
dplyr::summarise(quan = sum(quan), .groups = "drop") %>%
tidyr::pivot_wider(names_from = replicate, values_from = quan, names_prefix = "rep") %>%
dplyr::group_by(day, neuron, genotype) %>%
dplyr::summarise(
r12 = cor(rep1, rep2, use = "complete.obs"),
r13 = cor(rep1, rep3, use = "complete.obs"),
r23 = cor(rep2, rep3, use = "complete.obs"),
.groups = "drop"
) %>%
tidyr::pivot_longer(cols = c(r12, r13, r23), names_to = "pair", values_to = "correlation")
summary(rep_corr$correlation)
# Distribution plot
ReplicateDist <- ggplot(rep_corr, aes(x = correlation)) +
geom_histogram(bins = 30, fill = "dodgerblue4", color = "white") +
geom_vline(xintercept = 0.95, linetype = "dashed", color = "red") +
labs(x = "Pairwise replicate correlation (Pearson)", y = "Count") +
theme_bw(base_size = 6) +
theme(panel.grid = element_blank())
ReplicateDist2 <- ggplot(rep_corr, aes(x = correlation, fill = neuron, color = neuron)) +
geom_density(alpha = 0.5) +
scale_fill_manual(values = c("iN" = "grey50", "iDA" = "grey80")) +
scale_color_manual(values = c("iN" = "grey50", "iDA" = "grey80")) +
geom_vline(xintercept = 0.95, linetype = "dashed", color = "firebrick") +
labs(x = "Pairwise replicate correlation (Pearson)", y = "Density", fill = NULL, color = NULL) +
theme_bw(base_size = 6) +
theme(panel.grid = element_blank(), legend.position = "bottom")
ggsave(file.path(out_dir_d50_QC, "diff132_replicateDataDistribution.pdf"), ReplicateDist, width = 4, height = 4, device = cairo_pdf)
ggsave(file.path(out_dir_d50_QC, "diff132_replicateDataDistributionDensity.pdf"), ReplicateDist2, width = 4, height = 4, device = cairo_pdf)
# --- correlation matrix of samples
rep_wide_df <- cleaned_df %>%
dplyr::filter(neuron != "QC", !is.na(genotype)) %>%
dplyr::group_by(Genes, day, neuron, genotype, replicate) %>%
dplyr::summarise(quan = sum(quan), .groups = "drop") %>%
dplyr::mutate(sample_id = paste(genotype, neuron, day, replicate, sep = "_")) %>%
dplyr::arrange(day, neuron, genotype, replicate) %>%
dplyr::select(Genes, sample_id, quan) %>%
tidyr::pivot_wider(names_from = sample_id, values_from = quan)
sample_order <- cleaned_df %>%
dplyr::filter(neuron != "QC", !is.na(genotype)) %>%
dplyr::distinct(day, neuron, genotype, replicate) %>%
dplyr::arrange(day, neuron, genotype, replicate) %>%
dplyr::mutate(sample_id = paste(genotype, neuron, day, replicate, sep = "_")) %>%
dplyr::pull(sample_id)
rep_corr_matrix <- cor(rep_wide_df[, sample_order], use = "pairwise.complete.obs")
pdf(file.path(out_dir_d50_QC, "diff132_replicate_correlation_heatmap.pdf"), width = 12, height = 12)
pheatmap::pheatmap(
rep_corr_matrix,
color = colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdYlBu")))(100),
breaks = seq(0.5, 1, length.out = 101),
cluster_rows = FALSE,
cluster_cols = FALSE,
cellheight = 4,
cellwidth = 4,
border_color = NA,
fontsize = 4,
main = "Sample-wise Pearson Correlation"
)
dev.off()
# --- Within vs between replicate variance plots
rep_variance_per_gene <- cleaned_df %>%
dplyr::filter(neuron != "QC", genotype != "ctrl", !is.na(genotype)) %>%
dplyr::group_by(Genes, day, neuron, genotype) %>%
dplyr::filter(dplyr::n() == 3) %>%
dplyr::summarise(var_within = var(quan, na.rm = TRUE), .groups = "drop") %>%
dplyr::mutate(sample_id = paste(genotype, neuron, day, sep = "_"))
sample_ids <- unique(rep_variance_per_gene$sample_id)
sample_colors <- setNames(colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdYlBu")))(length(sample_ids)), sample_ids)
rep_variance_scatter <- ggplot(rep_variance_per_gene, aes(x = Genes, y = log2(var_within), color = sample_id)) +
geom_point(size = 0.1, alpha = 0.3) +
scale_color_manual(values = sample_colors) +
labs(x = "Proteins", y = "Within-replicate variance [log2]", color = NULL) +
theme_bw(base_size = 6) +
theme(
panel.grid = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "right"
)
ggsave(file.path(out_dir_d50_QC, "diff132_within_replicate_variance_scatter.pdf"), rep_variance_scatter, width = 8, height = 4, device = cairo_pdf)
```
# F-ratio variance analysis
```{r}
# --- Step 4c: F-ratio approach per annotation
# F-ratio = between-genotype variance / within-genotype variance
# MSbetween: variance between genotypes (based on all replicates)
# MSwithin: variance within genotypes (across replicates)
# If F > 1: biological signal is larger than technical noise = real differences
# If F < 1: technical noise is larger than biological signal = can't trust the differences
# also have log2 implementation and then it is F > 0 as cuttoff since log2(1) = 0
# --- Create replicate-level fold changes per organelle annotation
rep_fc_organelle <- cleaned_df %>%
dplyr::filter(neuron != "QC", !is.na(genotype)) %>%
dplyr::group_by(Genes, day, neuron, genotype, replicate) %>%
dplyr::summarise(quan = sum(quan), .groups = "drop") %>%
# pivot to get ctrl and KO quantities side by side
tidyr::pivot_wider(names_from = genotype, values_from = quan) %>%
# reshape back to long, keeping ctrl as reference
tidyr::pivot_longer(cols = -c(Genes, day, neuron, replicate, ctrl), names_to = "genotype", values_to = "ko_quan") %>%
dplyr::rename(ctrl_quan = ctrl) %>%
dplyr::filter(!is.na(ko_quan), !is.na(ctrl_quan)) %>%
# calculate log2 fold change vs ctrl
dplyr::mutate(log2fc = log2(ko_quan / ctrl_quan)) %>%
# add subcellular localization annotation
dplyr::left_join(long_subcell_df, by = "Genes") %>%
dplyr::filter(!is.na(Localization))
# --- Calculate F-ratio per gene per annotation
f_ratio_per_gene <- rep_fc_organelle %>%
dplyr::group_by(Localization, day, neuron, genotype, Genes) %>%
dplyr::filter(dplyr::n() == 3) %>%
dplyr::ungroup() %>%
dplyr::group_by(Localization, day, neuron, Genes) %>%
dplyr::summarise(
ms_within = mean(tapply(log2fc, genotype, var, na.rm = TRUE), na.rm = TRUE),
ms_between = var(log2fc, na.rm = TRUE),
f_ratio = ms_between / ms_within,
log2_f_ratio = log2(ms_between / ms_within),
.groups = "drop"
) %>%
dplyr::filter(is.finite(f_ratio), is.finite(log2_f_ratio))
# --- Summarize F-ratio per annotation x day x neuron
f_ratio_df <- f_ratio_per_gene %>%
dplyr::group_by(Localization, day, neuron) %>%
dplyr::summarise(
mean_f = mean(f_ratio, na.rm = TRUE),
sd_f = sd(f_ratio, na.rm = TRUE),
se_f = sd_f / sqrt(dplyr::n()),
mean_log2f = mean(log2_f_ratio, na.rm = TRUE),
sd_log2f = sd(log2_f_ratio, na.rm = TRUE),
se_log2f = sd_log2f / sqrt(dplyr::n()),
.groups = "drop"
)
write.csv(f_ratio_per_gene, file.path(out_dir_d50_QC, "diff132_f_ratio_per_gene_sourcedata.csv"), row.names = FALSE)
write.csv(f_ratio_df, file.path(out_dir_d50_QC, "diff132_f_ratio_summary_sourcedata.csv"), row.names = FALSE)
# --- Plot F-ratio (linear)
f_ratio_plot <- ggplot(f_ratio_df, aes(x = Localization, y = mean_f, color = neuron, fill = neuron, group = neuron)) +
geom_ribbon(aes(ymin = mean_f - se_f, ymax = mean_f + se_f), alpha = 0.2, color = NA) +
geom_line() +
geom_point(size = 0) +
geom_hline(yintercept = 1, linetype = "dashed", color = "firebrick") +
scale_color_manual(values = c("iN" = "grey50", "iDA" = "grey80")) +
scale_fill_manual(values = c("iN" = "grey50", "iDA" = "grey80")) +
facet_wrap(~ day, ncol = 2) +
labs(x = NULL, y = "F-ratio (between / within genotype variance)", color = NULL, fill = NULL,
caption = "F > 1: biological signal > technical noise\nF < 1: technical noise > biological signal") +
theme_bw(base_size = 6) +
theme(panel.grid = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
legend.position = "bottom", plot.caption = element_text(hjust = 0))
# --- Plot F-ratio (log2)
f_ratio_log2_plot <- ggplot(f_ratio_df, aes(x = Localization, y = mean_log2f, color = neuron, fill = neuron, group = neuron)) +
geom_ribbon(aes(ymin = mean_log2f - se_log2f, ymax = mean_log2f + se_log2f), alpha = 0.2, color = NA) +
geom_line() +
geom_point(size = 0) +
geom_hline(yintercept = 0, linetype = "dashed", color = "firebrick") +
scale_color_manual(values = c("iN" = "grey50", "iDA" = "grey80")) +
scale_fill_manual(values = c("iN" = "grey50", "iDA" = "grey80")) +
facet_wrap(~ day, ncol = 2) +
labs(x = NULL, y = "log2 F-ratio (between / within genotype variance)", color = NULL, fill = NULL,
caption = "F > 0: biological signal > technical noise\nF < 0: technical noise > biological signal") +
theme_bw(base_size = 6) +
theme(panel.grid = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
legend.position = "bottom", plot.caption = element_text(hjust = 0))
ggsave(file.path(out_dir_d50_QC, "diff132_f_ratio_variance.pdf"), f_ratio_plot, width = 4, height = 2, device = cairo_pdf)
ggsave(file.path(out_dir_d50_QC, "diff132_f_ratio_variance_log2.pdf"), f_ratio_log2_plot, width = 4, height = 2, device = cairo_pdf)
```
# PCA new
```{r}
# =============================================================================
# PCA plots
# =============================================================================
# --- Load data
lsd_pro_phen <- readr::read_csv("/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/tmpneurogit/proteome/lsd_pro_phen.csv", show_col_types = FALSE)
.quan_raw <- readr::read_csv("/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/tmpneurogit/proteome/lsd_pro_quan_na.csv", show_col_types = FALSE)
lsd_pro_quan_na <- as.matrix(.quan_raw[, -1]); rownames(lsd_pro_quan_na) <- .quan_raw[[1]]; rm(.quan_raw)
# --- Align phenotype to matrix column order
phen_ord <- lsd_pro_phen[match(colnames(lsd_pro_quan_na), lsd_pro_phen$sample), ]
is_qc <- is.na(phen_ord$genotype)
# --- PCA: all samples
mat_all <- lsd_pro_quan_na[complete.cases(lsd_pro_quan_na), ]
pca_all <- prcomp(t(mat_all), scale. = TRUE)
var_all <- round(100 * pca_all$sdev^2 / sum(pca_all$sdev^2), 1)
pc_df <- data.frame(PC1 = pca_all$x[, 1], PC2 = pca_all$x[, 2],
genotype = ifelse(is_qc, "", phen_ord$genotype),
neuron = ifelse(is_qc, "QC", phen_ord$neuron),
day = factor(phen_ord$day))
# --- PCA: noQC
mat_noqc <- lsd_pro_quan_na[, !is_qc]; mat_noqc <- mat_noqc[complete.cases(mat_noqc), ]
pca_noqc <- prcomp(t(mat_noqc), scale. = TRUE)
var_noqc <- round(100 * pca_noqc$sdev^2 / sum(pca_noqc$sdev^2), 1)
phen_noqc <- phen_ord[!is_qc, ]
pc_df_noqc <- data.frame(PC1 = pca_noqc$x[, 1], PC2 = pca_noqc$x[, 2],
genotype = phen_noqc$genotype, neuron = phen_noqc$neuron, day = factor(phen_noqc$day))
# --- Palettes
n_geno_nqc <- length(unique(pc_df_noqc$genotype))
n_days <- length(unique(pc_df$day))
n_days_nqc <- length(unique(pc_df_noqc$day))
geno_cols <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdYlBu"))(length(unique(pc_df$genotype)))
geno_cols_nqc <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdYlBu"))(n_geno_nqc)
day_cols <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdYlBu"))(n_days)
day_cols_nqc <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdYlBu"))(n_days_nqc)
neuron_cols <- c("iDA" = "#A50026", "iN" = "#313695", "QC" = "grey50")
lab_x <- paste0("PC1 (", var_all[1], "% variance)"); lab_y <- paste0("PC2 (", var_all[2], "% variance)")
lab_xn <- paste0("PC1 (", var_noqc[1], "% variance)"); lab_yn <- paste0("PC2 (", var_noqc[2], "% variance)")
bt <- theme_minimal() + theme(panel.border = element_rect(color = "black", fill = NA, linewidth = 1))
# --- All-sample plots
d50plot1 <- ggplot(pc_df, aes(x = PC1, y = PC2, color = genotype)) +
geom_point(size = 3, alpha = 0.8) +
geom_text_repel(data = dplyr::distinct(pc_df, genotype, .keep_all = TRUE), aes(label = genotype), size = 2, box.padding = 0.5, point.padding = 0.5, max.overlaps = 200) +
scale_color_manual(values = geno_cols) +
labs(title = "PCA of Gene Intensities Colored by Genotype", x = lab_x, y = lab_y) +
coord_fixed(ratio = 1.75) + bt
d50plot2 <- ggplot(pc_df, aes(x = PC1, y = PC2, color = neuron)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(aes(fill = neuron), geom = "polygon", alpha = 0.2, level = 0.95) +
scale_color_manual(values = neuron_cols) + scale_fill_manual(values = neuron_cols) +
labs(title = "PCA of Gene Intensities Colored by Neuron Type", x = lab_x, y = lab_y) +
coord_fixed(ratio = 1) + bt
d50plot3 <- ggplot(pc_df, aes(x = PC1, y = PC2, color = day)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(aes(fill = day), geom = "polygon", alpha = 0.2, level = 0.95) +
scale_color_manual(values = day_cols) + scale_fill_manual(values = day_cols) +
labs(title = "PCA of Gene Intensities Colored by Day of Differentiation", x = lab_x, y = lab_y) +
coord_fixed(ratio = 1) + bt
# --- noQC plots
d50plot1_noQC <- ggplot(pc_df_noqc, aes(x = PC1, y = PC2, color = genotype, alpha = day)) +
geom_point(size = 3) +
scale_alpha_manual(values = setNames(seq(0.4, 1, length.out = n_days_nqc), sort(unique(pc_df_noqc$day)))) +
geom_text_repel(data = dplyr::distinct(pc_df_noqc, genotype, .keep_all = TRUE), aes(label = genotype), color = "black", size = 6, box.padding = 0.5, point.padding = 0.5, max.overlaps = 1000) +
scale_color_manual(values = geno_cols_nqc) +
labs(title = "PCA of Gene Intensities Colored by Genotype", x = lab_xn, y = lab_yn) +
coord_fixed(ratio = 0.5) + bt + theme(legend.position = "bottom")
d50plot2_noQC <- ggplot(pc_df_noqc, aes(x = PC1, y = PC2, color = neuron)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(aes(fill = neuron), geom = "polygon", alpha = 0.2, level = 0.998) +
scale_color_manual(values = neuron_cols) + scale_fill_manual(values = neuron_cols) +
labs(title = "PCA of Gene Intensities Colored by Neuron Type", x = lab_xn, y = lab_yn) +
coord_fixed(ratio = 0.5) + bt
d50plot3_noQC <- ggplot(pc_df_noqc, aes(x = PC1, y = PC2, color = day)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(aes(fill = day), geom = "polygon", alpha = 0.2, level = 0.95) +
scale_color_manual(values = day_cols_nqc) + scale_fill_manual(values = day_cols_nqc) +
labs(title = "PCA of Gene Intensities Colored by Day of Differentiation", x = lab_xn, y = lab_yn) +
coord_fixed(ratio = 1) + bt
# --- Display
d50plot1; d50plot2; d50plot3
d50plot1_noQC; d50plot2_noQC; d50plot3_noQC
# --- Save
ggsave(file.path(out_dir_d50_QC, "diff132_d50_PCA_byGenotype.pdf"), d50plot1, width = 15, height = 15, units = "cm", dpi = 600)
ggsave(file.path(out_dir_d50_QC, "diff132_d50_PCA_byCelltype.pdf"), d50plot2, width = 15, height = 15, units = "cm", dpi = 600)
ggsave(file.path(out_dir_d50_QC, "diff132_d50_PCA_byDay.pdf"), d50plot3, width = 15, height = 15, units = "cm", dpi = 600)
ggsave(file.path(out_dir_d50_QC, "diff132_d50_noQC_PCA_byGenotype.pdf"), d50plot1_noQC, width = 15, height = 15, units = "cm", dpi = 600)
ggsave(file.path(out_dir_d50_QC, "diff132_d50_noQC_PCA_byCelltype.pdf"), d50plot2_noQC, width = 15, height = 15, units = "cm", dpi = 600)
ggsave(file.path(out_dir_d50_QC, "diff132_d50_noQC_PCA_byDay.pdf"), d50plot3_noQC, width = 15, height = 15, units = "cm", dpi = 600)
# --- export source data csvs
write.csv(pc_df, file.path(out_dir_d50_QC, "diff132_d50_PCA_allsamples_sourcedata.csv"), row.names = FALSE)
write.csv(pc_df_noqc, file.path(out_dir_d50_QC, "diff132_d50_PCA_noQC_sourcedata.csv"), row.names = FALSE)
```
# QC plot: Protein ID per sample x neuron
``` {r}
# --------------------------------------------------- #
# QC plot - Protein ID per sample x neuron
# --------------------------------------------------- #
n_genotypes_noQC <- n_geno_nqc
n_neurons_noQC <- length(unique(pc_df_noqc$neuron))
n_days_noQC <- n_days_nqc
# --- Step 1. Average number of ProteinIDs per type on neuron
gene_counts_neuron <- cleaned_df %>%
dplyr::group_by(neuron) %>%
dplyr::summarise(quantified_genes = n_distinct(Genes), .groups = "drop")
neuronquant <- ggplot(gene_counts_neuron, aes(x = neuron, y = quantified_genes, fill = neuron)) +
geom_bar(stat = "identity", width = 0.6) +
geom_text(aes(label = quantified_genes), vjust = -0.5, size = 3) +
scale_fill_manual(values = c("iN" = "grey60", "iDA" = "grey80", "QC" = "grey30")) +
labs(
title = "Total Number of Quantified ProteinIDs per Neuron Type",
y = "ProteinIDs", x = "Neuron"
) +
theme_bw() +
theme(
legend.position = "none",
panel.grid = element_blank()
)
neuronquant
# --- Step 2. plot protein IDs vs per genotype / neuron type ----
genotype_colors_noQC_2 <- colorRampPalette(brewer.pal(11, "RdYlBu"))(2*n_genotypes_noQC)
gene_seq_data <- diff132_d50_PCA_df %>%
dplyr::mutate(
genotype = ifelse(is.na(genotype) & neuron == "QC", "QC", genotype),
replicate = ifelse(genotype == "QC", NA, replicate)
) %>%
dplyr::filter(genotype != "QC") %>% # <-- this line removes QC
dplyr::group_by(Genes, genotype, neuron) %>%
dplyr::summarise(
total_sequences = sum(N.Sequences, na.rm = TRUE),
.groups = "drop"
) %>%
dplyr::mutate(group = interaction(genotype, neuron))
peptide_ID_QC <- ggplot(gene_seq_data, aes(x = Genes, y = log2(total_sequences), color = group)) +
geom_point(alpha = 0.7) +
scale_color_manual(values = genotype_colors_noQC_2) +
labs(
title = "Number N.Sequences vs IDs by Genotype × Neuron",
x = "Genes",
y = "Log2 Total N.Sequences"
) +
theme_bw() +
theme(
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.y = element_text(angle = 0, hjust = 1),
panel.grid = element_blank(),
legend.position = "bottom"
)
peptide_ID_QC
# --- Step 3. plot average protein IDs per genotype / neuron type ----
# First, compute quantified genes per replicate (if replicate column exists)
# calculate the average for all QC IDs
gene_counts_per_rep <- cleaned_df %>%
dplyr::mutate(genotype = ifelse(is.na(genotype) & neuron == "QC", "QC", genotype),
replicate = ifelse(genotype == "QC", NA, replicate)) %>%
dplyr::group_by(genotype, neuron, replicate) %>%
dplyr::summarise(quantified_genes = n_distinct(Genes), .groups = "drop")
# Separate QC and non-QC
qc_summary <- gene_counts_per_rep %>%
dplyr::filter(genotype == "QC") %>%
dplyr::group_by(genotype, neuron) %>%
dplyr::summarise(quantified_genes = sum(quantified_genes), .groups = "drop")
non_qc_summary <- gene_counts_per_rep %>%
dplyr::filter(genotype != "QC") %>%
dplyr::group_by(genotype, neuron) %>%
dplyr::summarise(quantified_genes = mean(quantified_genes), .groups = "drop")
# Combine
gene_counts_summary <- dplyr::bind_rows(qc_summary, non_qc_summary)
# Then, calculate mean and SD across replicates per genotype × neuron group
gene_counts_genotype <- gene_counts_summary %>%
dplyr::group_by(genotype, neuron) %>%
dplyr::summarise(
total_genes = sum(quantified_genes),
sd_genes = ifelse(n() == 1, 0, sd(quantified_genes)),
.groups = "drop"
) %>%
dplyr::mutate(group = interaction(genotype, neuron))
# add "QC" in neuron column for color matching
gene_counts_genotype <- gene_counts_genotype %>%
dplyr::mutate(genotype = replace_na(genotype, "QC"))
# Plot with SD error bars
ProteinID_genotype_neurontype <- ggplot(gene_counts_genotype, aes(x = group, y = total_genes, fill = neuron)) +
geom_bar(stat = "identity", width = 0.6) +
geom_errorbar(aes(ymin = total_genes - sd_genes, ymax = total_genes + sd_genes), width = 0.2) +
scale_fill_manual(values = c("iN" = "grey60", "iDA" = "grey80", "QC" = "grey30")) +
labs(
title = "Number of Quantified Genes per Genotype and Neuron Type",
y = "Number of ProteinIDs", x = "Genotype × Neuron"
) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 90, hjust = 1),
panel.grid = element_blank()
)
ProteinID_genotype_neurontype
# Compute average IDs per neuron type
avg_IDs_by_neuron <- gene_counts_summary %>%
dplyr::group_by(neuron) %>%
dplyr::summarise(mean_IDs = mean(quantified_genes, na.rm = TRUE)) %>%
arrange(desc(mean_IDs))
# Plot the average per neuron type
ProteinID_avg_neuron <- ggplot(avg_IDs_by_neuron, aes(x = neuron, y = mean_IDs, fill = neuron)) +
geom_col(width = 0.6) +
theme_bw() +
theme(
panel.grid = element_blank(),
legend.position = "none",) +
labs(title = "Avg. Protein IDs per Neuron Type", y = "Mean # Protein IDs", x = "") +
scale_fill_manual(values = c("iN" = "grey60", "iDA" = "grey80", "QC" = "grey30")) +
geom_text(aes(label = round(mean_IDs, 0)), vjust = -0.5, size = 3)
# --- Step 4. Save QC plots as PDFs in homedir of the project
ggsave(file.path(out_dir_d50_QC,"diff132_d50_QC_neuronquant.pdf"), neuronquant, width = 8, height = 10, units = "cm", dpi=600)
ggsave(file.path(out_dir_d50_QC,"diff132_d50_QC_ProteinID_genotype_neurontypee.pdf"), ProteinID_genotype_neurontype, width =15, height = 15, units = "cm", dpi=600)
ggsave(file.path(out_dir_d50_QC, "diff132_d50_QC_ProteinID_avg_by_neuron.pdf"),ProteinID_avg_neuron, width =12 , height = 6, units = "cm", dpi = 600)
# --- Step 5. Combine Plots
# save three main plots for QC under each other
# Force consistent theme and remove legends
p1 <- d50plot2 + theme(legend.position = "bottom")
p2 <- d50plot3 + theme(legend.position = "bottom")
p3 <- ProteinID_genotype_neurontype + theme(legend.position = "none")
p4 <- peptide_ID_QC + theme(legend.position = "none")
# Combine layout: (p1 | p2) / p3
combined_plot_QC <- ((p1 | p2) / p3 / p4) +
plot_layout(heights = c(1, 1,0.5)) &
theme(plot.margin = margin(0.01, 10, 0.01, 10)) # Top, Right, Bottom, Left padding
# save three main plots for QC under each other
# Force consistent theme and remove legends
p5 <- d50plot2_noQC + theme(legend.position = "bottom")
p6 <- d50plot3_noQC + theme(legend.position = "bottom")
p7 <- d50plot1_noQC + theme(legend.position = "bottom")
# Combine layout: (p5 | p6) / p7
combined_plot <- ((p5 | p6) / p7) +
plot_layout(heights = c(1, 1)) &
theme(plot.margin = margin(0.01, 10, 0.01, 10)) # Top, Right, Bottom, Left padding
# View or save
combined_plot
combined_plot_QC
# Save as PDF
ggsave(file.path(out_dir_d50_QC,"diff132_d50_QC_combined_plot.pdf"), combined_plot_QC, width = 12, height = 12, device = cairo_pdf)
ggsave(file.path(out_dir_d50_QC,"diff132_d50_combined_plot.pdf"), combined_plot, width = 12, height = 12, device = cairo_pdf)
```
# QC plot: RSD plot
```{r}
# --------------------------------------------------- #
# QC plot - RSD plot
# --------------------------------------------------- #
# Load the RData file containing rsd_q1_to_q3
load("/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/datasets/DataFormated.RData")
# Preprocess the RSD data:
# - Filter only Day 30 and Day 50
# - Calculate the IQR (Q3 - Q1)
# - Define lower and upper bounds as median ± half IQR
rsd_filtered <- rsd_q1_to_q3 %>%
dplyr::mutate(
Day = sub("^(\\d+)\\..*", "\\1", day_neuron_geno),
Neuron = sub("^\\d+\\.(\\w+)\\..*", "\\1", day_neuron_geno),
Genotype = sub("^\\d+\\.\\w+\\.(.*)", "\\1", day_neuron_geno)
) %>%
dplyr::filter(Day %in% c("30", "50")) %>%
dplyr::mutate(IQR = Q3 - Q1, lower = medianRSD - IQR / 2, upper = medianRSD + IQR / 2)
# --- 1. Function to generate a single plot for one neuron type and one day
make_day_plot <- function(df, neuron_type, day_label) {
# Subset the data for the chosen neuron type and day
df_day <- df %>% dplyr::filter(Neuron == neuron_type, Day == day_label)
ggplot(df_day, aes(x = Genotype)) +
# Semi-transparent ribbon showing the IQR band
geom_ribbon(aes(ymin = lower, ymax = upper, group = 1),
fill = "grey70", alpha = 0.3) +
# Dashed lines for lower and upper bounds
geom_line(aes(y = lower, group = 1),
color = "dodgerblue2", size = 0.8, linetype = "dashed") +
geom_line(aes(y = upper, group = 1),
color = "dodgerblue2", size = 0.8, linetype = "dashed") +
# Solid line for the median RSD
geom_line(aes(y = medianRSD, group = 1),
color = "dodgerblue3", size = 1.2, linetype = "solid") +
# Points marking the median values
geom_point(aes(y = medianRSD),
color = "dodgerblue3", size = 2) +
# Titles and axis labels
labs(title = paste("Day", day_label),
x = "Genotype", y = "Median RSD ± IQR") +
# Theme adjustments for cleaner look
theme_bw(base_size = 6) +
theme(
panel.grid = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 0.5)
)
}
# --- 2. Function to combine Day 30 and Day 50 plots for one neuron type
# - Day 30 panel narrower (rel_widths = c(1,4)) since fewer genotypes
plot_qc_rsd <- function(neuron_type, data) {
p30 <- make_day_plot(data, neuron_type, "30")
p50 <- make_day_plot(data, neuron_type, "50")
plot_grid(p30, p50, nrow = 1, rel_widths = c(1, 4))
}
# --- 3. Generate row plots for iN and iDA
plot_iN <- plot_qc_rsd("iN", rsd_filtered)
plot_iDA <- plot_qc_rsd("iDA", rsd_filtered)
# --- 4. Combine iN (top) and iDA (bottom) into a single figure
combined_plot_QC <- plot_grid(plot_iN, plot_iDA, ncol = 1, labels = c("iN", "iDA"))
# --- 5. Save the combined QC plot as a PDF
ggsave(file.path(out_dir_d50_QC, "diff132_d50_QC_combined_plot.pdf"), combined_plot_QC, width = 12, height = 4, device = cairo_pdf)
# --- 6. Export soucredata
write.csv(rsd_filtered, file.path(out_dir_d50_QC, "diff132_d50_QC_combined_plot_sourcedata.csv"), row.names = FALSE)
```
# =================================================
# Module 2: Summarized Data & Subcellular Analysis
# =================================================
# Read fold-change data (day 30 & day 50)
use the full dataset for plotting volcano and heatmaps etc
each measured protein has mean ctrl, mean ko and fold-change and p/q value associated with it
probably makes sense to have it in tibble for easy data manipulation
```{r}
# --------------------------------------------------- #
# iN/iDA: Day 50 Proteomics Heatmap Pipeline
# --------------------------------------------------- #
library(ComplexHeatmap)
detach("package:ComplexHeatmap", unload = TRUE, character.only = TRUE)
library(pheatmap)
library(RColorBrewer) # For color palettes
# --- Step 1. Load and Inspect Data
# Define the path to the input CSV file
input_LSD_file <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/datasets/lsd_fc_pv.csv"
# Read the CSV file into a data frame
diff132_d50_df <- read.csv(input_LSD_file, header = TRUE, sep = ",", stringsAsFactors = FALSE)
# Inspect the structure of the data
str(diff132_d50_df)
# --- Step 2. Data Cleaning
# Remove unnecessary columns: Protein.Group (1), Protein.Name (2), First.Protein.Description (4)
diff132_d50_df <- diff132_d50_df[, -c(1, 2, 4,14,15,16,17,18,19,21)]
# Filter out rows with missing or empty 'Genes' entries
clean_d50_df <- diff132_d50_df %>%
dplyr::filter(!is.na(Genes) & Genes != "") %>%
dplyr::mutate(Genes = sub(";.*", "", Genes)) # Keep first gene symbol only
# --- Step 3. Prepare Data for Heatmap
# Create a unique identifier for each sample by combining genotype and neuron
clean_d50_df<- clean_d50_df %>%
dplyr::mutate(sample_id = paste(genotype, neuron, sep = "_"))
# Pivot the data to have genes as rows and samples as columns
heatmap_d50_data <- clean_d50_df %>%
select(Genes, sample_id, fold_change) %>%
pivot_wider(
names_from = sample_id,
values_from = fold_change,
values_fill = list(fold_change = 0),
values_fn = list(fold_change = sum)
)
# Convert the data frame to a matrix and set row names to gene names
rownames_df <- heatmap_d50_data$Genes # Extract rownames before dropping the column
heatmap_matrix_d50 <- as.matrix(heatmap_d50_data[, -1])
rownames(heatmap_matrix_d50) <- rownames_df
heatmap_matrix_d50[is.na(heatmap_matrix_d50)] <- 0
# Remove rows/columns with zero variance
heatmap_matrix_d50 <- heatmap_matrix_d50[apply(heatmap_matrix_d50, 1, var) != 0, ]
heatmap_matrix_d50 <- heatmap_matrix_d50[, apply(heatmap_matrix_d50, 2, var) != 0]
# --- Step 4. Create Annotations for Heatmap & save as csv
# Extract sample information for annotations
sample_d50_info <- clean_d50_df %>%
dplyr::select(sample_id, genotype, neuron) %>%
distinct() %>%
column_to_rownames("sample_id")
df_day30 <- clean_d50_df %>% dplyr::filter(day == 30)
df_day50 <- clean_d50_df %>% dplyr::filter(day == 50)
write.csv(df_day30, file = file.path(out_dir_d50_dataframes, "df_day30_heatmap.csv"),row.names = FALSE)
write.csv(df_day50, file = file.path(out_dir_d50_dataframes, "df_day50_heatmap.csv"),row.names = FALSE)
# --- Step 5. Define a function to generate heatmaps
generate_heatmap <- function(df, day_label) {
heatmap_data <- df %>%
dplyr::mutate(sample_id = paste(genotype, neuron, sep = "_")) %>%
dplyr::select(Genes, sample_id, fold_change) %>%
pivot_wider(
names_from = sample_id,
values_from = fold_change,
values_fill = list(fold_change = 0),
values_fn = list(fold_change = sum)
)
# Extract matrix and rownames
rownames_mat <- heatmap_data$Genes
mat <- as.matrix(heatmap_data[, -1])
rownames(mat) <- rownames_mat
storage.mode(mat) <- "numeric"
mat[!is.finite(mat)] <- 0
mat <- mat[apply(mat, 1, var) != 0, , drop = FALSE]
mat <- mat[, apply(mat, 2, var) != 0, drop = FALSE]
# Annotations
sample_info <- df %>%
dplyr::mutate(sample_id = paste(genotype, neuron, sep = "_")) %>%
dplyr::select(sample_id, genotype, neuron) %>%
distinct() %>%
column_to_rownames("sample_id")
annotation_col <- data.frame(
Genotype = sample_info$genotype,
Neuron = sample_info$neuron
)
rownames(annotation_col) <- rownames(sample_info)
genotypes <- unique(sample_info$genotype)
neurons <- unique(sample_info$neuron)
ann_colors <- list(
Genotype = setNames(colorRampPalette(brewer.pal(9, "Set1"))(length(genotypes)), genotypes),
Neuron = setNames(colorRampPalette(c("grey30", "grey60", "grey80"))(length(neurons)), neurons)
)
pheatmap::pheatmap(
mat = mat,
scale = "row",
color = colorRampPalette(rev(brewer.pal(n = 11, name = "RdYlBu")))(100),
annotation_col = annotation_col,
annotation_colors = ann_colors,
breaks = seq(-4, 4, length.out = 101),
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = FALSE,
show_colnames = TRUE,
use_raster = FALSE,
angle_col = "90",
main = paste("Day", day_label, "iN / iDA Data"),
fontsize = 10,
border_color = NA
)
}
heatmap_d30 <- generate_heatmap(df_day30, 30)
heatmap_d50 <- generate_heatmap(df_day50, 50)
# --- Step 6. save single heatmaps
# Save Day 30 heatmap
pdf(file.path(out_dir_d50_heatmaps, "heatmap_day30.pdf"), width = 3, height = 6)
grid::grid.draw(heatmap_d30$gtable)
dev.off()
# Save Day 50 heatmap
pdf(file.path(out_dir_d50_heatmaps, "heatmap_day50.pdf"), width = 6, height = 6)
grid::grid.draw(heatmap_d50$gtable)
dev.off()
# --- Step 7. Combine Plots & save as PDF
# have to use gridExtra, since do not use ggplot objects, but pheatmap() returns a grid object, not a ggplot. To arrange multiple pheatmap() plots, use gridExtra::grid.arrange() instead:
library(gridExtra)
# Define the output file path
output_file_50 <- file.path(out_dir_d50_heatmaps, "diff132_d30_d50_combined_heatmap.pdf")
# Save the heatmap to a PDF file
pdf(output_file_50, width = 40 / 2.54, height = 30 / 2.54)
gridExtra::grid.arrange(
heatmap_d30$gtable,
heatmap_d50$gtable,
ncol = 2,
widths = c(1, 3) # Adjust for 25% / 75% layout
)
dev.off()
```
# Add subcell annotation & plot basic heatmaps
add subcell annotation to df and save
sub sub-matrix for d30 and d50
export heatmap as csv file
add all colums of interest: FC, mean intensity, q-value
```{r}
# --------------------------------------------------- #
# Subcell annotation to df for QC and Organelle Violin Plots
# add annotations do df and split by day and save as .csv files
# --------------------------------------------------- #
# --- Step 0: Define contaminant list
contaminant_genes <- c(
"KRT8", "KRT19", "KRT9", "KRT10", "KRT2", "KRT18", "KRT1", "KRT5",
"COL10A1", "COL13A1", "COL14A1", "COL18A1", "COL1A1", "COL1A2",
"COL26A1", "COL2A1", "COL4A1", "COL4A2", "COL5A1", "COL6A1", "COL6A2"
)
# --- Step 1. Transform columns
clean_d50_transformed_df <- clean_d50_df %>%
dplyr::mutate(
log2_fold_change = ifelse(fold_change > 0, log2(fold_change), NA),
neg_log10_q_value = ifelse(q_value > 0, -log10(q_value), NA),
neg_log10_p_value = ifelse(p_value > 0, -log10(p_value), NA)
) %>%
dplyr::filter(!Genes %in% contaminant_genes)
# --- Step 2: Load and prepare annotation matrix (only once)
cols_to_read <- c(1,2,3,4,5,6,7,9,10,12,15,16,17,18,19,20,21,22,23,25,26,27,28,29,30,31,32,33,34)
subcell_df <- read.csv("/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/tmpneurogit/proteome/SubCellAnnotation.csv", stringsAsFactors = FALSE)[, cols_to_read]
subcell_df[] <- lapply(subcell_df, as.character)
long_subcell_df <- subcell_df %>%
tidyr::pivot_longer(
cols = everything(),
names_to = "Localization",
values_to = "Genes"
) %>%
dplyr::filter(!is.na(Genes) & Genes != "") %>%
dplyr::mutate(Genes = trimws(Genes)) %>%
distinct()
binary_matrix <- long_subcell_df %>%
dplyr::mutate(value = TRUE) %>%
pivot_wider(
names_from = Localization,
values_from = value,
values_fill = list(value = FALSE)
)
# --- Step 3: Function to process each day independently
process_day_data <- function(df, day_value, save_path) {
df_day <- df %>% dplyr::filter(day == day_value)
complete_data <- df_day %>%
dplyr::select(Genes, sample_id, log2_fold_change, mean_ctrl, mean_ko, neg_log10_q_value, neg_log10_p_value) %>%
pivot_wider(
names_from = sample_id,
values_from = c(log2_fold_change, mean_ctrl, mean_ko, neg_log10_q_value, neg_log10_p_value),
names_glue = "{sample_id}_{.value}",
values_fill = list(
log2_fold_change = 0, mean_ctrl = 0, mean_ko = 0,
neg_log10_q_value = 0, neg_log10_p_value = 0
),
values_fn = list(
log2_fold_change = sum, mean_ctrl = mean, mean_ko = mean,
neg_log10_q_value = mean, neg_log10_p_value = mean
)
)
complete_data_annotated <- complete_data %>%
left_join(binary_matrix, by = "Genes")
complete_data_annotated[is.na(complete_data_annotated)] <- 0
# Save to CSV
write.csv(complete_data_annotated, save_path, row.names = FALSE)
return(complete_data_annotated)
}
# --- Step 4: Run for Day 30 and Day 50
complete_data_annotated_day30 <- process_day_data(
clean_d50_transformed_df, 30,
file.path(out_dir_d50_dataframes, "diff132_complete_annotated_day30.csv")
)
complete_data_annotated_day50 <- process_day_data(
clean_d50_transformed_df, 50,
file.path(out_dir_d50_dataframes, "diff132_complete_annotated_day50.csv")
)
# --- Step 5: Create log2FC_d50_annotated for downstream analysis (ternary, correlation, etc.)
# Reshape to have simple KO_neuron column names with log2fc values
log2fc_cols <- grep("_log2_fold_change$", colnames(complete_data_annotated_day50), value = TRUE)
annotation_cols <- colnames(binary_matrix)[-1] # exclude "Genes"
log2FC_d50_annotated <- complete_data_annotated_day50 %>%
dplyr::select(Genes, all_of(log2fc_cols), all_of(annotation_cols))
# Rename columns: remove "_log2_fold_change" suffix
colnames(log2FC_d50_annotated) <- gsub("_log2_fold_change$", "", colnames(log2FC_d50_annotated))
# --- Step 6: count sig values
diff132_log2fc_df <- clean_d50_df %>%
dplyr::mutate(log2fc = ifelse(fold_change > 0, log2(fold_change), NA))
diff132_log2fc_df %>%
dplyr::group_by(sample_id, day, neuron, genotype) %>%
dplyr::summarise(
total = n(),
n_na = sum(is.na(p_value)),
p001 = sum(p_value < 0.001, na.rm = TRUE),
p01 = sum(p_value < 0.01 & p_value >= 0.001, na.rm = TRUE),
p05 = sum(p_value < 0.05 & p_value >= 0.01, na.rm = TRUE),
ns = sum(p_value >= 0.05, na.rm = TRUE),
q001 = sum(q_value < 0.001, na.rm = TRUE),
q01 = sum(q_value < 0.01 & q_value >= 0.001, na.rm = TRUE),
q05 = sum(q_value < 0.05 & q_value >= 0.01, na.rm = TRUE),
q_ns = sum(q_value >= 0.05, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(day, neuron, genotype) %>%
print(n = Inf)
# save as csv
write.csv(diff132_log2fc_df, file.path(out_dir_d50_dataframes, "diff132_log2fc_stats.csv"), row.names = FALSE)
# Join diff132_log2fc_df with annotations
log2fc_annotated <- diff132_log2fc_df %>%
inner_join(long_subcell_df, by = "Genes")
# Filter significant (p < 0.05), summarize mean log2fc per annotation x sample_id
sig_annotation_heatmap <- log2fc_annotated %>%
dplyr::filter(p_value < 0.05) %>%
dplyr::group_by(Localization, sample_id, day, neuron) %>%
dplyr::summarise(mean_log2fc = mean(log2fc, na.rm = TRUE), .groups = "drop")
# Pivot to matrix: rows = annotations, columns = sample_id
sig_heatmap_matrix <- sig_annotation_heatmap %>%
dplyr::select(Localization, sample_id, mean_log2fc) %>%
pivot_wider(names_from = sample_id, values_from = mean_log2fc) %>%
column_to_rownames("Localization") %>%
as.matrix()
# Order columns by day and neuron
sample_order <- sig_annotation_heatmap %>%
distinct(sample_id, day, neuron) %>%
arrange(day, neuron, sample_id) %>%
pull(sample_id)
sig_heatmap_matrix <- sig_annotation_heatmap %>%
dplyr::select(Localization, sample_id, mean_log2fc) %>%
pivot_wider(names_from = sample_id, values_from = mean_log2fc, values_fn = mean) %>%
column_to_rownames("Localization") %>%
as.matrix()
# Plot
neuron_colors <- c("iN" = "grey50", "iDA" = "grey80")
# reorder columns
sig_heatmap_matrix <- sig_heatmap_matrix[, sample_order]
# filter for only day 50
sig_annotation_heatmap <- sig_annotation_heatmap %>%
dplyr::filter(day == 50)
# --- Column annotation
col_ann <- sig_annotation_heatmap %>%
dplyr::distinct(sample_id, day, neuron) %>%
dplyr::arrange(day, neuron, sample_id) %>%
tibble::column_to_rownames("sample_id") %>%
dplyr::mutate(day = factor(day))
# --- Color scale
lim <- max(abs(sig_heatmap_matrix), na.rm = TRUE)
breaks <- seq(-lim, lim, length.out = 101)
colors <- colorRampPalette(rev(brewer.pal(11,"RdYlBu")))(100)
# --- Plot
ann_colors <- list(neuron = c("iN" = "grey50", "iDA" = "grey80"))
sig_heatmap_matrix[is.na(sig_heatmap_matrix)] <- 0
ph_sig <- pheatmap::pheatmap(sig_heatmap_matrix,
color = colors,
breaks = breaks,
annotation_col = col_ann,
annotation_colors = ann_colors,
cluster_rows = TRUE,
cluster_cols = TRUE,
clustering_method = "ward.D2",
border_color = NA,
cellheight = 8,
cellwidth = 8,
fontsize = 6,
angle_col = 90,
na_col = "grey90",
main = "Mean log2FC by Localization (p < 0.05)")
ggplot2::ggsave(file.path(out_dir_d50_heatmaps, "diff132_d50_subcell_heatmap_pval.pdf"), ph_sig, width = 14, height = 10, device = grDevices::cairo_pdf)
dev.off()
```
# Euler plots for overlap annotation
```{r}
# --------------------------------------------------- #
# Euler plots of Subcellular Annotation Overlap
# --------------------------------------------------- #
# Load eulerr if not already loaded
library(eulerr)
# -- 1.1: Euler Diagram for Top 10 Vesicular Categories --
# Define terms to exclude from the top-N plot (housekeeping or broad compartments)
excluded_terms <- c("Nucleus", "Mito", "MitoIMS", "MitoMatrix", "MitoMIM", "MitoMOM", "OXPHOS", "mtComplexI",
"Golgi", "Lysosome", "EndoLyso", "ER", "Cytoplasm", "Autophagy", "NeuroDev",
"EndoDetailLoac", "RecyclingEndosome", "EarlyEndosome", "Autophagy_det")
# Select all included terms
included_terms <- setdiff(names(subcell_df), excluded_terms)
# Create named list of gene sets from annotations
gene_sets <- subcell_df %>%
dplyr::select(all_of(included_terms)) %>%
tidyr::pivot_longer(cols = everything(), names_to = "Category", values_to = "Gene") %>%
dplyr::filter(!is.na(Gene) & Gene != "") %>%
dplyr::mutate(Gene = trimws(Gene)) %>%
dplyr::group_by(Category) %>%
dplyr::summarise(GeneSet = list(unique(Gene)), .groups = "drop") %>%
deframe()
# Limit to top 10 categories by set size
top_terms <- names(sort(sapply(gene_sets, length), decreasing = TRUE))[1:10]
gene_sets_top <- gene_sets[top_terms]
gene_sets_top <- gene_sets_top[!is.na(names(gene_sets_top)) & names(gene_sets_top) != ""]
# Fit and plot Euler diagram for top vesicular categories
fit <- euler(gene_sets_top)
pdf(file.path(out_dir_d50_euler, "subcell_eulerdiagram_top10.pdf"), width = 8, height = 8)
plot(fit,
fills = list(fill = "dodgerblue3", alpha = 0.2),
labels = list(font = 1),
main = "Euler Diagram: Top 10 Subcellular Annotations")
dev.off()
# -- 3.2: Euler Diagram for Endosomal/Lysosomal Categories --
endo_terms <- c("Endo_iN_curated_Hundley", "Lysosome", "EarlyEndosome", "RecyclingEndosome")
# Create named list of gene sets
gene_sets_endo <- subcell_df %>%
dplyr::select(all_of(endo_terms)) %>%
tidyr::pivot_longer(cols = everything(), names_to = "Category", values_to = "Gene") %>%
dplyr::filter(!is.na(Gene) & Gene != "") %>%
dplyr::mutate(Gene = trimws(Gene)) %>%
dplyr::group_by(Category) %>%
dplyr::summarise(GeneSet = list(unique(Gene)), .groups = "drop") %>%
deframe()
# Fit and plot Euler diagram for endosomal/lysosomal terms
fit_endo <- euler(gene_sets_endo)
pdf(file.path(out_dir_d50_euler, "subcell_eulerdiagram_endo-lysosome.pdf"), width = 8, height = 8)
plot(fit_endo,
fills = list(fill = "firebrick", alpha = 0.5),
labels = list(font = 1),
main = "Endosomal/Lysosomal Annotations")
dev.off()
# -- 3.3: Euler Diagram for Synaptic Categories --
synaptic_terms <- c("SynGO", "SynapseSVs", "Presynaptic", "Postsynaptic",
"SVfusion", "Svexocytosis", "Svendocytosis")
# Create named list of gene sets
gene_sets_synaptic <- subcell_df %>%
dplyr::select(all_of(synaptic_terms)) %>%
tidyr::pivot_longer(cols = everything(), names_to = "Category", values_to = "Gene") %>%
dplyr::filter(!is.na(Gene) & Gene != "") %>%
dplyr::mutate(Gene = trimws(Gene)) %>%
dplyr::group_by(Category) %>%
dplyr::summarise(GeneSet = list(unique(Gene)), .groups = "drop") %>%
deframe()
# Fit and plot Euler diagram for synaptic annotations
fit_syn <- euler(gene_sets_synaptic)
pdf(file.path(out_dir_d50_euler, "subcell_eulerdiagram_synaptic.pdf"), width = 8, height = 8)
plot(fit_syn,
fills = list(fill = "dodgerblue3", alpha = 0.5),
labels = list(font = 1),
main = "Synaptic Annotations")
dev.off()
```
# LSD protein abundance across genotypes
plot bar graph of abundance of all LSDs found in the dataset
plot a bar graph of all LSD-KO matched pairs in the dataset (iN and iDA)
plot a heatmap of all LSD proteins of interest across the dataset -> check that these are KOs
```{r}
# --------------------------------------------------- #
# Plot LSD protein abundance in their respective KOs
# barplots
# --------------------------------------------------- #
# ---- Define LSD Gene List ----
LSDgenes <- c(
"AGA", "ARSA", "ARSB", "ASAH1", "ATP13A2", "CLN3", "CLN5", "CLN6", "CLN8", "CTNS", "CTSA", "CTSD", "CTSF", "DNAJC5", "FUCA", "GAA", "GALC", "GALNS", "GBA1", "GLA", "GLB1", "GM2A", "GNPTAB", "GNPTG", "GNS", "GRN", "GUSB", "HEXA", "HEXB", "HGSNAT", "HYAL1", "IDS", "IDUA", "KCTD7", "LAMP2", "LIPA", "MAN2B1", "MANBA", "MCOLN1", "MFSD8", "NAGA", "NAGLU", "NEU1", "NPC1", "NPC2", "PPT1", "PSAP", "SCARB2", "SGSH", "SLC17A5", "SMPD1", "SUMF1", "TPP1"
)
# ---- Prepare Long Format Data ----
# Filter for LSD genes and pivot to long format
df_clean_heatmap <- clean_d50_transformed_df %>%
dplyr::filter(Genes %in% LSDgenes) %>%
dplyr::select(Genes, genotype, neuron, day, mean_ctrl, mean_ko) %>%
tidyr::pivot_longer(cols = c(mean_ctrl, mean_ko), names_to = "condition", values_to = "value") %>%
dplyr::mutate(day = as.integer(day))
# ---- Filter for Matching Gene-Genotype Combinations at Day 50 ----
# Only keep rows where the gene matches the genotype (e.g., GBA1 expression in GBA1 KO)
plot_df <- df_clean_heatmap %>%
dplyr::filter(Genes == genotype, day == 50)
# ---- Replace NAs with 0 ----
plot_df <- plot_df %>%
dplyr::mutate(value = ifelse(is.na(value), 0, value))
# ---- Create DataFrame for Annotating Undetected Proteins ----
label_df <- plot_df %>%
dplyr::group_by(Genes) %>%
dplyr::mutate(max_y = max(value, na.rm = TRUE)) %>%
dplyr::filter(condition == "mean_ko" & value == 0) %>%
dplyr::mutate(
label = "not detected",
y_label = max_y * 0.5, # Position label at 50% of max height
x_shift = as.numeric(factor(neuron)) + 0.2 # Shift x position slightly for clarity
)
# ---- Plot Bar Plot with Facets per Gene ----
LSDgene_Abundance <- ggplot(plot_df, aes(x = neuron, y = value, fill = condition)) +
geom_bar(stat = "summary", fun = mean, position = position_dodge(width = 0.6), width = 0.6) +
geom_text( # Add labels for undetected mean_ko
data = label_df,
aes(x = x_shift, y = y_label, label = label),
inherit.aes = FALSE,
angle = 90,
size = 3,
color = "red"
) +
facet_wrap(~Genes, scales = "free_y") +
scale_fill_manual(values = c("mean_ctrl" = "grey70", "mean_ko" = "dodgerblue3")) +
scale_color_manual(values = c("mean_ctrl" = "grey40", "mean_ko" = "dodgerblue4")) +
labs(
title = "LSD Protein Abundance (ctrl vs KO)",
x = "Neuron Type",
y = "Abundance (raw)",
fill = "Condition"
) +
theme_bw() +
theme(
legend.position = "bottom",
panel.grid = element_blank()
)
# ---- Plot LSD gene abundance in control background ----
plot_df_ctrl <- df_clean_heatmap %>%
dplyr::filter(condition == "mean_ctrl", day == 50) %>%
distinct(Genes, neuron, value) %>%
dplyr::mutate(value = ifelse(is.na(value), 0, value))
LSDgene_Abundance_ctrl <- ggplot(plot_df_ctrl, aes(x = reorder(Genes, -value), y = value, fill = neuron)) +
geom_bar(stat = "summary", fun = mean, position = position_dodge(width = 0.7), width = 0.7) +
scale_fill_manual(values = c("iN" = "grey70", "iDA" = "dodgerblue3")) +
labs(title = "LSD Protein Abundance in Control Neurons", x = "LSD Gene", y = "Abundance (raw)", fill = "Neuron Type") +
theme_bw() +
theme(legend.position = "bottom",
panel.grid = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6))
# ---- Display Plot ----
LSDgene_Abundance
LSDgene_Abundance_ctrl
# ---- Save Plot and Data ----
ggsave(file.path(out_dir_d50_QC,"diff132_LSDgene_Abundance.pdf"), LSDgene_Abundance, width = 12, height = 12, device = cairo_pdf)
ggsave(file.path(out_dir_d50_QC,"diff132_LSDgeneAll_Abundance_inCtrl.pdf"), LSDgene_Abundance_ctrl, width = 4, height = 4, device = cairo_pdf)
write.csv(df_clean_heatmap, file = file.path(out_dir_d50_QC, "LSD_sorted.csv"), row.names = FALSE)
# --------------------------------------------------- #
# Plot LSD protein abundance in their respective KOs
# heatmaps
# --------------------------------------------------- #
# --- 0. Define expected genotype list (from your LSD gene list)
expected_genes <- c( "ASAH1", "ATP13A2", "CLN3", "CLN5", "CLN6", "CLN8", "CTSD", "CTSF",
"DNAJC5", "GAA", "GBA1", "GRN", "HEXA", "HEXB", "LIPA", "MCOLN1",
"MFSD8", "NPC1", "NPC2", "PPT1", "PSAP", "SMPD1", "TPP1")
# --- 1. Construct all combinations of neuron × condition × gene
all_rows <- expand.grid(
neuron = c("iN", "iDA"),
condition = c("mean_ctrl", "mean_ko"),
Genes = expected_genes,
stringsAsFactors = FALSE
) %>%
dplyr::mutate(RowID = paste(condition, neuron, sep = "_"))
# --- 2. Get actual values from df_clean_heatmap
values_df <- df_clean_heatmap %>%
dplyr::filter(day == 50, Genes == genotype) %>%
dplyr::group_by(Genes, neuron, condition) %>%
dplyr::summarise(value = mean(value, na.rm = TRUE), .groups = "drop") %>%
dplyr::mutate(RowID = paste(condition, neuron, sep = "_"))
# --- 4. Build matrix
heatmap_matrix_abundance <- all_rows %>%
left_join(values_df, by = c("Genes","neuron","condition","RowID")) %>%
dplyr::select(Genes, RowID, value) %>%
pivot_wider(names_from = Genes, values_from = value) %>%
column_to_rownames("RowID") %>%
as.matrix()
# --- 5. Clip high values for better contrast at low end
cap_value <- 10000000
heatmap_matrix_capped <- pmin(heatmap_matrix_abundance, cap_value)
# --- 6. Save CSV copies (raw and capped)
readr::write_csv(
heatmap_matrix_abundance %>% as.data.frame() %>% rownames_to_column("RowID"),
file.path(out_dir_d50_QC, "LSDgene_Abundance_Heatmap_matrix_raw.csv")
)
readr::write_csv(
heatmap_matrix_capped %>% as.data.frame() %>% rownames_to_column("RowID"),
file.path(out_dir_d50_QC, "LSDgene_Abundance_Heatmap_matrix_cilpped.csv")
)
# --- 7. Plot heatmap with clipped range
min_val <- suppressWarnings(min(heatmap_matrix_capped, na.rm = TRUE))
bk <- seq(min_val, cap_value, length.out = 100)
pdf(file.path(out_dir_d50_QC, "LSDgene_Abundance_Heatmap_RdYlBu_capped.pdf"),
width = 8, height = 6)
pheatmap::pheatmap(
heatmap_matrix_capped,
color = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(length(bk) - 1),
breaks = bk,
cluster_rows = FALSE,
cluster_cols = FALSE,
na_col = "grey70",
cellheight = 15,
cellwidth = 15,
border_color = "black",
legend_labels = NA,
main = "LSD Protein Abundance (Control vs KO) at Day 50 - capped"
)
dev.off()
# --- 8. Implement log2 converted heatmap here
heatmap_matrix_log2 <- all_rows %>%
left_join(values_df, by = c("Genes","neuron","condition","RowID")) %>%
dplyr::select(Genes, RowID, value) %>%
pivot_wider(names_from = Genes, values_from = value) %>%
column_to_rownames("RowID") %>%
# log2 transform (add pseudocount to avoid log2(0))
apply(2, function(x) log2(x + 1)) %>%
as.matrix()
# --- 9. Save CSV
readr::write_csv(heatmap_matrix_log2 %>% as.data.frame() %>% rownames_to_column("RowID"), file.path(out_dir_d50_QC,"LSDgene_Abundance_log2_matrix.csv"))
# --- 10. Determine breaks for color scale
min_val_log2 <- suppressWarnings(min(heatmap_matrix_log2, na.rm = TRUE))
max_val_log2 <- suppressWarnings(max(heatmap_matrix_log2, na.rm = TRUE))
bk_log2 <- seq(min_val_log2, max_val_log2, length.out = 100)
# --- 11. Plot heatmap log2 scaled
pdf(file.path(out_dir_d50_QC, "LSDgene_Abundance_Heatmap_RdYlBu_log2.pdf"),
width = 8, height = 6)
pheatmap::pheatmap(
heatmap_matrix_log2,
color = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(length(bk_log2) - 1),
breaks = bk_log2,
cluster_rows = FALSE,
cluster_cols = FALSE,
na_col = "grey70",
cellheight = 15,
cellwidth = 15,
border_color = "black",
main = "LSD Protein Abundance (Control vs KO) at Day 50 - log2 transformed"
)
dev.off()
dev.off()
```
# Sorted heatmap of neuro-related annotations
```{r}
# --------------------------------------------------- #
# Heatmap of neuro-related annotations
# plot by genotype and split for iN and iDA
# --------------------------------------------------- #
loc_ann <- c("Lysosome","RecyclingEndosome", "Endo_iN_curated_Hundley","EarlyEndosome",
"SynapseSVs","Presynaptic","Postsynaptic",
"SVfusion","Svendocytosis","Svexocytosis","vATPase")
# --- 1. helper function for one neuron type
make_ann_heatmap <- function(df, neuron_type, loc_ann, out_file) {
long_df <- df %>%
tidyr::pivot_longer(
cols = matches(paste0("_", neuron_type, "_log2_fold_change$")),
names_to = "Condition",
values_to = "log2FC"
) %>%
dplyr::mutate(genotype = sub("_.*", "", Condition)) %>%
tidyr::pivot_longer(cols = all_of(loc_ann),
names_to = "Annotation", values_to = "is_annotated") %>%
dplyr::filter(is_annotated == 1)
ann_matrix <- long_df %>%
dplyr::group_by(genotype, Annotation) %>%
dplyr::summarise(mean_log2FC = mean(log2FC, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(names_from = Annotation, values_from = mean_log2FC) %>%
column_to_rownames("genotype")
ann_mat <- as.matrix(ann_matrix)
ann_mat <- ann_mat[, loc_ann, drop = FALSE]
ann_mat[!is.finite(ann_mat)] <- NA
ann_mat <- ann_mat[rowSums(!is.na(ann_mat)) > 0, ]
# call heatmap function and save row order from pheatmap object
sorted_heatmap_iNiDA <- pheatmap::pheatmap(
mat = ann_mat,
scale = "row",
color = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(100),
cluster_rows = TRUE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
border_color = NA,
fontsize = 6,
cellheight = 8,
cellwidth = 8,
main = paste("Mean log2FC per Genotype - ", neuron_type),
filename = out_file
)
# return row order (genotype names)
return(rownames(ann_mat)[sorted_heatmap_iNiDA$tree_row$order])
}
# --- 2. Call the function and plot heatmaps for both iN and iDA
# make_ann_heatmap(
# complete_data_annotated_day50,
# "iN", loc_ann,
# file.path(out_dir_d50_ascd_avg, "heatmap_d50_iN_rowZ_byGenotype.pdf")
# )
#
# make_ann_heatmap(
# complete_data_annotated_day50,
# "iDA", loc_ann,
# file.path(out_dir_d50_ascd_avg, "heatmap_d50_iDA_rowZ_byGenotype.pdf")
# )
row_order_iN <- make_ann_heatmap(
complete_data_annotated_day50,
"iN", loc_ann,
file.path(out_dir_d50_ascd_avg, "heatmap_d50_iN_rowZ_byGenotype.pdf")
)
row_order_iDA <- make_ann_heatmap(
complete_data_annotated_day50,
"iDA", loc_ann,
file.path(out_dir_d50_ascd_avg, "heatmap_d50_iDA_rowZ_byGenotype.pdf")
)
# Sig-filtered heatmaps (p and q versions)
# --------------------------------------------------- #
# Same as make_ann_heatmap but masks non-significant
# log2FC values to NA before computing means
# sig_type: "p" or "q"
# --------------------------------------------------- #
make_ann_heatmap_sig <- function(df, neuron_type, loc_ann, out_file,
sig_type = "q", sig_threshold = 0.05) {
val_suffix <- if (sig_type == "q") "neg_log10_q_value" else "neg_log10_p_value"
neg_log10_thresh <- -log10(sig_threshold)
long_fc <- df %>%
tidyr::pivot_longer(cols = matches(paste0("_", neuron_type, "_log2_fold_change$")),
names_to = "Condition", values_to = "log2FC") %>%
dplyr::mutate(genotype = sub("_.*", "", Condition))
long_sig <- df %>%
tidyr::pivot_longer(cols = matches(paste0("_", neuron_type, "_", val_suffix, "$")),
names_to = "Condition_sig", values_to = "sig_val") %>%
dplyr::mutate(genotype = sub("_.*", "", Condition_sig)) %>%
dplyr::select(Genes, genotype, sig_val)
long_df <- long_fc %>%
left_join(long_sig, by = c("Genes", "genotype")) %>%
dplyr::mutate(log2FC = ifelse(!is.na(sig_val) & sig_val >= neg_log10_thresh, log2FC, NA)) %>%
tidyr::pivot_longer(cols = all_of(loc_ann), names_to = "Annotation", values_to = "is_annotated") %>%
dplyr::filter(is_annotated == 1)
ann_matrix <- long_df %>%
dplyr::group_by(genotype, Annotation) %>%
dplyr::summarise(mean_log2FC = mean(log2FC, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(names_from = Annotation, values_from = mean_log2FC) %>%
column_to_rownames("genotype")
ann_mat <- as.matrix(ann_matrix)
ann_mat <- ann_mat[, loc_ann, drop = FALSE]
ann_mat[!is.finite(ann_mat)] <- NA
ann_mat <- ann_mat[rowSums(!is.na(ann_mat)) > 0, , drop = FALSE]
ph <- pheatmap::pheatmap(
mat = ann_mat,
scale = "row",
color = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(100),
cluster_rows = TRUE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
border_color = NA,
fontsize = 6,
cellheight = 8,
cellwidth = 8,
na_col = "grey90",
main = paste0("Mean log2FC (", sig_type, " < ", sig_threshold, ") - ", neuron_type),
filename = out_file
)
readr::write_csv(tibble::rownames_to_column(as.data.frame(ann_mat), "genotype"),
sub("\\.pdf$", "_sourcedata.csv", out_file))
return(rownames(ann_mat)[ph$tree_row$order])
}
# --- p-value filtered
make_ann_heatmap_sig(complete_data_annotated_day50, "iN", loc_ann,
file.path(out_dir_d50_ascd_avg, "heatmap_d50_iN_rowZ_byGenotype_psig.pdf"), sig_type = "p")
make_ann_heatmap_sig(complete_data_annotated_day50, "iDA", loc_ann,
file.path(out_dir_d50_ascd_avg, "heatmap_d50_iDA_rowZ_byGenotype_psig.pdf"), sig_type = "p")
# --- q-value filtered
make_ann_heatmap_sig(complete_data_annotated_day50, "iN", loc_ann,
file.path(out_dir_d50_ascd_avg, "heatmap_d50_iN_rowZ_byGenotype_qsig.pdf"), sig_type = "q")
make_ann_heatmap_sig(complete_data_annotated_day50, "iDA", loc_ann,
file.path(out_dir_d50_ascd_avg, "heatmap_d50_iDA_rowZ_byGenotype_qsig.pdf"), sig_type = "q")
```
# Pvalue heatmaps for annotation KO vs Control
Keep same clustering / order from mean Log2FC heatmap
```{r}
# --------------------------------------------------- #
# Heatmap of neuro-related annotations (p-values)
# Same layout as VId, but cells encode signed −log10(q)
# --------------------------------------------------- #
# --- 0. helper function for stats
combine_p_fisher <- function(p_vec) {
p_use <- pmin(pmax(p_vec, .Machine$double.xmin), 1)
stat <- -2 * sum(log(p_use))
pchisq(stat, df = 2 * length(p_use), lower.tail = FALSE)
}
# --- 1. P-value heatmap with fixed row order and CSV export
make_ann_heatmap_pvals <- function(df_day50, annot_wide, neuron_type, loc_ann, out_file,
out_csv_base = NULL, row_order = NULL) {
# 1.1 defaults
if (is.null(out_csv_base)) out_csv_base <- sub("\\.pdf$", "", out_file)
# 1.2 annotation flags in long form, keep only TRUE/1, deduplicate
ann_flags_long <- annot_wide %>%
dplyr::select(Genes, dplyr::all_of(loc_ann)) %>%
tidyr::pivot_longer(cols = dplyr::all_of(loc_ann),
names_to = "Annotation", values_to = "is_annotated") %>%
dplyr::mutate(is_annotated = is_annotated %in% c(TRUE, 1)) %>%
dplyr::filter(is_annotated) %>%
dplyr::distinct(Genes, Annotation)
# 1.3 restrict to neuron, join flags (many-to-many is expected)
df_neu <- df_day50 %>%
dplyr::filter(neuron == neuron_type) %>%
dplyr::inner_join(ann_flags_long, by = "Genes", relationship = "many-to-many")
# 1.4 combine p-values per genotype x annotation; sign by mean fold_change
comb_tbl <- df_neu %>%
dplyr::filter(is.finite(p_value)) %>%
dplyr::group_by(genotype, Annotation) %>%
dplyr::summarise(
n_genes = dplyr::n(),
mean_fc = mean(fold_change, na.rm = TRUE),
p_comb = if (dplyr::n() >= 3) combine_p_fisher(p_value) else NA_real_,
.groups = "drop"
) %>%
dplyr::mutate(q_comb = p.adjust(p_comb, method = "BH")) %>%
dplyr::mutate(score = sign(replace_na(mean_fc, 0)) * (-log10(q_comb)))
# 1.5 wide matrices in requested column order
ann_matrix <- comb_tbl %>%
dplyr::select(genotype, Annotation, score, q_comb) %>%
tidyr::pivot_wider(names_from = Annotation, values_from = c(score, q_comb)) %>%
dplyr::arrange(genotype)
score_cols <- paste0("score_", loc_ann)
q_cols <- paste0("q_comb_", loc_ann)
missing_sc <- setdiff(score_cols, names(ann_matrix))
missing_q <- setdiff(q_cols, names(ann_matrix))
if (length(missing_sc)) ann_matrix[, missing_sc] <- NA_real_
if (length(missing_q)) ann_matrix[, missing_q] <- NA_real_
score_mat <- ann_matrix %>%
dplyr::select(dplyr::all_of(score_cols)) %>%
as.matrix()
rownames(score_mat) <- ann_matrix$genotype
colnames(score_mat) <- loc_ann
q_mat <- ann_matrix %>%
dplyr::select(dplyr::all_of(q_cols)) %>%
as.matrix()
rownames(q_mat) <- ann_matrix$genotype
colnames(q_mat) <- loc_ann
# 1.6 cap Inf from q=0
if (any(is.infinite(score_mat), na.rm = TRUE)) {
finite_max <- max(score_mat[is.finite(score_mat)], na.rm = TRUE)
score_mat[is.infinite(score_mat)] <- finite_max + 1
}
# 1.7 mask non-significant (q >= 0.05) or missing to NA
mask <- (q_mat >= 0.05) | is.na(q_mat)
score_mat[mask] <- NA
# 1.8 enforce external row order exactly; do not drop rows after masking
if (!is.null(row_order)) {
idx <- match(row_order, rownames(score_mat)) # preserves supplied order
idx <- idx[!is.na(idx)]
score_mat <- score_mat[idx, , drop = FALSE]
q_mat <- q_mat[idx, colnames(score_mat), drop = FALSE]
}
# 1.9 CSV exports
comb_tbl %>%
dplyr::arrange(genotype, Annotation) %>%
write.csv(file = paste0(out_csv_base, "_", neuron_type, "_summary_long.csv"),
row.names = FALSE, quote = FALSE)
tmp_scores <- cbind(genotype = rownames(score_mat), as.data.frame(score_mat))
write.csv(tmp_scores,
file = paste0(out_csv_base, "_", neuron_type, "_heatmap_scores.csv"),
row.names = FALSE, quote = FALSE)
tmp_q <- cbind(genotype = rownames(q_mat), as.data.frame(q_mat))
write.csv(tmp_q,
file = paste0(out_csv_base, "_", neuron_type, "_heatmap_qvals.csv"),
row.names = FALSE, quote = FALSE)
# 1.10 plot (no clustering to preserve order)
piyg_cols <- colorRampPalette(brewer.pal(11, "PiYG"))(201)
pheatmap::pheatmap(
mat = score_mat,
color = piyg_cols,
na_col = "grey70",
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
border_color = NA,
fontsize = 6,
cellheight = 8,
cellwidth = 8,
main = paste("Annotation p-value (signed -log10 q) -", neuron_type),
filename = out_file
)
invisible(list(summary_long = comb_tbl, score_mat = score_mat, q_mat = q_mat))
}
# --- 2. Use the row orders returned by make_ann_heatmap() above
make_ann_heatmap_pvals(
df_day50 = df_day50,
annot_wide = complete_data_annotated_day50,
neuron_type = "iN",
loc_ann = loc_ann,
out_file = file.path(out_dir_d50_heatmaps, "heatmap_d50_iN_pvals_byGenotype.pdf"),
out_csv_base = file.path(out_dir_d50_heatmaps, "heatmap_d50"),
row_order = row_order_iN
)
make_ann_heatmap_pvals(
df_day50 = df_day50,
annot_wide = complete_data_annotated_day50,
neuron_type = "iDA",
loc_ann = loc_ann,
out_file = file.path(out_dir_d50_heatmaps, "heatmap_d50_iDA_pvals_byGenotype.pdf"),
out_csv_base = file.path(out_dir_d50_heatmaps, "heatmap_d50"),
row_order = row_order_iDA
)
```
# Stouffer's method correction
```{r}
# Weighted Stouffer’s method (combine z scores with weights like √n or 1/SE). This balances smaller and larger sets.
### Plot heatmaps of pvalues for annotation KO vs Control comparison. Keep same clustering / order from mean Log2FC heatmap in VId.
# --------------------------------------------------- #
# Heatmap of neuro-related annotations (p-values)
# Same layout as VId, but cells encode signed −log10(q)
# --------------------------------------------------- #
# --- 0. helper function for stats
combine_p_fisher <- function(p_vec) {
p_use <- pmin(pmax(p_vec, .Machine$double.xmin), 1)
stat <- -2 * sum(log(p_use))
pchisq(stat, df = 2 * length(p_use), lower.tail = FALSE)
}
# --- helper: weighted Stouffer’s method
combine_p_stouffer <- function(p_vec, w_vec = NULL) {
# clip to avoid Inf
p_use <- pmin(pmax(p_vec, .Machine$double.xmin), 1)
z_vec <- qnorm(p_use, lower.tail = FALSE)
if (is.null(w_vec)) {
w_vec <- rep(1, length(z_vec)) # equal weights if none supplied
}
z_comb <- sum(w_vec * z_vec, na.rm = TRUE) / sqrt(sum(w_vec^2, na.rm = TRUE))
pnorm(z_comb, lower.tail = FALSE)
}
# --- 1. P-value heatmap with fixed row order and CSV export
make_ann_heatmap_pvals <- function(df_day50, annot_wide, neuron_type, loc_ann, out_file,
out_csv_base = NULL, row_order = NULL) {
# 1.1 defaults
if (is.null(out_csv_base)) out_csv_base <- sub("\\.pdf$", "", out_file)
# 1.2 annotation flags in long form, keep only TRUE/1, deduplicate
ann_flags_long <- annot_wide %>%
dplyr::select(Genes, dplyr::all_of(loc_ann)) %>%
tidyr::pivot_longer(cols = dplyr::all_of(loc_ann),
names_to = "Annotation", values_to = "is_annotated") %>%
dplyr::mutate(is_annotated = is_annotated %in% c(TRUE, 1)) %>%
dplyr::filter(is_annotated) %>%
dplyr::distinct(Genes, Annotation)
# 1.3 restrict to neuron, join flags (many-to-many is expected)
df_neu <- df_day50 %>%
dplyr::filter(neuron == neuron_type) %>%
dplyr::inner_join(ann_flags_long, by = "Genes", relationship = "many-to-many")
# 1.4 combine p-values per genotype x annotation; sign by mean fold_change
comb_tbl <- df_neu %>%
dplyr::filter(is.finite(p_value)) %>%
dplyr::group_by(genotype, Annotation) %>%
dplyr::summarise(
n_genes = dplyr::n(),
mean_fc = mean(fold_change, na.rm = TRUE),
p_comb = if (n_genes >= 3) combine_p_stouffer(p_value, w_vec = sqrt(n_genes)) else NA_real_,
.groups = "drop"
) %>%
dplyr::mutate(q_comb = p.adjust(p_comb, method = "BH")) %>%
dplyr::mutate(score = sign(replace_na(mean_fc, 0)) * (-log10(q_comb)))
# 1.5 wide matrices in requested column order
ann_matrix <- comb_tbl %>%
dplyr::select(genotype, Annotation, score, q_comb) %>%
tidyr::pivot_wider(names_from = Annotation, values_from = c(score, q_comb)) %>%
dplyr::arrange(genotype)
score_cols <- paste0("score_", loc_ann)
q_cols <- paste0("q_comb_", loc_ann)
missing_sc <- setdiff(score_cols, names(ann_matrix))
missing_q <- setdiff(q_cols, names(ann_matrix))
if (length(missing_sc)) ann_matrix[, missing_sc] <- NA_real_
if (length(missing_q)) ann_matrix[, missing_q] <- NA_real_
score_mat <- ann_matrix %>%
dplyr::select(dplyr::all_of(score_cols)) %>%
as.matrix()
rownames(score_mat) <- ann_matrix$genotype
colnames(score_mat) <- loc_ann
q_mat <- ann_matrix %>%
dplyr::select(dplyr::all_of(q_cols)) %>%
as.matrix()
rownames(q_mat) <- ann_matrix$genotype
colnames(q_mat) <- loc_ann
# 1.6 cap Inf from q=0
if (any(is.infinite(score_mat), na.rm = TRUE)) {
finite_max <- max(score_mat[is.finite(score_mat)], na.rm = TRUE)
score_mat[is.infinite(score_mat)] <- finite_max + 1
}
# 1.7 mask non-significant (q >= 0.05) or missing to NA
mask <- (q_mat >= 0.05) | is.na(q_mat)
score_mat[mask] <- NA
# 1.8 enforce external row order exactly; do not drop rows after masking
if (!is.null(row_order)) {
idx <- match(row_order, rownames(score_mat)) # preserves supplied order
idx <- idx[!is.na(idx)]
score_mat <- score_mat[idx, , drop = FALSE]
q_mat <- q_mat[idx, colnames(score_mat), drop = FALSE]
}
# 1.9 CSV exports
comb_tbl %>%
dplyr::arrange(genotype, Annotation) %>%
write.csv(file = paste0(out_csv_base, "_", neuron_type, "_summary_long_Stouffer.csv"),
row.names = FALSE, quote = FALSE)
tmp_scores <- cbind(genotype = rownames(score_mat), as.data.frame(score_mat))
write.csv(tmp_scores,
file = paste0(out_csv_base, "_", neuron_type, "_heatmap_scores_Stouffer.csv"),
row.names = FALSE, quote = FALSE)
tmp_q <- cbind(genotype = rownames(q_mat), as.data.frame(q_mat))
write.csv(tmp_q,
file = paste0(out_csv_base, "_", neuron_type, "_heatmap_qvals_Stouffer.csv"),
row.names = FALSE, quote = FALSE)
# 1.10 plot (no clustering to preserve order)
piyg_cols <- colorRampPalette(brewer.pal(11, "PiYG"))(201)
pheatmap::pheatmap(
mat = score_mat,
color = piyg_cols,
na_col = "grey70",
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
border_color = NA,
fontsize = 6,
cellheight = 8,
cellwidth = 8,
main = paste("Annotation p-value (signed -log10 q) -", neuron_type),
filename = out_file
)
invisible(list(summary_long = comb_tbl, score_mat = score_mat, q_mat = q_mat))
}
# --- 2. Use the row orders returned by make_ann_heatmap() above
make_ann_heatmap_pvals(
df_day50 = df_day50,
annot_wide = complete_data_annotated_day50,
neuron_type = "iN",
loc_ann = loc_ann,
out_file = file.path(out_dir_d50_ascd_avg, "heatmap_d50_iN_pvals_byGenotype_Stouffer.pdf"),
out_csv_base = file.path(out_dir_d50_ascd_avg, "heatmap_d50"),
row_order = row_order_iN
)
make_ann_heatmap_pvals(
df_day50 = df_day50,
annot_wide = complete_data_annotated_day50,
neuron_type = "iDA",
loc_ann = loc_ann,
out_file = file.path(out_dir_d50_ascd_avg, "heatmap_d50_iDA_pvals_byGenotypeStouffer.pdf"),
out_csv_base = file.path(out_dir_d50_ascd_avg, "heatmap_d50"),
row_order = row_order_iDA
)
```
# Ternary plots from log2FC per annotation
```{r}
library(Ternary)
# --- Step 1: create out dir for function call
out_dir_d50_ternaryplots <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/TernaryPlots_log2fc"
dir.create(out_dir_d50_ternaryplots, recursive = TRUE, showWarnings = FALSE)
# --- Step 2: Define function to loop over log2Fc and calculate ternary plots
plot_ternary_log2fc <- function(log2fc_mat, disease_class_df, out_dir, disease_class_palette,
tern_classes = c("SynGO", "Mito", "EndoLyso")) {
# Create output directory
tern_out_dir <- file.path(out_dir, "plots")
dir.create(tern_out_dir, recursive = TRUE, showWarnings = FALSE)
tern_suffix <- paste(tern_classes, collapse = "_")
# Identify KO columns (format: GENE_neuron, exclude annotation cols and _mean_ctrl)
ko_cols <- grep("^[A-Z0-9]+_(iN|iDA)$", colnames(log2fc_mat), value = TRUE)
ko_cols <- ko_cols[!grepl("_mean_ctrl", ko_cols)]
# Check that all tern_classes exist as columns
missing <- setdiff(tern_classes, colnames(log2fc_mat))
if (length(missing) > 0) stop("Missing annotation columns: ", paste(missing, collapse = ", "))
# Prepare disease class mapping
disease_map <- disease_class_df |>
dplyr::mutate(sample = toupper(sample)) |>
dplyr::select(sample, DiseaseClass) |>
dplyr::distinct()
# For each KO column, compute mean log2FC for genes in each annotation
results <- lapply(ko_cols, function(ko_col) {
parts <- strsplit(ko_col, "_")[[1]]
ko_gene <- parts[1]
neuron <- parts[2]
# Get log2FC values and annotation membership
df_sub <- log2fc_mat[, c("Genes", ko_col, tern_classes), drop = FALSE]
colnames(df_sub)[2] <- "log2fc"
# Compute mean log2FC for genes in each annotation, separately for pos and neg
means_pos <- sapply(tern_classes, function(tc) {
in_annot <- (df_sub[[tc]] == TRUE | df_sub[[tc]] == 1) & df_sub$log2fc > 0
if (sum(in_annot, na.rm = TRUE) == 0) return(NA_real_)
mean(df_sub$log2fc[in_annot], na.rm = TRUE)
})
means_neg <- sapply(tern_classes, function(tc) {
in_annot <- (df_sub[[tc]] == TRUE | df_sub[[tc]] == 1) & df_sub$log2fc < 0
if (sum(in_annot, na.rm = TRUE) == 0) return(NA_real_)
mean(abs(df_sub$log2fc[in_annot]), na.rm = TRUE)
})
data.frame(genotype = ko_col, KO = ko_gene, neuron = neuron,
v1_pos = means_pos[1], v2_pos = means_pos[2], v3_pos = means_pos[3],
v1_neg = means_neg[1], v2_neg = means_neg[2], v3_neg = means_neg[3],
stringsAsFactors = FALSE)
}) |> dplyr::bind_rows()
# Add disease class
results <- results |>
dplyr::left_join(disease_map, by = c("KO" = "sample"), relationship = "many-to-many") |>
dplyr::filter(!is.na(DiseaseClass)) |>
dplyr::distinct(genotype, KO, neuron, v1_pos, v2_pos, v3_pos, v1_neg, v2_neg, v3_neg, DiseaseClass)
# Compute proportions for positive log2FC (accumulation)
results_pos <- results |>
dplyr::filter(!is.na(v1_pos), !is.na(v2_pos), !is.na(v3_pos)) |>
dplyr::mutate(total = v1_pos + v2_pos + v3_pos, p_1 = v1_pos / total, p_2 = v2_pos / total, p_3 = v3_pos / total) |>
dplyr::filter(total > 0)
# Compute proportions for negative log2FC (vulnerability)
results_neg <- results |>
dplyr::filter(!is.na(v1_neg), !is.na(v2_neg), !is.na(v3_neg)) |>
dplyr::mutate(total = v1_neg + v2_neg + v3_neg, p_1 = v1_neg / total, p_2 = v2_neg / total, p_3 = v3_neg / total) |>
dplyr::filter(total > 0)
message("log2FC ternary data: ", nrow(results_pos), " pos, ", nrow(results_neg), " neg KO x neuron combinations")
# Color palettes for contours
blue_contour <- function(n) colorRampPalette(c("grey90", "dodgerblue", "dodgerblue3"))(n)
red_contour <- function(n) colorRampPalette(c("grey90", "salmon", "firebrick"))(n)
# --- Helper function to plot one ternary ---
plot_ternary_single <- function(dat, suffix, title_extra, contour_col) {
if (nrow(dat) < 3) return(invisible(NULL))
# Plot 1: All genotypes, colored by disease class, dots + labels
pdf(file.path(tern_out_dir, paste0("ternary_log2fc_combined_", suffix, "_", tern_suffix, ".pdf")), width = 6, height = 6)
TernaryPlot(alab = tern_classes[1], blab = tern_classes[2], clab = tern_classes[3],
lab.cex = 0.8, grid.lines = 5, grid.lty = "dashed", grid.minor.lines = 0,
axis.cex = 0.6, main = paste0("log2FC composition - ", title_extra))
cols <- disease_class_palette[as.character(dat$DiseaseClass)]
pch_vec <- ifelse(dat$neuron == "iN", 19, 17)
TernaryPoints(dat[, c("p_1", "p_2", "p_3")], col = cols, pch = pch_vec, cex = 1.2)
TernaryText(dat[, c("p_1", "p_2", "p_3")], labels = dat$genotype, col = "black", cex = 0.4, pos = 3, offset = 0.3)
legend("topright", legend = names(disease_class_palette), fill = disease_class_palette, cex = 0.5, bty = "n", title = "Disease Class")
legend("bottomright", legend = c("iN", "iDA"), pch = c(19, 17), cex = 0.6, bty = "n", title = "Neuron")
dev.off()
# Plot 2: Faceted by neuron and disease class
disease_class_counts <- dat |>
dplyr::group_by(DiseaseClass, neuron) |>
dplyr::summarise(n_geno = dplyr::n(), .groups = "drop") |>
dplyr::filter(n_geno >= 3) |>
dplyr::distinct(DiseaseClass) |>
dplyr::pull(DiseaseClass)
disease_classes <- intersect(unique(dat$DiseaseClass), disease_class_counts)
neuron_types <- c("iN", "iDA")
n_dc <- length(disease_classes)
if (n_dc > 0) {
pdf(file.path(tern_out_dir, paste0("ternary_log2fc_faceted_", suffix, "_", tern_suffix, ".pdf")), width = 3 * n_dc, height = 6)
par(mfrow = c(2, n_dc), mar = c(1, 1, 2, 1))
for (nt in neuron_types) {
for (dc in disease_classes) {
dd <- dat |> dplyr::filter(neuron == nt, DiseaseClass == dc)
TernaryPlot(alab = tern_classes[1], blab = tern_classes[2], clab = tern_classes[3],
lab.cex = 0.5, grid.lines = 4, grid.lty = "dashed", grid.minor.lines = 0,
axis.cex = 0.5, main = paste0(dc, " (", nt, ") - ", title_extra))
if (nrow(dd) >= 3) {
coords <- as.matrix(dd[, c("p_1", "p_2", "p_3")])
tryCatch({
TernaryDensityContour(coords, resolution = 50, col = contour_col, filled = TRUE, nlevels = 10, lwd = 0.25, labcex = 0.3)
}, error = function(e) NULL)
pt_col <- disease_class_palette[dc]
TernaryPoints(coords, pch = 19, cex = 0.8, col = pt_col)
TernaryText(coords, labels = dd$KO, cex = 0.5, pos = 3, offset = 0.3, col = "black")
}
}
}
dev.off()
}
}
# Generate plots for positive (accumulation) and negative (vulnerability)
plot_ternary_single(results_pos, "pos", "Accumulation (pos log2FC)", red_contour)
plot_ternary_single(results_neg, "neg", "Vulnerability (neg log2FC)", blue_contour)
# Save summary tables
write.csv(results_pos, file.path(tern_out_dir, paste0("ternary_log2fc_pos_summary_", tern_suffix, ".csv")), row.names = FALSE)
write.csv(results_neg, file.path(tern_out_dir, paste0("ternary_log2fc_neg_summary_", tern_suffix, ".csv")), row.names = FALSE)
invisible(list(pos = results_pos, neg = results_neg))
}
# --- Step 3: Call for main annotation triplet
plot_ternary_log2fc(
log2fc_mat = log2FC_d50_annotated,
disease_class_df = disease_class_df,
out_dir = out_dir_d50_ternaryplots,
disease_class_palette = disease_class_palette,
tern_classes = c("SynapseSVs", "OXPHOS", "Postsynaptic")
)
# c("SynapseSV", "OXPHOS", "Postsynaptic")
# c("Endo_iN_curated_Hundley", "OXPHOS","Lysosome")
#Sphingolipidoses iN:
# c("RecyclingEndosome", "Mito", "SynapseSVs")
#Sphingolipidoses iDA:
# c("RecyclingEndosome", "Svexocytosis", "SVfusion")
#Integral Membrane Protein Disorders iN:
# c("ER", "Autophagy", "EarlyEndosome")
#Integral Membrane Protein Disorders iDA:
# c("ER", "Presynaptic", "Postsynaptic")
#NCL iN:
# c("SVfusion", "Svendocytosis", "Lysosome")
#NCL iDA:
# c("ER", "EarlyEndosome", "Mito")
###### Ternary / Triangle Plots across all annotations (log2FC)
# Generates ternary plots showing log2FC composition across ALL ANNOTATIONS / SELECT PLURALITY
# 3-way combinations of specified annotation tags.
# Split by positive (accumulation) and negative (vulnerability) log2FC
# --- Step 4: Define function that plots for any given annotations the ternary plot
plot_all_ternary_log2fc_combinations <- function(log2fc_mat, disease_class_df, out_dir, disease_class_palette, tags = NULL) {
# --- Setup ---
if (is.null(tags)) {
ko_cols <- grep("^[A-Z0-9]+_(iN|iDA)$", colnames(log2fc_mat), value = TRUE)
tags <- setdiff(colnames(log2fc_mat), c("Genes", ko_cols))
tags <- tags[!grepl("_mean_ctrl", tags)]
}
# Create subfolder for all combinations
combo_out_dir <- file.path(out_dir, "plots/all_combinations")
dir.create(combo_out_dir, recursive = TRUE, showWarnings = FALSE)
# Generate all 3-way combinations
combos <- combn(tags, 3, simplify = FALSE)
n_combos <- length(combos)
message("Generating ", n_combos, " ternary plots from ", length(tags), " tags...")
# Identify KO columns
ko_cols <- grep("^[A-Z0-9]+_(iN|iDA)$", colnames(log2fc_mat), value = TRUE)
ko_cols <- ko_cols[!grepl("_mean_ctrl", ko_cols)]
# Prepare disease class mapping
disease_map <- disease_class_df |>
dplyr::mutate(sample = toupper(sample)) |>
dplyr::select(sample, DiseaseClass) |>
dplyr::distinct()
# Color palettes
blue_contour <- function(n) colorRampPalette(c("grey90", "dodgerblue", "dodgerblue3"))(n)
red_contour <- function(n) colorRampPalette(c("grey90", "salmon", "firebrick"))(n)
# --- Loop through each 3-way combination ---
for (i in seq_along(combos)) {
tern_classes <- combos[[i]]
tern_suffix <- paste(tern_classes, collapse = "_")
message("[", i, "/", n_combos, "] ", tern_suffix)
tryCatch({
# Check that all tern_classes exist as columns
missing <- setdiff(tern_classes, colnames(log2fc_mat))
if (length(missing) > 0) {
message(" Skipped: missing columns ", paste(missing, collapse = ", "))
next
}
# For each KO column, compute mean log2FC for genes in each annotation
results <- lapply(ko_cols, function(ko_col) {
parts <- strsplit(ko_col, "_")[[1]]
ko_gene <- parts[1]
neuron <- parts[2]
df_sub <- log2fc_mat[, c("Genes", ko_col, tern_classes), drop = FALSE]
colnames(df_sub)[2] <- "log2fc"
means_pos <- sapply(tern_classes, function(tc) {
in_annot <- (df_sub[[tc]] == TRUE | df_sub[[tc]] == 1) & df_sub$log2fc > 0
if (sum(in_annot, na.rm = TRUE) == 0) return(NA_real_)
mean(df_sub$log2fc[in_annot], na.rm = TRUE)
})
means_neg <- sapply(tern_classes, function(tc) {
in_annot <- (df_sub[[tc]] == TRUE | df_sub[[tc]] == 1) & df_sub$log2fc < 0
if (sum(in_annot, na.rm = TRUE) == 0) return(NA_real_)
mean(abs(df_sub$log2fc[in_annot]), na.rm = TRUE)
})
data.frame(genotype = ko_col, KO = ko_gene, neuron = neuron,
v1_pos = means_pos[1], v2_pos = means_pos[2], v3_pos = means_pos[3],
v1_neg = means_neg[1], v2_neg = means_neg[2], v3_neg = means_neg[3],
stringsAsFactors = FALSE)
}) |> dplyr::bind_rows()
# Add disease class
results <- results |>
dplyr::left_join(disease_map, by = c("KO" = "sample"), relationship = "many-to-many") |>
dplyr::filter(!is.na(DiseaseClass)) |>
dplyr::distinct(genotype, KO, neuron, v1_pos, v2_pos, v3_pos, v1_neg, v2_neg, v3_neg, DiseaseClass)
# Compute proportions for pos and neg
results_pos <- results |>
dplyr::filter(!is.na(v1_pos), !is.na(v2_pos), !is.na(v3_pos)) |>
dplyr::mutate(total = v1_pos + v2_pos + v3_pos, p_1 = v1_pos / total, p_2 = v2_pos / total, p_3 = v3_pos / total) |>
dplyr::filter(total > 0)
results_neg <- results |>
dplyr::filter(!is.na(v1_neg), !is.na(v2_neg), !is.na(v3_neg)) |>
dplyr::mutate(total = v1_neg + v2_neg + v3_neg, p_1 = v1_neg / total, p_2 = v2_neg / total, p_3 = v3_neg / total) |>
dplyr::filter(total > 0)
# Skip if insufficient data
if (nrow(results_pos) < 3 && nrow(results_neg) < 3) {
message(" Skipped: insufficient data")
next
}
# --- Helper to generate one faceted plot ---
plot_faceted <- function(dat, suffix, title_extra, contour_col) {
if (nrow(dat) < 3) return(invisible(NULL))
disease_class_counts <- dat |>
dplyr::group_by(DiseaseClass, neuron) |>
dplyr::summarise(n_geno = dplyr::n(), .groups = "drop") |>
dplyr::filter(n_geno >= 3) |>
dplyr::distinct(DiseaseClass) |>
dplyr::pull(DiseaseClass)
disease_classes <- intersect(unique(dat$DiseaseClass), disease_class_counts)
neuron_types <- c("iN", "iDA")
n_dc <- length(disease_classes)
if (n_dc == 0) return(invisible(NULL))
pdf(file.path(combo_out_dir, paste0("ternary_log2fc_", suffix, "_", tern_suffix, ".pdf")), width = 3 * n_dc, height = 7)
par(mfrow = c(2, n_dc), mar = c(1, 1, 2, 1), oma = c(0, 0, 3, 0))
for (nt in neuron_types) {
for (dc in disease_classes) {
dd <- dat |> dplyr::filter(neuron == nt, DiseaseClass == dc)
TernaryPlot(alab = tern_classes[1], blab = tern_classes[2], clab = tern_classes[3],
lab.cex = 0.5, grid.lines = 4, grid.lty = "dashed", grid.col = "skyblue",
grid.lwd = 0.25, grid.minor.lines = 0, axis.cex = 0.5,
main = paste0(dc, " (", nt, ")"))
if (nrow(dd) >= 3) {
coords <- as.matrix(dd[, c("p_1", "p_2", "p_3")])
tryCatch({
TernaryDensityContour(coords, resolution = 50, col = contour_col, filled = TRUE, nlevels = 10, lwd = 0.25, labcex = 0.3)
}, error = function(e) NULL)
pt_col <- disease_class_palette[dc]
TernaryPoints(coords, pch = 19, cex = 0.8, col = pt_col)
TernaryText(coords, labels = dd$genotype, cex = 0.5, pos = 3, offset = 0.3, col = "black")
}
}
}
mtext(paste0("Each KO positioned by % |log2FC| - ", title_extra), outer = TRUE, cex = 0.6, line = 0.5)
dev.off()
}
# Generate both pos and neg plots
plot_faceted(results_pos, "pos", "Accumulation (pos log2FC)", red_contour)
plot_faceted(results_neg, "neg", "Vulnerability (neg log2FC)", blue_contour)
}, error = function(e) {
message(" Error: ", e$message)
})
}
message("Done. ", n_combos, " combinations processed. Output: ", combo_out_dir)
}
# --- Step 5: Call function to iterate over the plotting call
plot_all_ternary_log2fc_combinations(
log2fc_mat = log2FC_d50_annotated,
disease_class_df = disease_class_df,
out_dir = out_dir_d50_ternaryplots,
disease_class_palette = disease_class_palette,
tags = c("Golgi", "EarlyEndosome", "RecyclingEndosome", "Endo_iN_curated_Hundley", "OXPHOS", "MitoIMS", "MitoMatrix", "MitoMIM", "mtComplexI", "Lysosome", "Presynaptic", "Postsynaptic", "SynapseSVs", "Svendocytosis", "Svexocytosis", "SVfusion")
)
```
# Half-circos plot of mean log2FC per annotation per KO
```{r}
# Split by neuron type (iN left, iDA right), sorted by overall mean log2FC
# Approach: duplicate data to fill full circle, then only render real sectors in top half
library(circlize)
# --- Step 1: create out dir
out_dir_d50_semicircos <- "/Users/felix/Documents/PostDoc/Harvard/03_LSD-PD/Proteomics/d50_diff132/day50/halfcircos_log2fc"
dir.create(out_dir_d50_semicircos, recursive = TRUE, showWarnings = FALSE)
# --- Step 2: Define function
plot_halfcircos_log2fc <- function(log2fc_mat, disease_class_df, out_dir,
filename = "halfcircos_log2fc",
disease_class_palette,
celltype_colors = c(iDA = "grey80", iN = "grey40"),
annotations = c("Lysosome", "Mito", "OXPHOS", "Golgi",
"Endo_iN_curated_Hundley", "EarlyEndosome",
"RecyclingEndosome", "SynapseSVs", "ER", "Autophagy")) {
# Identify KO columns
ko_cols <- grep("^[A-Z0-9]+_(iN|iDA)$", colnames(log2fc_mat), value = TRUE)
ko_cols <- ko_cols[!grepl("_mean_ctrl", ko_cols)]
# Check annotations exist
missing <- setdiff(annotations, colnames(log2fc_mat))
if (length(missing) > 0) stop("Missing annotation columns: ", paste(missing, collapse = ", "))
# Prepare disease class mapping
disease_map <- disease_class_df |>
dplyr::mutate(sample = toupper(sample)) |>
dplyr::select(sample, DiseaseClass) |>
dplyr::distinct()
# Compute mean log2FC per KO per annotation
summary_list <- lapply(ko_cols, function(ko_col) {
parts <- strsplit(ko_col, "_")[[1]]
ko_gene <- parts[1]
neuron <- parts[2]
means <- sapply(annotations, function(ann) {
in_annot <- log2fc_mat[[ann]] == TRUE
if (sum(in_annot, na.rm = TRUE) == 0) return(NA_real_)
mean(log2fc_mat[[ko_col]][in_annot], na.rm = TRUE)
})
data.frame(KO = ko_gene, neuron = neuron, sample_id = ko_col, t(means), stringsAsFactors = FALSE)
})
summary_df <- dplyr::bind_rows(summary_list)
# Add disease class
summary_df <- summary_df |>
dplyr::left_join(disease_map, by = c("KO" = "sample"), relationship = "many-to-many") |>
dplyr::filter(!is.na(DiseaseClass)) |>
dplyr::distinct(sample_id, .keep_all = TRUE)
# Compute overall mean for sorting
summary_df$overall_mean <- rowMeans(summary_df[, annotations], na.rm = TRUE)
# Save summary
write.csv(summary_df, file.path(out_dir, paste0(filename, "_summary.csv")), row.names = FALSE)
message("Summary saved: ", nrow(summary_df), " KO x neuron combinations")
# Sort: within each neuron type, sort by overall_mean
df_iN <- summary_df |> dplyr::filter(neuron == "iN") |> dplyr::arrange(desc(overall_mean))
df_iDA <- summary_df |> dplyr::filter(neuron == "iDA") |> dplyr::arrange(overall_mean)
# Combine real data - iN first (left side), then iDA (right side)
plot_df_real <- rbind(df_iN, df_iDA)
plot_df_real$sector_id <- paste0(plot_df_real$sample_id, "_", seq_len(nrow(plot_df_real)))
plot_df_real$is_dummy <- FALSE
# Create dummy data (same number of sectors, will be in bottom half - not visible)
plot_df_dummy <- plot_df_real
plot_df_dummy$sector_id <- paste0("dummy_", seq_len(nrow(plot_df_dummy)))
plot_df_dummy$is_dummy <- TRUE
# Combine: real data first (top half), then dummy (bottom half)
plot_df <- rbind(plot_df_real, plot_df_dummy)
n_iN <- nrow(df_iN)
n_iDA <- nrow(df_iDA)
n_real <- nrow(plot_df_real)
n_total <- nrow(plot_df)
# Color scale for log2FC (spectral)
log2fc_range <- range(unlist(plot_df[, annotations]), na.rm = TRUE)
log2fc_max <- max(abs(log2fc_range))
#spectral_cols <- colorRampPalette(rev(RColorBrewer::brewer.pal(11, "Spectral")))(100)
spectral_cols <- colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(100)
get_color <- function(val) {
if (is.na(val)) return("grey70")
idx <- round((val + log2fc_max) / (2 * log2fc_max) * 99) + 1
idx <- max(1, min(100, idx))
spectral_cols[idx]
}
# --- Plot ---
pdf(file.path(out_dir, paste0(filename, ".pdf")), width = 14, height = 20)
circos.clear()
# gap.after: one value per sector
# Structure: [iN real] gap [iDA real] gap [iN dummy] [iDA dummy]
gap_after <- c(
rep(0, n_iN - 1), 5,
rep(0, n_iDA - 1), 5,
rep(0, n_iN - 1), 0,
rep(0, n_iDA - 1), 0
)
circos.par(
canvas.xlim = c(-1, 1),
canvas.ylim = c(-1, 0), # show bottom half (which contains our real data)
start.degree = 180, # flip so real data appears at top visually
gap.after = gap_after,
cell.padding = c(0, 0, 0, 0)
)
message("Total sectors: ", n_total, ", gap.after length: ", length(gap_after))
# Initialize sectors
circos.initialize(
factors = factor(plot_df$sector_id, levels = plot_df$sector_id),
xlim = c(0, 1)
)
# Helper to check if sector is real (not dummy)
is_real_sector <- function(sector.index) {
!grepl("^dummy_", sector.index)
}
# Add "iN" and "iDA" labels at the sides
text(-1.05, -0.5, "iN", cex = 1.2, font = 2)
text(1.05, -0.5, "iDA", cex = 1.2, font = 2)
# Track 1: KO labels (outermost) - only for real sectors
circos.track(ylim = c(0, 1), track.height = 0.05, bg.border = NA, panel.fun = function(x, y) {
sector.index <- CELL_META$sector.index
if (!is_real_sector(sector.index)) return()
ko_name <- plot_df$KO[plot_df$sector_id == sector.index]
circos.text(0.5, 0.5, ko_name, facing = "clockwise", niceFacing = TRUE, cex = 0.4)
})
# Track 2: Neuron type - only for real sectors
circos.track(ylim = c(0, 1), track.height = 0.04, bg.border = NA, panel.fun = function(x, y) {
sector.index <- CELL_META$sector.index
if (!is_real_sector(sector.index)) return()
neuron <- plot_df$neuron[plot_df$sector_id == sector.index]
circos.rect(0, 0, 1, 1, col = celltype_colors[neuron], border = NA)
})
# Track 3: Disease class - only for real sectors
circos.track(ylim = c(0, 1), track.height = 0.04, bg.border = NA, panel.fun = function(x, y) {
sector.index <- CELL_META$sector.index
if (!is_real_sector(sector.index)) return()
dc <- plot_df$DiseaseClass[plot_df$sector_id == sector.index]
circos.rect(0, 0, 1, 1, col = disease_class_palette[dc], border = NA)
})
# Tracks 4+: Annotation heatmaps - only for real sectors
for (ann in annotations) {
circos.track(ylim = c(0, 1), track.height = 0.05, bg.border = NA, panel.fun = function(x, y) {
sector.index <- CELL_META$sector.index
if (!is_real_sector(sector.index)) return()
val <- plot_df[[ann]][plot_df$sector_id == sector.index]
circos.rect(0, 0, 1, 1, col = get_color(val), border = NA)
})
}
# Add legends
circos.clear()
# Add track labels on the right side
track_labels <- c("Neuron", "Disease", annotations)
track_y <- seq(0.95, 0.95 - 0.07 * length(track_labels), length.out = length(track_labels))
for (i in seq_along(track_labels)) {
text(1.15, track_y[i], track_labels[i], adj = 0, cex = 0.5)
}
# Neuron legend
legend("topleft", legend = names(celltype_colors), fill = celltype_colors, title = "Neuron", cex = 0.6, bty = "n")
# Disease class legend
legend("topright", legend = names(disease_class_palette), fill = disease_class_palette, title = "Disease Class", cex = 0.5, bty = "n")
# Color scale legend
legend("bottomright", legend = c(sprintf("%.1f", log2fc_max), "0", sprintf("%.1f", -log2fc_max)),
fill = spectral_cols[c(100, 50, 1)], title = "log2FC", cex = 0.5, bty = "n")
dev.off()
message("Plot saved: ", file.path(out_dir, "halfcircos_log2fc.pdf"))
invisible(summary_df)
}
# --- Step 3: Call function
plot_halfcircos_log2fc(
log2fc_mat = log2FC_d50_annotated,
disease_class_df = disease_class_df,
out_dir = out_dir_d50_semicircos,
filename = "halfcircos_log2fc",
disease_class_palette = disease_class_palette,
celltype_colors = c(iDA = "grey80", iN = "grey40"),
annotations = c("Lysosome", "Mito", "OXPHOS", "Golgi", "Endo_iN_curated_Hundley",
"EarlyEndosome", "RecyclingEndosome", "SynapseSVs", "ER", "Autophagy")
)
## --- q-value filtered half-circos
# mask non-significant values to NA
fc_cols_circos <- grep("^[A-Z0-9]+_(iN|iDA)$", colnames(log2FC_d50_annotated), value = TRUE)
q_sig_pairs <- diff132_log2fc_df %>%
dplyr::filter(day == 50, q_value < 0.05) %>%
dplyr::select(Genes, sample_id)
log2FC_d50_qsig <- log2FC_d50_annotated
for (col in fc_cols_circos) {
sig_genes <- q_sig_pairs$Genes[q_sig_pairs$sample_id == col]
log2FC_d50_qsig[[col]][!log2FC_d50_qsig$Genes %in% sig_genes] <- NA
}
plot_halfcircos_log2fc(
log2fc_mat = log2FC_d50_qsig,
filename = "halfcircos_log2fc_qsig",
disease_class_df = disease_class_df,
out_dir = out_dir_d50_semicircos,
disease_class_palette = disease_class_palette,
celltype_colors = c(iDA = "grey80", iN = "grey40"),
annotations = c("Lysosome", "Mito", "OXPHOS", "Golgi", "Endo_iN_curated_Hundley",
"EarlyEndosome", "RecyclingEndosome", "SynapseSVs", "ER", "Autophagy")
)
```
# Barplots of # of IDs sig up / down or -0.5 log2FC per KO/neuron
```{r}
# Count genes with log2FC <= -0.5 per sample
sample_cols_lostIDs <- grep("_i(N|DA)$", colnames(log2FC_d50_annotated), value = TRUE)
lost_counts <- sapply(sample_cols_lostIDs, function(col) {
sum(log2FC_d50_annotated[[col]] <= -1, na.rm = TRUE)
})
# Create df for plotting
lostIDs_df <- data.frame(
Sample = names(lost_counts),
Count = as.numeric(lost_counts)
) %>%
dplyr::mutate(
Genotype = gsub("_i(N|DA)$", "", Sample),
Neuron = gsub(".*_(i[NDA]+)$", "\\1", Sample)
)
# Order genotypes by total lost
genotype_order <- lostIDs_df %>%
dplyr::group_by(Genotype) %>%
dplyr::summarise(Total = sum(Count), .groups = "drop") %>%
dplyr::arrange(dplyr::desc(Total)) %>%
dplyr::pull(Genotype)
lostIDs_df$Genotype <- factor(lostIDs_df$Genotype, levels = genotype_order)
celltype_colors <- c(iDA = "grey80", iN = "grey40")
p_IDlost_stack <- ggplot2::ggplot(lostIDs_df, ggplot2::aes(x = Genotype, y = Count, fill = Neuron)) +
ggplot2::geom_col() +
ggplot2::scale_fill_manual(values = celltype_colors) +
ggplot2::labs(
x = NULL,
y = "Proteins lost (log2FC =< -1)",
fill = NULL
) +
ggplot2::theme_bw(base_size = 6) +
ggplot2::theme(
panel.grid = ggplot2::element_blank(),
axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5)
)
ggplot2::ggsave(file.path(out_dir_d50_heatmaps, "lost_proteins_stacked_barplot.pdf"), p_IDlost_stack, width = 3.5, height = 3)
## based on sig levels p and q levels
# Significance thresholds
q_threshold <- -log10(0.05) # 1.3
p_threshold <- -log10(0.05) # 1.3
# Function to count significant proteins by direction
count_sig_proteins <- function(data, sig_type = "q") {
if (sig_type == "q") {
fc_cols <- grep("_log2_fold_change$", colnames(data), value = TRUE)
sig_cols <- grep("_neg_log10_q_value$", colnames(data), value = TRUE)
threshold <- q_threshold
} else {
fc_cols <- grep("_log2_fold_change$", colnames(data), value = TRUE)
sig_cols <- grep("_neg_log10_p_value$", colnames(data), value = TRUE)
threshold <- p_threshold
}
results <- list()
for (i in seq_along(fc_cols)) {
fc_col <- fc_cols[i]
genotype_neuron <- gsub("_log2_fold_change$", "", fc_col)
sig_col <- paste0(genotype_neuron, "_neg_log10_", sig_type, "_value")
if (sig_col %in% sig_cols) {
is_sig <- data[[sig_col]] >= threshold & !is.na(data[[sig_col]])
fc_vals <- data[[fc_col]]
n_up <- sum(is_sig & fc_vals > 0, na.rm = TRUE)
n_down <- sum(is_sig & fc_vals < 0, na.rm = TRUE)
neuron <- gsub(".*_(iDA|iN)$", "\\1", genotype_neuron)
genotype <- gsub("_(iDA|iN)$", "", genotype_neuron)
results[[length(results) + 1]] <- data.frame(Genotype = genotype, Neuron = neuron, Direction = "Up", Count = n_up, stringsAsFactors = FALSE)
results[[length(results) + 1]] <- data.frame(Genotype = genotype, Neuron = neuron, Direction = "Down", Count = n_down, stringsAsFactors = FALSE)
}
}
do.call(rbind, results)
}
# Count for q-values
sig_df_q <- count_sig_proteins(complete_data_annotated_day50, sig_type = "q")
sig_df_q_iDA <- sig_df_q %>% dplyr::filter(Neuron == "iDA")
sig_df_q_iN <- sig_df_q %>% dplyr::filter(Neuron == "iN")
# Count for p-values
sig_df_p <- count_sig_proteins(complete_data_annotated_day50, sig_type = "p")
sig_df_p_iDA <- sig_df_p %>% dplyr::filter(Neuron == "iDA")
sig_df_p_iN <- sig_df_p %>% dplyr::filter(Neuron == "iN")
# Order genotypes by total
genotype_order_q_iDA <- sig_df_q_iDA %>% dplyr::group_by(Genotype) %>% dplyr::summarise(Total = sum(Count), .groups = "drop") %>% dplyr::arrange(dplyr::desc(Total)) %>% dplyr::pull(Genotype)
sig_df_q_iDA$Genotype <- factor(sig_df_q_iDA$Genotype, levels = genotype_order_q_iDA)
genotype_order_q_iN <- sig_df_q_iN %>% dplyr::group_by(Genotype) %>% dplyr::summarise(Total = sum(Count), .groups = "drop") %>% dplyr::arrange(dplyr::desc(Total)) %>% dplyr::pull(Genotype)
sig_df_q_iN$Genotype <- factor(sig_df_q_iN$Genotype, levels = genotype_order_q_iN)
genotype_order_p_iDA <- sig_df_p_iDA %>% dplyr::group_by(Genotype) %>% dplyr::summarise(Total = sum(Count), .groups = "drop") %>% dplyr::arrange(dplyr::desc(Total)) %>% dplyr::pull(Genotype)
sig_df_p_iDA$Genotype <- factor(sig_df_p_iDA$Genotype, levels = genotype_order_p_iDA)
genotype_order_p_iN <- sig_df_p_iN %>% dplyr::group_by(Genotype) %>% dplyr::summarise(Total = sum(Count), .groups = "drop") %>% dplyr::arrange(dplyr::desc(Total)) %>% dplyr::pull(Genotype)
sig_df_p_iN$Genotype <- factor(sig_df_p_iN$Genotype, levels = genotype_order_p_iN)
# Colors
direction_colors_iDA <- c(Down = "#9ECAE1", Up = "#F28B72")
direction_colors_iN <- c(Down = "#2A77AD", Up = "#CC4221")
## Plots
# Plot q-value iDA
p_sig_q_iDA <- ggplot2::ggplot(sig_df_q_iDA, ggplot2::aes(x = Genotype, y = Count, fill = Direction)) +
ggplot2::geom_col() + ggplot2::scale_fill_manual(values = direction_colors_iDA) +
ggplot2::labs(x = NULL, y = "Significant proteins (q < 0.05)", fill = NULL, title = "iDA") +
ggplot2::theme_bw(base_size = 6) +
ggplot2::theme(
panel.grid = ggplot2::element_blank(),
axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5))
ggplot2::ggsave(file.path(out_dir_d50_heatmaps, "sig_proteins_qvalue_iDA_barplot.pdf"), p_sig_q_iDA, width = 3.5, height = 3)
#ggplot2::ggsave(file.path(out_dir_d50_heatmaps, "sig_proteins_qvalue_iDA_barplot.pdf"), p_sig_q_iDA, width = 3.5, height = 1.5)
# Plot q-value iN
p_sig_q_iN <- ggplot2::ggplot(sig_df_q_iN, ggplot2::aes(x = Genotype, y = Count, fill = Direction)) +
ggplot2::geom_col() + ggplot2::scale_fill_manual(values = direction_colors_iN) +
ggplot2::labs(x = NULL, y = "Significant proteins (q < 0.05)", fill = NULL, title = "iN") +
ggplot2::theme_bw(base_size = 6) +
ggplot2::theme(
panel.grid = ggplot2::element_blank(),
axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5))
ggplot2::ggsave(file.path(out_dir_d50_heatmaps, "sig_proteins_qvalue_iN_barplot.pdf"), p_sig_q_iN, width = 3.5, height = 3)
#ggplot2::ggsave(file.path(out_dir_d50_heatmaps, "sig_proteins_qvalue_iN_barplot.pdf"), p_sig_q_iN, width = 3.5, height = 1.5)
# Plot p-value iDA
p_sig_p_iDA <- ggplot2::ggplot(sig_df_p_iDA, ggplot2::aes(x = Genotype, y = Count, fill = Direction)) +
ggplot2::geom_col() + ggplot2::scale_fill_manual(values = direction_colors_iDA) +
ggplot2::labs(x = NULL, y = "Significant proteins (p < 0.05)", fill = NULL, title = "iDA") +
ggplot2::theme_bw(base_size = 6) +
ggplot2::theme(
panel.grid = ggplot2::element_blank(),
axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5))
ggplot2::ggsave(file.path(out_dir_d50_heatmaps, "sig_proteins_pvalue_iDA_barplot.pdf"), p_sig_p_iDA, width = 3.5, height = 3)
# Plot p-value iN
p_sig_p_iN <- ggplot2::ggplot(sig_df_p_iN, ggplot2::aes(x = Genotype, y = Count, fill = Direction)) +
ggplot2::geom_col() + ggplot2::scale_fill_manual(values = direction_colors_iN) +
ggplot2::labs(x = NULL, y = "Significant proteins (p < 0.05)", fill = NULL, title = "iN") +
ggplot2::theme_bw(base_size = 6) + ggplot2::theme(
panel.grid = ggplot2::element_blank(),
axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5))
ggplot2::ggsave(file.path(out_dir_d50_heatmaps, "sig_proteins_pvalue_iN_barplot.pdf"), p_sig_p_iN, width = 3.5, height = 3)
```
# Sorted lineplots of neuro-related annotations
```{r}
# --------------------------------------------------- #
# Line plots of Sph mutants of interest
# plot by genotype and split for iN and iDA
# --------------------------------------------------- #
# --- 0. Define input groups for selective highlights on plot
highlight_blue <- c("SMPD1","ASAH1","GBA1","PSAP","HEXA","HEXB")
highlight_thick <- c("SMPD1","ASAH1","GBA1","CTSD","CLN6","ATP13A2","PPT1")
# --- 1. define plotting function for user-select annotations
make_lineplot_all <- function(df, neuron_type, loc_ann, out_file) {
# 1.1 reshape to long
long_df <- df %>%
tidyr::pivot_longer(
cols = matches(paste0("_", neuron_type, "_log2_fold_change$")),
names_to = "Condition",
values_to = "log2FC"
) %>%
dplyr::mutate(genotype = sub("_.*", "", Condition))
# 1.2 attach annotations
long_df <- long_df %>%
tidyr::pivot_longer(cols = all_of(loc_ann),
names_to = "Annotation", values_to = "is_annotated") %>%
filter(is_annotated == 1)
# 1.3 summarise mean log2FC per genotype × annotation
plot_df <- long_df %>%
dplyr::group_by(genotype, Annotation) %>%
dplyr::summarise(mean_log2FC = mean(log2FC, na.rm = TRUE), .groups = "drop") %>%
dplyr::mutate(Annotation = factor(Annotation, levels = loc_ann))
# 1.4 flag highlight sets
plot_df <- plot_df %>%
dplyr::mutate(
highlight_color = case_when(
genotype %in% highlight_blue ~ "dodgerblue",
TRUE ~ "grey80"
),
highlight_size = case_when(
genotype %in% highlight_thick ~ 1.2,
TRUE ~ 0.4
),
highlight_label = ifelse(genotype %in% union(highlight_blue, highlight_thick),
genotype, NA)
)
# 1.5 plot
p <- ggplot(plot_df, aes(x = Annotation, y = mean_log2FC, group = genotype)) +
geom_line(aes(color = highlight_color, size = highlight_size)) +
scale_color_identity() +
scale_size_identity() +
geom_text_repel(
data = plot_df %>% filter(!is.na(highlight_label)) %>% group_by(genotype) %>% slice_tail(n = 1),
aes(label = highlight_label, color = highlight_color),
size = 2, nudge_x = 0.5, direction = "y", hjust = 0
) +
theme_bw(base_size = 6) +
theme(
panel.grid = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)
) +
labs(
title = paste("Mean log2FC across annotations –", neuron_type),
x = "Annotation",
y = expression("Mean log"[2]*"(KO / Ctrl)")
)
ggsave(out_file, p, width = 4, height = 4, device = cairo_pdf)
return(p)
}
# --- 2. Call function and run for both iN and iDA
p_iN <- make_lineplot_all(
complete_data_annotated_day50, "iN", loc_ann,
file.path(out_dir_d50_ascd_avg, "lineplot_allGenotypes_iN.pdf")
)
p_iDA <- make_lineplot_all(
complete_data_annotated_day50, "iDA", loc_ann,
file.path(out_dir_d50_ascd_avg, "lineplot_allGenotypes_iDA.pdf")
)
```
# Key driver of SVfusion & exocytosis in iDA
```{r}
# --------------------------------------------------- #
# Heatmap of keydrivers of synaptic anotations in iDA neurons
# --------------------------------------------------- #
# --- 1. Parameters
n_top_drivers <- 40
out_dir_key <- out_dir_d50_ascd_avg_SV
# --- 2. Genotype groups
svexo_group1 <- c("CLN6","CTSF","MFSD8","PPT1","MCOLN1","CLN3","NPC2","CLN5","GRN")
svexo_group2 <- c("ASAH1","GBA1","CLN8","NPC1","DNAJC5","TPP1","ATP13A2","SMPD1","PSAP","HEXA","HEXB","GAA","CTSD","LIPA")
svfusion_group1 <- c("CLN6","CTSF","MFSD8","PPT1","MCOLN1","CLN3","NPC2","CLN5","GRN","NPC1","DNAJC5","TPP1","GAA","CTSD","LIPA")
svfusion_group2 <- c("ASAH1","CLN8","GBA1","ATP13A2","PSAP","SMPD1","HEXA","HEXB")
# --- 3. Helper: make ranked table + heatmap for one annotation ---
make_topvar_heatmap_iDA <- function(annotation_col, group1_vec, group2_vec, n_top = 40, out_prefix = "SV") {
# 3.1 subset to annotation == 1 and iDA log2FC columns
long_df <- complete_data_annotated_day50 %>%
dplyr::filter(.data[[annotation_col]] == 1) %>%
tidyr::pivot_longer(
cols = matches("_iDA_log2_fold_change$"),
names_to = "Condition",
values_to = "log2FC"
) %>%
dplyr::mutate(genotype = sub("_.*$", "", Condition)) %>%
dplyr::select(Genes, genotype, log2FC)
# 3.2 restrict to requested genotypes that actually exist
requested_genos <- unique(c(group1_vec, group2_vec))
present_genos <- intersect(requested_genos, unique(long_df$genotype))
if (length(present_genos) == 0) {
message(sprintf("No requested genotypes found for %s. Skipping.", annotation_col))
return(invisible(NULL))
}
long_df <- long_df %>% dplyr::filter(genotype %in% present_genos)
# 3.3 variability per protein across present genotypes
var_tbl <- long_df %>%
dplyr::group_by(Genes) %>%
dplyr::summarise(
n_genotypes = n_distinct(genotype),
mean_log2FC = mean(log2FC, na.rm = TRUE),
sd_log2FC = sd(log2FC, na.rm = TRUE),
var_log2FC = var(log2FC, na.rm = TRUE),
mad_log2FC = mad(log2FC, na.rm = TRUE),
.groups = "drop"
) %>%
dplyr::filter(n_genotypes >= 2) %>%
dplyr::arrange(desc(sd_log2FC))
# write full ranking
write_csv(var_tbl, file.path(out_dir_key, paste0(out_prefix, "_top_variable_genes_iDA_ranked.csv")))
# 3.4 select top N and build matrix
top_genes <- var_tbl %>% slice_head(n = n_top) %>% pull(Genes)
wide_mat <- long_df %>%
dplyr::filter(Genes %in% top_genes) %>%
dplyr::mutate(Genes = factor(Genes, levels = top_genes)) %>%
tidyr::pivot_wider(names_from = genotype, values_from = log2FC) %>%
dplyr::arrange(Genes)
# order columns: group1 then group2
g1_order <- intersect(group1_vec, colnames(wide_mat))
g2_order <- intersect(group2_vec, colnames(wide_mat))
col_order <- c(g1_order, g2_order)
wide_mat <- wide_mat[, c("Genes", col_order), drop = FALSE]
# to matrix
mat_vals <- wide_mat %>%
column_to_rownames("Genes") %>%
as.matrix()
# save the actual heatmap matrix (gene names are row.names)
write.csv(mat_vals,
file = file.path(out_dir_key, paste0(out_prefix, "_heatmap_matrix.csv")),
row.names = TRUE)
# 3.5) heat range + palette
rng_min <- suppressWarnings(min(mat_vals, na.rm = TRUE))
rng_max <- suppressWarnings(max(mat_vals, na.rm = TRUE))
if (!is.finite(rng_min) || !is.finite(rng_max)) {
message(sprintf("All values are NA for %s. Skipping heatmap.", annotation_col))
return(invisible(NULL))
}
if (isTRUE(all.equal(rng_min, rng_max))) {
rng_min <- rng_min - 0.05
rng_max <- rng_max + 0.05
}
bk <- seq(rng_min, rng_max, length.out = 100)
pal <- colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(length(bk) - 1)
# 3.6 plot heatmap
pdf(file.path(out_dir_key, paste0(out_prefix, "_top_variable_genes_iDA_heatmap.pdf")), width = 6, height = 8)
pheatmap::pheatmap(
mat_vals,
color = pal,
breaks = bk,
cluster_rows = FALSE,
cluster_cols = TRUE,
na_col = "grey90",
border_color = NA,
show_rownames = TRUE,
show_colnames = TRUE,
fontsize = 6,
cellheight = 8,
cellwidth = 8,
main = paste0("Top variable ", out_prefix, " proteins in iDA (log2FC)")
)
dev.off()
}
# --- 4. Run for both annotations & save outputs in specified outdir
make_topvar_heatmap_iDA(
annotation_col = "Svexocytosis",
group1_vec = svexo_group1,
group2_vec = svexo_group2,
n_top = n_top_drivers,
out_prefix = "SVexocytosis"
)
make_topvar_heatmap_iDA(
annotation_col = "Svendocytosis",
group1_vec = svfusion_group1,
group2_vec = svfusion_group2,
n_top = n_top_drivers,
out_prefix = "Svendocytosis"
)
make_topvar_heatmap_iDA(
annotation_col = "vATPase",
group1_vec = svfusion_group1,
group2_vec = svfusion_group2,
n_top = n_top_drivers,
out_prefix = "vATPase"
)
make_topvar_heatmap_iDA(
annotation_col = "SVfusion",
group1_vec = svfusion_group1,
group2_vec = svfusion_group2,
n_top = n_top_drivers,
out_prefix = "SVfusion"
)
dev.off()
dev.off()
```
# PCA on iDA annotation subsets
#### PCA on iDA Svendocytosis","vATPase","SVfusion","Svexocytosis" in separate iterations
```{r}
# --------------------------------------------------- #
# --- PCA drivers for synaptic-angle annotations in iDA
# --------------------------------------------------- #
library(tidyverse)
library(ggplot2)
# --- 1. define genotypes to highlight in PCA plots => mostly Sph-Mutants
highlight_genos <- c("ASAH1","CLN8","GBA1","ATP13A2","PSAP","SMPD1","HEXA","HEXB")
# --- 2. Define output dir and also user-select annotations which IDs will be PCA run on
# plot PC1 - PC5 in all pairwise permutations --> facet-wrap them
pca_out_dir <- out_dir_d50_ascd_avg_SV
candidate_annotations <- unique(c("Svendocytosis","vATPase","SVfusion","Svexocytosis"))
annotations_to_run <- intersect(candidate_annotations, colnames(complete_data_annotated_day50))
n_pcs_plot <- 5
# --- 3. define PCA function for iDA protein abundance
run_synaptic_pca_iDA <- function(annotation_col, n_pcs = 5) {
ann_long <- complete_data_annotated_day50 %>%
dplyr::filter(.data[[annotation_col]] == 1) %>%
dplyr::select(Genes, matches("_iDA_log2_fold_change$")) %>%
tidyr::pivot_longer(cols = -Genes, names_to = "Condition", values_to = "log2FC") %>%
dplyr::mutate(genotype = sub("_.*$", "", Condition)) %>%
dplyr::select(genotype, Genes, log2FC)
if (nrow(ann_long) == 0) {
message("No rows after filtering for ", annotation_col, ". Skipping.")
return(invisible(NULL))
}
wide_geno_by_gene <- ann_long %>%
tidyr::pivot_wider(names_from = Genes, values_from = log2FC) %>%
dplyr::arrange(genotype)
geno_ids <- wide_geno_by_gene$genotype
mat <- wide_geno_by_gene %>% select(-genotype) %>% as.data.frame()
# 3.1 drop genes with any NA across genotypes
mat <- mat[, colSums(is.na(mat)) == 0, drop = FALSE]
# 3.2 drop zero-variance genes (constant across genotypes)
sds <- apply(mat, 2, sd)
const_genes <- names(sds)[!is.finite(sds) | sds == 0]
if (length(const_genes) > 0) {
# save the dropped list
readr::write_csv(tibble(Genes = const_genes),
file.path(pca_out_dir, paste0(annotation_col, "_iDA_PCA_constant_genes_removed.csv")))
mat <- mat[, setdiff(colnames(mat), const_genes), drop = FALSE]
}
# 3.3 keep only genotypes with at least one non-NA left (should already hold)
keep_rows <- rowSums(is.na(mat)) < ncol(mat)
mat <- mat[keep_rows, , drop = FALSE]
geno_ids <- geno_ids[keep_rows]
if (ncol(mat) < 2 || nrow(mat) < 2) {
message("Not enough variable genes or genotypes for PCA in ", annotation_col, ". Skipping.")
return(invisible(NULL))
}
# 3.4 center PCA matrix & scale
pca_fit <- prcomp(mat, center = TRUE, scale. = TRUE)
pcs_to_use <- min(n_pcs, ncol(pca_fit$x))
scores_df <- as.data.frame(pca_fit$x[, 1:pcs_to_use, drop = FALSE]) %>%
dplyr::mutate(genotype = geno_ids) %>%
relocate(genotype)
loadings_df <- as.data.frame(pca_fit$rotation[, 1:pcs_to_use, drop = FALSE]) %>%
rownames_to_column("Genes")
# 3.5 facet all PC pairings among first n PCs
pc_pairs <- combn(paste0("PC", 1:pcs_to_use), 2, simplify = FALSE)
plot_df <- purrr::map_dfr(pc_pairs, function(p) {
tibble(
genotype = scores_df$genotype,
x = scores_df[[p[1]]],
y = scores_df[[p[2]]],
panel = paste(p[1], "vs", p[2]),
highlight = genotype %in% highlight_genos
)
})
# 3.6 define colors for highlighed genotypes and background genotypes
geno_colors <- setNames(
ifelse(unique(plot_df$genotype) %in% highlight_genos, "dodgerblue", "grey60"),
unique(plot_df$genotype)
)
# 3.7 plot PCA facetplot
p_facets <- ggplot(plot_df, aes(x = x, y = y)) +
geom_point(aes(color = genotype), size = 1.6, alpha = 0.85) +
geom_text_repel(
data = subset(plot_df, highlight),
aes(x = x, y = y, label = genotype, color = genotype),
size = 2.5,
max.overlaps = Inf,
min.segment.length = 0,
box.padding = 0.3,
segment.color = "grey50"
) +
scale_color_manual(values = geno_colors) +
facet_wrap(~ panel, scales = "free") +
labs(x = NULL, y = NULL,
title = paste0(annotation_col, " iDA PCA: pairwise PCs (top ", pcs_to_use, ")"),
color = "Genotype") +
theme_bw() +
theme(panel.grid = element_blank())
# 3.7 save PCA scores, loadings and plot
readr::write_csv(scores_df, file.path(pca_out_dir, paste0(annotation_col, "_iDA_PCA_scores.csv")))
readr::write_csv(loadings_df, file.path(pca_out_dir, paste0(annotation_col, "_iDA_PCA_loadings.csv")))
ggsave(file.path(pca_out_dir, paste0(annotation_col, "_iDA_PCA_PCpairs_facet.pdf")), p_facets, width = 7.5, height = 6.0, units = "in")
invisible(list(scores = scores_df, loadings = loadings_df))
}
# --- 4. run for all available annotations (use walk() from purr)
walk(annotations_to_run, ~ run_synaptic_pca_iDA(.x, n_pcs = n_pcs_plot))
```
#### PCA on iDA union annotation
```{r}
highlight_genos <- c("ASAH1","CLN8","GBA1","ATP13A2","PSAP","SMPD1","HEXA","HEXB")
annotations_to_combine <- c("Svexocytosis","Svendocytosis","vATPase","SVfusion")
# --- 1. get combined gene set
combined_genes <- complete_data_annotated_day50 %>%
dplyr::filter(rowSums(select(., all_of(annotations_to_combine))) > 0) %>%
pull(Genes) %>%
unique()
# --- 2. subset iDA log2FC data
ida_df <- complete_data_annotated_day50 %>%
dplyr::filter(Genes %in% combined_genes) %>%
dplyr::select(Genes, matches("_iDA_log2_fold_change$")) %>%
tidyr::pivot_longer(cols = -Genes, names_to = "Condition", values_to = "log2FC") %>%
dplyr::mutate(genotype = sub("_.*$", "", Condition)) %>%
tidyr::pivot_wider(names_from = Genes, values_from = log2FC) %>%
dplyr::arrange(genotype)
geno_ids <- ida_df$genotype
mat <- ida_df %>% dplyr::select(-genotype) %>% as.data.frame()
# 2.1 Ensure all features are numeric
mat <- mat %>% mutate(across(everything(), as.numeric))
# 2.2 Remove genes with any NA across genotypes
mat <- mat[, colSums(is.na(mat)) == 0, drop = FALSE]
# 2.3 Calculate SDs
sds <- apply(mat, 2, sd, na.rm = TRUE)
# 2.4 Drop zero-variance (constant) genes
mat <- mat[, sds > 0, drop = FALSE]
# --- 3. run PCA
pca_combined <- prcomp(mat, center = TRUE, scale. = TRUE)
# --- 4. PC permutations plot
# 4.1 scores for genotypes
scores_df <- as.data.frame(pca_combined$x) %>%
dplyr::mutate(genotype = geno_ids) %>%
relocate(genotype)
pcs_to_use <- min(5, ncol(pca_combined$x))
pc_names <- paste0("PC", 1:pcs_to_use)
# 4.2 build all PC pair panels with x = first, y = second
pc_pairs <- combn(pc_names, 2, simplify = FALSE)
plot_df <- purrr::map_dfr(pc_pairs, function(p) {
tibble(
genotype = scores_df$genotype,
x = scores_df[[p[1]]],
y = scores_df[[p[2]]],
panel = paste0(p[1], " vs ", p[2]),
highlight = genotype %in% highlight_genos
)
})
# 4.3 color map: highlights in dodgerblue, others grey60
geno_levels <- unique(scores_df$genotype)
geno_colors <- setNames(
ifelse(geno_levels %in% highlight_genos, "dodgerblue", "grey60"),
geno_levels
)
# 4.4 define plotting function for facets
p_facets <- ggplot(plot_df, aes(x = x, y = y)) +
geom_point(aes(color = genotype), size = 1.6, alpha = 0.85) +
geom_text_repel(
data = dplyr::filter(plot_df, highlight),
aes(label = genotype, color = genotype),
size = 2.5,
max.overlaps = Inf,
min.segment.length = 0,
box.padding = 0.3,
segment.color = "grey50"
) +
scale_color_manual(values = geno_colors) +
facet_wrap(~ panel, scales = "free") +
labs(
x = NULL, y = NULL,
title = "iDA PCA on union of SV annotations (PC1–PC5 pairings)",
color = "Genotype"
) +
theme_bw() +
theme(panel.grid = element_blank())
# 4.5 save PDF
ggsave(file.path(pca_out_dir, "iDA_union_annotations_PCA_PCpairs_facet.pdf"), p_facets, width = 7.5, height = 6.0, units = "in")
# --- 5. Save main drivers for each PC behind GBA1 and ASAH1
scores_df <- as.data.frame(pca_combined$x) %>% mutate(genotype = geno_ids)
# --- 6. Call function for ASAH1, PC of interest = PC where ASAH1 has sig. score (check on PCA plots)
pc_for_asah1 <- which.max(abs(scores_df %>% dplyr::filter(genotype == "ASAH1") %>% dplyr::select(PC3)))
top_asah1_drivers <- as.data.frame(pca_combined$rotation)[, pc_for_asah1, drop = FALSE] %>%
tibble::rownames_to_column("Genes") %>%
dplyr::rename(loading = 2) %>%
dplyr::arrange(desc(abs(loading))) %>%
dplyr::slice_head(n = 20)
write_csv(top_asah1_drivers, file.path(pca_out_dir, paste0("combined_PCA_top20_drivers_", "ASAH1", "_PC", pc_for_asah1, ".csv")))
# --- 7. Repeat for GBA1
pc_for_gba1 <- which.max(abs(scores_df %>% dplyr::filter(genotype == "GBA1") %>% dplyr::select(PC3)))
top_gba1_drivers <- as.data.frame(pca_combined$rotation)[, pc_for_gba1, drop = FALSE] %>%
tibble::rownames_to_column("Genes") %>%
dplyr::rename(loading = 2) %>%
dplyr::arrange(desc(abs(loading))) %>%
dplyr::slice_head(n = 20)
write_csv(top_gba1_drivers, file.path(pca_out_dir, paste0("combined_PCA_top20_drivers_", "GBA1", "_PC", pc_for_gba1, ".csv")))
```
# Cell type comparison: iN vs iDA (ctrl, day 50)
```{r}
# -------------------------------------------------- #
# Build celltypecompare_df
# -------------------------------------------------- #
# --- Step 1: Extract ctrl neurons at day 50 from cleaned_df (raw quan)
ctrl_d50 <- cleaned_df %>%
dplyr::filter(genotype == "ctrl", day == 50, neuron %in% c("iN", "iDA")) %>%
dplyr::select(Genes, neuron, replicate, quan)
# --- Step 2: Mean quan per gene per cell type
mean_by_celltype <- ctrl_d50 %>%
dplyr::group_by(Genes, neuron) %>%
dplyr::summarise(mean_quan = mean(quan, na.rm = TRUE), .groups = "drop") %>%
tidyr::pivot_wider(names_from = neuron, values_from = mean_quan, names_prefix = "mean_")
# --- Step 3: Two-sample t-test per gene (iN replicates vs iDA replicates, raw quan)
# Only genes detected in both cell types are tested.
gene_list_ct <- intersect(
ctrl_d50 %>% dplyr::filter(neuron == "iN") %>% pull(Genes) %>% unique(),
ctrl_d50 %>% dplyr::filter(neuron == "iDA") %>% pull(Genes) %>% unique()
)
ttest_results_ct <- lapply(gene_list_ct, function(g) {
iN_vals <- ctrl_d50 %>% dplyr::filter(Genes == g, neuron == "iN") %>% pull(quan)
iDA_vals <- ctrl_d50 %>% dplyr::filter(Genes == g, neuron == "iDA") %>% pull(quan)
if (length(iN_vals) < 2 | length(iDA_vals) < 2) {
return(data.frame(Genes = g, ttest_p = NA_real_))
}
tt <- suppressWarnings(t.test(iN_vals, iDA_vals, var.equal = FALSE))
data.frame(Genes = g, ttest_p = tt$p.value)
})
ttest_ct_df <- dplyr::bind_rows(ttest_results_ct)
ttest_ct_df$ttest_padj <- p.adjust(ttest_ct_df$ttest_p, method = "BH")
# --- Step 3b: Linear model per gene – quan ~ CellType (iDA reference)
# coefficient CellTypeiN = iN - iDA effect; p-value from model summary
lm_results_ct <- lapply(gene_list_ct, function(g) {
df_lm <- ctrl_d50 %>%
dplyr::filter(Genes == g) %>%
dplyr::mutate(CellType = factor(neuron, levels = c("iDA", "iN")))
if (nrow(df_lm) < 4) return(NULL)
fit <- tryCatch(lm(quan ~ CellType, data = df_lm), error = function(e) NULL)
if (is.null(fit)) return(NULL)
s <- summary(fit)
ct <- coef(fit)["CellTypeiN"]
p <- coef(s)["CellTypeiN", "Pr(>|t|)"]
data.frame(Genes = g, lm_coef_iN = as.numeric(ct), lm_p = p)
})
lm_ct_df <- dplyr::bind_rows(lm_results_ct)
lm_ct_df$lm_padj <- p.adjust(lm_ct_df$lm_p, method = "BH")
# --- Step 4: Ratio iN / iDA (keep only proteins quantified in both; avoid Inf/NaN)
ratio_ct_df <- mean_by_celltype %>%
dplyr::filter(!is.na(mean_iN) & !is.na(mean_iDA) & mean_iN > 0 & mean_iDA > 0) %>%
dplyr::mutate(
ratio_iN_iDA = mean_iN / mean_iDA,
log2_ratio_iN_iDA = log2(mean_iN / mean_iDA)
)
# --- Step 5: Join stats, annotate significance and direction (both methods)
# method 1: ttest + log2FC >= 1
# method 2: linear model p < 0.05 + log2FC >= 1
celltypecompare_base <- ratio_ct_df %>%
dplyr::left_join(ttest_ct_df, by = "Genes") %>%
dplyr::left_join(lm_ct_df, by = "Genes") %>%
dplyr::mutate(
sig_ttest = !is.na(ttest_p) & ttest_p < 0.05,
sig_fc = abs(log2_ratio_iN_iDA) >= 1,
direction = dplyr::case_when(
sig_ttest & sig_fc & ratio_iN_iDA > 1 ~ "sig_enriched_iN",
sig_ttest & sig_fc & ratio_iN_iDA < 1 ~ "sig_enriched_iDA",
TRUE ~ "unchanged"
),
direction_lm = dplyr::case_when(
!is.na(lm_p) & lm_p < 0.05 & sig_fc & lm_coef_iN > 0 ~ "sig_enriched_iN",
!is.na(lm_p) & lm_p < 0.05 & sig_fc & lm_coef_iN < 0 ~ "sig_enriched_iDA",
TRUE ~ "unchanged"
)
)
# --- Step 6: Add binary subcellular annotation columns
celltypecompare_df <- celltypecompare_base %>%
dplyr::left_join(binary_matrix, by = "Genes")
# --- Step 7: Save dataframe
write.csv(celltypecompare_df, file = file.path(out_dir_d50_dataframes, "celltypecompare_iN_iDA.csv"), row.names = FALSE)
# -------------------------------------------------- #
# Barplot: relative enrichment per annotation group
# mean log2(iN/iDA) per subcellular annotation
# -------------------------------------------------- #
annot_cols_ct <- colnames(binary_matrix)[-1]
barplot_ct_data <- lapply(annot_cols_ct, function(ann) {
genes_in_ann <- binary_matrix %>%
dplyr::filter(.data[[ann]] == TRUE) %>%
pull(Genes)
sub_df <- celltypecompare_df %>%
dplyr::filter(Genes %in% genes_in_ann, is.finite(log2_ratio_iN_iDA))
if (nrow(sub_df) == 0) return(NULL)
data.frame(
Annotation = ann,
mean_log2ratio = mean(sub_df$log2_ratio_iN_iDA, na.rm = TRUE),
sem = sd(sub_df$log2_ratio_iN_iDA, na.rm = TRUE) / sqrt(nrow(sub_df)),
n_proteins = nrow(sub_df)
)
})
barplot_ct_data <- dplyr::bind_rows(barplot_ct_data)
barplot_ct_data <- barplot_ct_data %>%
dplyr::arrange(mean_log2ratio) %>%
dplyr::mutate(Annotation = factor(Annotation, levels = Annotation))
p_ct_barplot <- ggplot(barplot_ct_data, aes(x = Annotation, y = mean_log2ratio)) +
geom_col(fill = "grey75", width = 0.7) +
geom_errorbar(aes(ymin = mean_log2ratio - sem, ymax = mean_log2ratio + sem),
width = 0.3, linewidth = 0.3) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black", linewidth = 0.3) +
labs(x = "Annotation", y = "Mean log2(iN / iDA)",
title = "Relative Enrichment of Subcellular Annotations: iN vs iDA (ctrl, day 50)") +
theme_bw(base_size = 6) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), panel.grid = element_blank())
ggsave(file.path(out_dir_d50_QC, "celltype_barplot_annotation_iN_iDA.pdf"), p_ct_barplot, width = 6, height = 4, device = cairo_pdf)
readr::write_csv(barplot_ct_data, file.path(out_dir_d50_QC, "celltype_barplot_annotation_iN_iDA_sourcedata.csv"))
# -------------------------------------------------- #
# Upset plots
# 1. ttest + log2FC method
# 2. linear model method
# 3. combined – overlap between both methods (4 sets)
# -------------------------------------------------- #
library(UpSetR)
# gene sets: ttest + log2FC method
upset_ttest_list <- list(
unchanged = celltypecompare_df %>% dplyr::filter(direction == "unchanged") %>% pull(Genes),
sig_enriched_iN = celltypecompare_df %>% dplyr::filter(direction == "sig_enriched_iN") %>% pull(Genes),
sig_enriched_iDA = celltypecompare_df %>% dplyr::filter(direction == "sig_enriched_iDA") %>% pull(Genes)
)
# gene sets: linear model method
upset_lm_list <- list(
unchanged = celltypecompare_df %>% dplyr::filter(direction_lm == "unchanged") %>% pull(Genes),
sig_enriched_iN = celltypecompare_df %>% dplyr::filter(direction_lm == "sig_enriched_iN") %>% pull(Genes),
sig_enriched_iDA = celltypecompare_df %>% dplyr::filter(direction_lm == "sig_enriched_iDA") %>% pull(Genes)
)
# gene sets: combined (4 sets for cross-method overlap)
upset_combined_list <- list(
ttest_fc_iN = celltypecompare_df %>% dplyr::filter(direction == "sig_enriched_iN") %>% pull(Genes),
ttest_fc_iDA = celltypecompare_df %>% dplyr::filter(direction == "sig_enriched_iDA") %>% pull(Genes),
lm_iN = celltypecompare_df %>% dplyr::filter(direction_lm == "sig_enriched_iN") %>% pull(Genes),
lm_iDA = celltypecompare_df %>% dplyr::filter(direction_lm == "sig_enriched_iDA") %>% pull(Genes)
)
# plot 1: ttest + log2FC
pdf(file.path(out_dir_d50_QC, "celltype_upset_ttest_fc.pdf"), width = 5, height = 4)
UpSetR::upset(UpSetR::fromList(upset_ttest_list),
sets = c("unchanged", "sig_enriched_iN", "sig_enriched_iDA"),
order.by = "freq", mainbar.y.label = "Protein Count", sets.x.label = "Set Size",
main.bar.color = "grey75", sets.bar.color = "grey75", text.scale = 1)
grid::grid.text("iN vs iDA – t-test + log2FC filter", x = 0.65, y = 0.97,
gp = grid::gpar(fontsize = 7))
dev.off()
# plot 2: linear model
pdf(file.path(out_dir_d50_QC, "celltype_upset_lm.pdf"), width = 5, height = 4)
UpSetR::upset(UpSetR::fromList(upset_lm_list),
sets = c("unchanged", "sig_enriched_iN", "sig_enriched_iDA"),
order.by = "freq", mainbar.y.label = "Protein Count", sets.x.label = "Set Size",
main.bar.color = "grey75", sets.bar.color = "grey75", text.scale = 1)
grid::grid.text("iN vs iDA – linear model", x = 0.65, y = 0.97,
gp = grid::gpar(fontsize = 7))
dev.off()
# plot 3: combined cross-method overlap
pdf(file.path(out_dir_d50_QC, "celltype_upset_combined.pdf"), width = 6, height = 4)
UpSetR::upset(UpSetR::fromList(upset_combined_list),
sets = c("ttest_fc_iN", "ttest_fc_iDA", "lm_iN", "lm_iDA"),
order.by = "freq", mainbar.y.label = "Protein Count", sets.x.label = "Set Size",
main.bar.color = "grey75", sets.bar.color = "grey75", text.scale = 1)
grid::grid.text("iN vs iDA – method comparison (ttest+FC vs lm)", x = 0.65, y = 0.97,
gp = grid::gpar(fontsize = 7))
dev.off()
# source files
readr::write_csv(
celltypecompare_df %>% dplyr::select(Genes, mean_iN, mean_iDA, log2_ratio_iN_iDA, ttest_p, ttest_padj, direction),
file.path(out_dir_d50_QC, "celltype_upset_ttest_fc_sourcedata.csv")
)
readr::write_csv(
celltypecompare_df %>% dplyr::select(Genes, mean_iN, mean_iDA, log2_ratio_iN_iDA, lm_coef_iN, lm_p, lm_padj, direction_lm),
file.path(out_dir_d50_QC, "celltype_upset_lm_sourcedata.csv")
)
# -------------------------------------------------- #
# Stacked barplots: count + fraction per annotation
# run for both methods via helper function
# -------------------------------------------------- #
direction_colors_ct <- c(sig_enriched_iN = "grey45", unchanged = "grey80", sig_enriched_iDA = "dodgerblue3")
direction_labels_ct <- c(sig_enriched_iN = "Enriched iN", unchanged = "Unchanged", sig_enriched_iDA = "Enriched iDA")
all_dirs_ct <- c("sig_enriched_iDA", "unchanged", "sig_enriched_iN")
make_stacked_ct <- function(celltypecompare_df, direction_col, method_label, out_dir) {
stacked_list <- lapply(annot_cols_ct, function(ann) {
genes_in_ann <- binary_matrix %>% dplyr::filter(.data[[ann]] == TRUE) %>% pull(Genes)
sub_df <- celltypecompare_df %>% dplyr::filter(Genes %in% genes_in_ann)
if (nrow(sub_df) == 0) return(NULL)
n_total <- nrow(sub_df)
counts <- table(factor(sub_df[[direction_col]], levels = all_dirs_ct))
data.frame(Annotation = ann, direction = names(counts),
n = as.integer(counts), n_total = n_total,
fraction = as.numeric(counts) / n_total)
})
stacked_data <- dplyr::bind_rows(stacked_list)
annot_order <- stacked_data %>%
dplyr::filter(direction == "sig_enriched_iN") %>%
dplyr::arrange(desc(fraction)) %>%
pull(Annotation)
stacked_data <- stacked_data %>%
dplyr::mutate(Annotation = factor(Annotation, levels = annot_order),
direction = factor(direction, levels = all_dirs_ct))
# absolute barplot numbers
p_count <- ggplot(stacked_data, aes(x = Annotation, y = n, fill = direction)) +
geom_col(width = 0.95, color = "black", linewidth = 0.25) +
scale_fill_manual(values = direction_colors_ct, labels = direction_labels_ct, name = NULL) +
labs(x = NULL, y = "Protein count",
title = paste0("Cell type proteins per annotation – count (", method_label, ")")) +
theme_bw(base_size = 6) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
panel.grid = element_blank(), legend.position = "bottom")
# fraction barplot
p_frac <- ggplot(stacked_data, aes(x = Annotation, y = fraction, fill = direction)) +
geom_col(width = 0.95, color = "black", linewidth = 0.25) +
geom_text(aes(label = dplyr::if_else(n > 0, paste0(n, "\n(", round(fraction * 100, 1), "%)"), ""),
group = direction),
position = position_stack(vjust = 0.5), size = 1.8, color = "black", angle = 90) +
scale_fill_manual(values = direction_colors_ct, labels = direction_labels_ct, name = NULL) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(x = NULL, y = "Fraction of proteins",
title = paste0("Cell type proteins per annotation – fraction (", method_label, ")")) +
theme_bw(base_size = 6) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
panel.grid = element_blank(), legend.position = "bottom")
ggsave(file.path(out_dir, paste0("celltype_stacked_barplot_count_", method_label, ".pdf")), p_count, width = 7, height = 4, device = cairo_pdf)
ggsave(file.path(out_dir, paste0("celltype_stacked_barplot_fraction_", method_label, ".pdf")), p_frac, width = 7, height = 4, device = cairo_pdf)
# summary source file csv
readr::write_csv(stacked_data, file.path(out_dir, paste0("celltype_stacked_barplot_", method_label, "_sourcedata.csv")))
# gene-level source file
stacked_genes_list <- lapply(annot_cols_ct, function(ann) {
genes_in_ann <- binary_matrix %>% dplyr::filter(.data[[ann]] == TRUE) %>% pull(Genes)
celltypecompare_df %>%
dplyr::filter(Genes %in% genes_in_ann) %>%
dplyr::select(Genes, mean_iN, mean_iDA, log2_ratio_iN_iDA, ttest_p, ttest_padj,
lm_coef_iN, lm_p, lm_padj, all_of(direction_col)) %>%
dplyr::mutate(Annotation = ann)
})
stacked_genes_df <- dplyr::bind_rows(stacked_genes_list) %>%
dplyr::rename(direction = all_of(direction_col)) %>%
dplyr::select(Annotation, Genes, direction, mean_iN, mean_iDA, log2_ratio_iN_iDA,
ttest_p, ttest_padj, lm_coef_iN, lm_p, lm_padj)
readr::write_csv(stacked_genes_df, file.path(out_dir, paste0("celltype_stacked_barplot_", method_label, "_genes_sourcedata.csv")))
}
# run for ttest + log2FC method
make_stacked_ct(celltypecompare_df, "direction","ttest_fc", out_dir_d50_QC)
# run for linear model method
make_stacked_ct(celltypecompare_df, "direction_lm", "lm", out_dir_d50_QC)
# gene vectors for downstream use
genes_ct_sig_iN <- celltypecompare_df %>% dplyr::filter(direction == "sig_enriched_iN") %>% pull(Genes)
genes_ct_sig_iDA <- celltypecompare_df %>% dplyr::filter(direction == "sig_enriched_iDA") %>% pull(Genes)
# -------------------------------------------------- #
# Annotation-level analysis: pie chart + GO enrichment
# direction_lm used for significance; sig_enriched_iN and sig_enriched_iDA only
# GO background: all proteins in celltypecompare_df
# call: plot_annotation_analysis("Mito", ...) or any binary_matrix annotation
# -------------------------------------------------- #
library(clusterProfiler)
library(org.Hs.eg.db)
plot_annotation_analysis <- function(annotation, celltypecompare_df, binary_matrix, out_dir) {
# proteins in this annotation
genes_in_ann <- binary_matrix %>%
dplyr::filter(.data[[annotation]] == TRUE) %>%
pull(Genes)
# sig proteins only (direction_lm)
ann_sig_df <- celltypecompare_df %>%
dplyr::filter(Genes %in% genes_in_ann,
direction_lm %in% c("sig_enriched_iN", "sig_enriched_iDA"))
if (nrow(ann_sig_df) == 0) {
message("No significant proteins for annotation: ", annotation)
return(invisible(NULL))
}
# --- Pie chart
pie_data <- ann_sig_df %>%
dplyr::count(direction_lm) %>%
dplyr::mutate(
pct = round(n / sum(n) * 100, 1),
pie_label = paste0(n, "\n(", pct, "%)"),
direction_lm = factor(direction_lm, levels = c("sig_enriched_iN", "sig_enriched_iDA"))
)
pie_colors_ann <- c(sig_enriched_iN = "grey45", sig_enriched_iDA = "dodgerblue3")
p_pie <- ggplot(pie_data, aes(x = "", y = n, fill = direction_lm)) +
geom_col(width = 1, color = "black", linewidth = 0.25) +
coord_polar(theta = "y") +
geom_text(aes(label = pie_label), position = position_stack(vjust = 0.5),
size = 2, color = "white") +
scale_fill_manual(values = pie_colors_ann,
labels = c(sig_enriched_iN = "Enriched iN", sig_enriched_iDA = "Enriched iDA"),
name = NULL) +
labs(title = paste0(annotation, " – significant proteins (lm, log2FC \u2265 1)")) +
theme_void(base_size = 6) +
theme(legend.position = "bottom", aspect.ratio = 1)
ggsave(file.path(out_dir, paste0("celltype_pie_", annotation, "_lm.pdf")), p_pie, width = 4, height = 4, device = cairo_pdf)
readr::write_csv(pie_data %>% dplyr::select(direction_lm, n, pct), file.path(out_dir, paste0("celltype_pie_", annotation, "_lm_sourcedata.csv"))
)
# --- GO enrichment per direction
bg_genes_ann <- unique(celltypecompare_df$Genes)
run_go_ann <- function(gene_set, label) {
if (length(gene_set) < 5) {
message("Too few genes for GO (", length(gene_set), "): ", label)
return(invisible(NULL))
}
for (ont in c("BP", "MF")) {
ego <- clusterProfiler::enrichGO(
gene = gene_set,
universe = bg_genes_ann,
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
ont = ont,
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.2,
readable = TRUE
)
if (is.null(ego) || nrow(as.data.frame(ego)) == 0) {
message("No GO ", ont, " results for: ", label)
next
}
readr::write_csv(as.data.frame(ego),
file.path(out_dir, paste0("GO_", ont, "_", label, "_sourcedata.csv")))
p_go <- clusterProfiler::dotplot(ego, showCategory = 5, font.size = 6) +
scale_color_gradient(low = "dodgerblue4", high = "lightblue", name = "p.adjust") +
labs(title = paste0("GO ", ont, " | ", label)) +
theme_bw(base_size = 6) +
theme(panel.grid = element_blank(),
axis.text.y = element_text(size = 5),
aspect.ratio = 1)
ggsave(file.path(out_dir, paste0("GO_", ont, "_", label, "_dotplot.pdf")), p_go, width = 4, height = 4, device = cairo_pdf)
}
}
genes_ann_iN <- ann_sig_df %>% dplyr::filter(direction_lm == "sig_enriched_iN") %>% pull(Genes)
genes_ann_iDA <- ann_sig_df %>% dplyr::filter(direction_lm == "sig_enriched_iDA") %>% pull(Genes)
run_go_ann(genes_ann_iN, paste0(annotation, "_sig_iN"))
run_go_ann(genes_ann_iDA, paste0(annotation, "_sig_iDA"))
}
# run for selected annotations
plot_annotation_analysis("Mito", celltypecompare_df, binary_matrix, out_dir_d50_QC)
plot_annotation_analysis("SynapseSVs", celltypecompare_df, binary_matrix, out_dir_d50_QC)
```
# =================================================
# Module 3: Neuro QC & Synaptic Markers
# =================================================
Aimed to plot neuro markers to plot differences between neuron types
# Plot pre/post synaptic marker
```{r}
# --------------------------------------------------- #
# neuro QC plots
# plot pre / post synaptic markers of ctrl iN vs iDA & plot
# --------------------------------------------------- #
# --- Step 1. add annotation to cleaned_df (from PCA calculation, since has raw intensities for each gene in it)
PCAclean_d50_transformed_df <- cleaned_df %>%
dplyr::filter(!Genes %in% contaminant_genes)
# --- Step 2. Join localization annotations into the filtered df
PCAclean_d50_annotated <- PCAclean_d50_transformed_df %>%
left_join(binary_matrix, by = "Genes")
# --- Step 3. Fill NAs with FALSE for localization columns
PCAclean_d50_annotated[is.na(PCAclean_d50_annotated)] <- FALSE
# --- Step 4. Dunction for plotting proteins of Interest
plot_multiple_genes <- function(genes, data = PCAclean_d50_annotated, filename = NULL) {
gene_df <- data %>%
dplyr::filter(Genes %in% genes, genotype == "ctrl", neuron %in% c("iN", "iDA"), day == 50)
if (!is.null(filename)) {
write.csv(gene_df, file = filename, row.names = FALSE)
}
ggplot(gene_df, aes(x = neuron, y = quan, fill = neuron)) +
geom_bar(stat = "summary", fun = mean, position = position_dodge(width = 0.6), width = 0.5) +
geom_jitter(aes(color = neuron), width = 0.15, size = 1.8, alpha = 0.6) +
geom_errorbar(
stat = "summary",
fun.data = mean_se,
width = 0.2,
position = position_dodge(width = 0.6)
) +
scale_fill_manual(values = c("iN" = "grey60", "iDA" = "grey80")) +
scale_color_manual(values = c("iN" = "grey30", "iDA" = "grey10")) +
facet_wrap(~Genes, scales = "free_y") +
labs(
title = "Expression of Selected Proteins in Ctrl Neurons day 50",
x = "Neuron Type",
y = "Quantified Intensity (quan)"
) +
theme_bw() +
theme(
legend.position = "none",
panel.grid = element_blank()
)
}
# --- Step 5. Call function and plot proteins of Interest
Neuro_QC_Proteins <- plot_multiple_genes(c("SYN1", "NEFH"), filename = file.path(out_dir_d50_QC, "diff132_d50_Neuro_QC_Proteins_sourcedata.csv"))
ggsave(file.path(out_dir_d50_QC, "diff132_d50_Neuro_QC_Proteins.pdf"), Neuro_QC_Proteins, width = 4, height = 4, device = cairo_pdf)
Neuro_QC_Proteins <- plot_multiple_genes(c("TH", "KCNJ6"), filename = file.path(out_dir_d50_QC, "diff132_d50_Neuro_QC_Proteins_DA_sourcedata.csv"))
ggsave(file.path(out_dir_d50_QC, "diff132_d50_Neuro_QC_Proteins_DA.pdf"), Neuro_QC_Proteins, width = 4, height = 4, device = cairo_pdf)
```
# Neuro annotation plot by neuron x genotype
Plotted in ascending order
```{r}
# --------------------------------------------------- #
# neuro QC plots
# Function for all annotations. By neuron type, in ascending order
# --------------------------------------------------- #
# --- Define annotation list ---
# all_ann <- c("Mito", "MitoIMS", "MitoMatrix", "MitoMIM", "MitoMOM", "OXPHOS",
# "mtComplexI", "Golgi", "Lysosome", "ER", "Cytoplasm", "Nucleus",
# "Autophagy", "SynGO", "Endo_iN_curated_Hundley", "SynapseSVs",
# "Presynaptic", "Postsynaptic", "NeuroDev", "EndoLyso", "EarlyEndosome",
# "RecyclingEndosome", "SVfusion", "Svendocytosis", "Svexocytosis", "vATPase")
all_ann <- c("Svexocytosis")
# --- 1. Reshape log2FC data into long format ---
ascd_ann_long_df <- complete_data_annotated_day50 %>%
tidyr::pivot_longer(
cols = matches("_log2_fold_change$"),
names_to = "Condition",
values_to = "log2FC"
) %>%
dplyr::mutate(
neuron = ifelse(grepl("_iDA_", Condition), "iDA", "iN"),
genotype = sub("_.*", "", Condition), # everything before first "_" is genotype
group = paste(genotype, neuron, sep = "_")
) %>%
# pivot annotations
tidyr::pivot_longer(cols = all_of(all_ann), names_to = "Annotation", values_to = "is_annotated") %>%
dplyr::filter(is_annotated == 1)
# --- 2. Compute average log2FC per annotation per group ---
ascd_avg_ann_fc_summary <- ascd_ann_long_df %>%
dplyr::group_by(group, neuron, Annotation) %>%
dplyr::summarise(
avg_log2FC = mean(log2FC, na.rm = TRUE),
sd_log2FC = sd(log2FC, na.rm = TRUE),
.groups = "drop")
# --- 3. Order groups within each neuron ---
ascd_avg_ann_fc_summary <- ascd_avg_ann_fc_summary %>%
dplyr::arrange(neuron, avg_log2FC) %>%
group_by(neuron) %>%
dplyr::mutate(group = factor(group, levels = unique(group)))
# --- 4. Merge ordering back into main df ---
ascd_ann_long_df <- ascd_ann_long_df %>%
inner_join(ascd_avg_ann_fc_summary %>% dplyr::select(group, neuron),
by = c("group", "neuron")) %>%
dplyr::mutate(group = factor(group, levels = levels(ascd_avg_ann_fc_summary$group)))
# --- 5. Plot: mean ± SD from summary ---
d50_ascd_ann_plot <- ggplot(ascd_avg_ann_fc_summary, aes(x = group, y = avg_log2FC)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_linerange(aes(ymin = avg_log2FC - sd_log2FC,
ymax = avg_log2FC + sd_log2FC),
color = "grey70", linewidth = 1.2, alpha = 0.8) +
geom_line(aes(group = 1), color = "red", linewidth = 0.9) +
geom_point(color = "red", size = 0.8) +
facet_wrap(~ neuron, scales = "free_x") +
theme_bw(base_size = 6) +
theme(
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
panel.grid = element_blank()
) +
labs(
title = "Functional Annotations: mean ± SD by Genotype x Neuron",
x = "Genotype x Neuron (sorted by mean log2FC)",
y = expression(log[2]*"(KO / Ctrl)")
)
d50_ascd_ann_plot
# --- 6. Save PDF ---
ggsave(file.path(out_dir_d50_ascd_avg, "diff132_d50_ascending_Avg_Svexocytosis.pdf"), d50_ascd_ann_plot, width = 4, height = 3, device = cairo_pdf)
# --- 7. Zoom into first 6 groups per neuron ---
ascd_top6_summary <- ascd_avg_ann_fc_summary %>%
dplyr::group_by(neuron) %>%
slice_head(n = 6) %>%
ungroup() %>%
dplyr::mutate(group = factor(group, levels = levels(ascd_avg_ann_fc_summary$group)))
# --- 8. Plot top 6
d50_ascd_ann_plot_top6 <- ggplot(ascd_top6_summary, aes(x = group, y = avg_log2FC)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_linerange(aes(ymin = avg_log2FC - sd_log2FC,
ymax = avg_log2FC + sd_log2FC),
color = "grey70", linewidth = 1.2, alpha = 0.8) +
geom_line(aes(group = 1), color = "red", linewidth = 0.9) +
geom_point(color = "red", size = 0.8) +
facet_wrap(~ neuron, scales = "free_x") +
theme_bw(base_size = 6) +
theme(
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
panel.grid = element_blank()
) +
labs(
title = "Top Genotypes: mean ± SD",
x = "Genotype x Neuron (Top 6)",
y = expression(log[2]*"(KO / Ctrl)")
)
# --- 9. Save second PDF ---
ggsave(file.path(out_dir_d50_ascd_avg, "diff132_d50_Top5_Avg_Svexocytosis.pdf"), d50_ascd_ann_plot_top6, width = 2, height = 3, device = cairo_pdf)
```
# Plot differentiation markers for all KOs
```{r}
# =================================================
# Differentiation QC
# =================================================
# --- Setup output directory
out_dir_diffQC <- file.path(out_dir_d50_QC, "DifferentiationQC")
dir.create(out_dir_diffQC, recursive = TRUE, showWarnings = FALSE)
# --- Create differentiation_df from clean_d50_df
# Goal: reshape data so ctrl is a genotype alongside KOs for barplot comparison
# clean_d50_df has mean_ctrl/sd_ctrl (same for all KOs) and mean_ko/sd_ko per condition
# We extract ctrl as its own rows, then bind with KO rows
# Ctrl: take one row per gene x neuron (mean_ctrl is identical across KO rows)
ctrl_df <- clean_d50_df %>%
dplyr::filter(day == 50) %>%
distinct(Genes, neuron, .keep_all = TRUE) %>%
transmute(Genes, neuron, genotype = "ctrl", mean_intensity = mean_ctrl, sd_intensity = sd_ctrl)
# KO: rename mean_ko/sd_ko to match ctrl structure
ko_df <- clean_d50_df %>%
dplyr::filter(day == 50) %>%
transmute(Genes, neuron, genotype, mean_intensity = mean_ko, sd_intensity = sd_ko)
# Set genotype factor with ctrl first, then alphabetical
all_genotypes <- sort(unique(c("ctrl", ko_df$genotype)))
all_genotypes <- c("ctrl", all_genotypes[all_genotypes != "ctrl"])
differentiation_df <- bind_rows(ctrl_df, ko_df) %>%
dplyr::mutate(genotype = factor(genotype, levels = all_genotypes))
# --- Gene lists
neuron_drivers <- c("BDNF", "CAMK2B", "CRTC1", "DCX", "JUN", "MAP2", "NCAM1", "NEFH", "NEFL", "NEFM", "NES", "POU3F2", "SLC17A7", "SYN1", "SYP", "TUBB3", "TH", "BSN", "SYNJ1", "GAP43", "SYN2", "SYN3")
presyn_genes <- c("BSN", "CPLX1", "CPLX2", "PCLO", "RAB3A", "SNAP25", "STX1A", "STX1B", "SV2A", "SV2B", "SV2C", "SYN1", "SYN2", "SYP", "SYT1", "SYT2", "SYT7", "VAMP1", "VAMP2")
postsyn_genes <- c("DLG4", "DLG3", "DLG2", "SHANK1", "SHANK2", "SHANK3", "HOMER1", "HOMER2", "HOMER3", "GPHN", "SYNGAP1", "GRIN1", "GRIN2A", "GRIN2B", "GRIA1", "GRIA2", "GRIA3", "GRIA4", "CAMK2A", "CAMK2B", "NRGN", "NPTX2")
# --- Plot function for marker genes (individual facets)
plot_diff_markers <- function(genes, data, title_text) {
plot_df <- data %>% dplyr::filter(Genes %in% genes)
if (nrow(plot_df) == 0) return(NULL)
ggplot(plot_df, aes(x = genotype, y = mean_intensity, fill = neuron)) +
geom_bar(stat = "identity", position = position_dodge(width = 1), width = 0.7) +
geom_errorbar(aes(ymin = mean_intensity - sd_intensity, ymax = mean_intensity + sd_intensity), position = position_dodge(width = 1), width = 0.2) +
scale_fill_manual(values = c("iN" = "grey60", "iDA" = "grey80")) +
facet_wrap(~Genes, scales = "free_y") +
labs(title = title_text, x = "Genotype", y = "Mean Intensity") +
theme_bw(base_size = 6) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), panel.grid = element_blank())
}
# --- Plot 1: Neuronal markers (individual)
p_neuro <- plot_diff_markers(neuron_drivers, differentiation_df, "Neuronal Markers across all KOs")
ggsave(file.path(out_dir_diffQC, "diff132_d50_Neuronal_markers.pdf"), p_neuro, width = 12, height = 10, device = cairo_pdf)
# --- Plot 2: TH for iDA only (log2)
th_df <- differentiation_df %>%
dplyr::filter(Genes == "TH", neuron == "iDA") %>%
dplyr::mutate(log2_intensity = log2(mean_intensity), log2_sd_upper = log2(mean_intensity + sd_intensity) - log2_intensity, log2_sd_lower = log2_intensity - log2(pmax(mean_intensity - sd_intensity, 1)))
p_th <- ggplot(th_df, aes(x = genotype, y = log2_intensity, fill = genotype)) +
geom_bar(stat = "identity", width = 0.7) +
geom_errorbar(aes(ymin = log2_intensity - log2_sd_lower, ymax = log2_intensity + log2_sd_upper), width = 0.2) +
scale_fill_manual(values = c("ctrl" = "grey50", setNames(rep("grey80", length(all_genotypes) - 1), all_genotypes[all_genotypes != "ctrl"]))) +
labs(title = "TH Expression in iDA Neurons", x = "Genotype", y = expression(log[2](Mean~Intensity))) +
theme_bw(base_size = 6) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), panel.grid = element_blank(), legend.position = "none")
ggsave(file.path(out_dir_diffQC, "diff132_d50_TH_iDA.pdf"), p_th, width = 6, height = 4, device = cairo_pdf)
######
# --- Helper function: calculate annotation average with log2 transform and SEM
calc_annotation_avg_log2 <- function(genes, data) {
data %>%
dplyr::filter(Genes %in% genes) %>%
dplyr::group_by(genotype, neuron) %>%
dplyr::summarise(avg_intensity = mean(mean_intensity, na.rm = TRUE), sem_intensity = stats::sd(mean_intensity, na.rm = TRUE) / sqrt(dplyr::n()), .groups = "drop") %>%
dplyr::mutate(log2_intensity = log2(avg_intensity), log2_sem_upper = log2(avg_intensity + sem_intensity) - log2_intensity, log2_sem_lower = log2_intensity - log2(pmax(avg_intensity - sem_intensity, 1)))
}
# --- Helper function: plot annotation average (log2)
plot_annotation_avg_log2 <- function(data, title_text) {
ggplot(data, aes(x = genotype, y = log2_intensity, fill = neuron)) +
geom_bar(stat = "identity", position = position_dodge(width = 1), width = 0.7) +
geom_errorbar(aes(ymin = log2_intensity - log2_sem_lower, ymax = log2_intensity + log2_sem_upper), position = position_dodge(width = 1), width = 0.2) +
scale_fill_manual(values = c("iN" = "grey60", "iDA" = "grey80")) +
labs(title = title_text, x = "Genotype", y = expression(log[2](Mean~Intensity))) +
theme_bw(base_size = 6) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), panel.grid = element_blank())
}
# --- Plot 3: Presynaptic average (log2)
presyn_avg <- calc_annotation_avg_log2(presyn_genes, differentiation_df)
p_presyn <- plot_annotation_avg_log2(presyn_avg, "Presynaptic Markers (Average)")
ggsave(file.path(out_dir_diffQC, "diff132_d50_Presynaptic_avg.pdf"), p_presyn, width = 6, height = 4, device = cairo_pdf)
# --- Plot 4: Postsynaptic average (log2)
postsyn_avg <- calc_annotation_avg_log2(postsyn_genes, differentiation_df)
p_postsyn <- plot_annotation_avg_log2(postsyn_avg, "Postsynaptic Markers (Average)")
ggsave(file.path(out_dir_diffQC, "diff132_d50_Postsynaptic_avg.pdf"), p_postsyn, width = 6, height = 4, device = cairo_pdf)
# --- Plot 5: Neuronal markers average (log2)
neuro_avg <- calc_annotation_avg_log2(neuron_drivers, differentiation_df)
p_neuro_avg <- plot_annotation_avg_log2(neuro_avg, "Neuronal Markers (Average)")
ggsave(file.path(out_dir_diffQC, "diff132_d50_Neuronal_avg.pdf"), p_neuro_avg, width = 6, height = 4, device = cairo_pdf)
```
# vATPase by neuron x genotype
Plotted in ascending order
```{r}
# --------------------------------------------------- #
# neuro QC plots
# vATPase plots
# --------------------------------------------------- #
vatpase_genes <- c("ATP6V1G1", "ATP6AP2", "ATP6V1B2", "ATP6V1C1", "ATP6V1E1", "ATP6V1A", "ATP6V0D1",
"ATP6AP1", "ATP6V1F", "ATP6V0A1", "ATP6V0A2", "ATP6V1H", "ATP6V1D", "ROGDI",
"SYP", "DMXL1", "DMXL2", "PODXL", "PODXL2")
# --- 1. Filter and compute log2FC ---
vatpase_long_df <- clean_d50_transformed_df %>%
dplyr::filter(Genes %in% vatpase_genes, day == 50) %>%
dplyr::mutate(log2FC = log2(mean_ko / mean_ctrl)) %>%
dplyr::mutate(group = paste(genotype, neuron, sep = "_"))
# --- 2. Compute average log2FC per group ---
avg_fc_summary <- vatpase_long_df %>%
dplyr::group_by(group, neuron) %>%
dplyr::summarise(avg_log2FC = mean(log2FC, na.rm = TRUE), .groups = "drop")
# --- 3. Order by increasing avg log2FC ---
avg_fc_summary <- avg_fc_summary %>%
dplyr::arrange(neuron, avg_log2FC) %>%
dplyr::group_by(neuron) %>%
dplyr::mutate(group = factor(group, levels = unique(group)))
# --- 4. Merge ordering into main data frame ---
vatpase_long_df <- vatpase_long_df %>%
inner_join(avg_fc_summary %>% dplyr::select(group, neuron), by = c("group", "neuron")) %>%
dplyr::mutate(group = factor(group, levels = levels(avg_fc_summary$group)))
# --- 5. Plot ---
d50_vATPase_plot <- ggplot(vatpase_long_df, aes(x = group, y = log2FC, group = Genes)) +
geom_line(color = "grey70", size = 0.3, alpha = 0.8) +
geom_hline(yintercept =0, linetype = "dashed", color = "black") +
geom_point(color = "grey70", size = 0.5, alpha = 0.8) +
geom_line(data = avg_fc_summary,
aes(x = group, y = avg_log2FC, group = 1),
color = "red", size = 1.2, inherit.aes = FALSE) +
facet_wrap(~neuron, scales = "free_x") +
theme_bw(base_size = 6) +
theme(
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
panel.grid = element_blank()
) +
labs(
title = "v-ATPase Genes: log2FC Summary by Genotype × Neuron",
x = "Genotype × Neuron (sorted by mean log2FC)",
y = expression(log[2]*"(KO / Ctrl)")
)
d50_vATPase_plot
# --- 6. save pdf ---
ggsave(file.path(out_dir_d50_ascd_avg, "diff132_d50_Neuro_vATPase.pdf"), d50_vATPase_plot, width = 4, height = 3, device = cairo_pdf)
# --------------------------------------------------- #
# Function: vATPase color table by genotype (iN and iDA)
# --------------------------------------------------- #
make_vatpase_color_table <- function(
user_genotype,
input_df = clean_d50_transformed_df,
vatpase_gene_set = c(
# V1 sector
"ATP6V1A",
"ATP6V1B1","ATP6V1B2",
"ATP6V1C1","ATP6V1C2",
"ATP6V1D",
"ATP6V1E1","ATP6V1E2",
"ATP6V1F",
"ATP6V1G1","ATP6V1G2","ATP6V1G3",
"ATP6V1H",
# V0 sector (includes full c-ring and small subunits)
"ATP6V0A1","ATP6V0A2","ATP6V0A3","ATP6V0A4",
"ATP6V0C", "ATP6V0B",
"ATP6V0D1","ATP6V0D2",
"ATP6V0E1","ATP6V0E2",
"RNASEK",
# Accessory subunits
"ATP6AP1","ATP6AP2",
# Assembly or regulatory factors
"DMXL1","DMXL2","WDR7","ROGDI"
),
out_dir = out_dir_d50_ascd_avg,
pseudocount = 0
) {
stopifnot(dir.exists(out_dir))
# helper to compute log2FC robustly
safe_log2fc <- function(ko, ctrl, pc = 0) {
x <- log2((ko + pc) / (ctrl + pc))
x[!is.finite(x)] <- NA_real_
x
}
# 6WM2 chain→gene mapping (human subunits only; SidK chains X/Y/Z omitted)
chain_gene_map <- dplyr::bind_rows(
tibble::tibble(chain = "0", gene = "ATP6V0B"),
tibble::tibble(chain = as.character(1:9), gene = "ATP6V0C"),
tibble::tibble(chain = c("A","B","C"), gene = "ATP6V1A"),
tibble::tibble(chain = c("D","E","F"), gene = "ATP6V1B2"),
tibble::tibble(chain = "G", gene = "ATP6V1D"),
tibble::tibble(chain = c("H","I","J"), gene = "ATP6V1E1"),
tibble::tibble(chain = c("K","L","M"), gene = "ATP6V1G1"),
tibble::tibble(chain = "N", gene = "ATP6V1F"),
tibble::tibble(chain = "O", gene = "ATP6V1C1"),
tibble::tibble(chain = "P", gene = "ATP6V1H"),
tibble::tibble(chain = "Q", gene = "ATP6V0D1"),
tibble::tibble(chain = "R", gene = "ATP6V0A1"),
tibble::tibble(chain = "S", gene = "ATP6V0E1"),
tibble::tibble(chain = "T", gene = "RNASEK"),
tibble::tibble(chain = "U", gene = "ATP6V1S1"),
tibble::tibble(chain = "V", gene = "ATP6AP2")
)
# --- 1. Build long table with log2FC at day 50
global_long_df <- input_df %>%
dplyr::filter(Genes %in% vatpase_gene_set, day == 50) %>%
dplyr::transmute(
Genes,
genotype,
neuron,
log2FC = safe_log2fc(mean_ko, mean_ctrl, pseudocount)
)
# --- 2. Genotype-only symmetric domain across iN and iDA
geno_vals_vec <- global_long_df %>%
dplyr::filter(genotype == user_genotype, neuron %in% c("iN","iDA")) %>%
dplyr::pull(log2FC)
geno_vals_vec <- geno_vals_vec[is.finite(geno_vals_vec)]
if (length(geno_vals_vec) == 0) stop("No finite log2FC values for this genotype.")
geno_lim_val <- max(abs(range(geno_vals_vec, na.rm = TRUE)))
domain_vec <- c(-geno_lim_val, geno_lim_val)
# --- 3. Color mapper blue white red using genotype-only domain
col_fn_local <- scales::col_numeric(
palette = colorRampPalette(c("#2c7bb6", "white", "#d7191c"))(201),
domain = domain_vec,
na.color = "#B0B0B0"
)
# --- 3a. Save color scale CSV
colorscale_df <- tibble::tibble(
value = seq(domain_vec[1], domain_vec[2], length.out = 201),
color = toupper(col_fn_local(value))
)
colorscale_file <- file.path(out_dir, paste0("vATPase_colorscale_", user_genotype, ".csv"))
write.csv(colorscale_df, colorscale_file, row.names = FALSE)
# --- 3b. Save color scale PDF
cs_plot_df <- tibble::tibble(
x = seq(domain_vec[1], domain_vec[2], length.out = 500),
y = 1,
z = x
)
colorscale_pdf <- file.path(out_dir, paste0("vATPase_colorscale_", user_genotype, ".pdf"))
grDevices::cairo_pdf(colorscale_pdf, width = 4, height = 0.6)
print(
ggplot(cs_plot_df, aes(x = x, y = y, fill = z)) +
geom_raster() +
scale_fill_gradient2(
low = "#2c7bb6", mid = "white", high = "#d7191c",
midpoint = 0, limits = domain_vec, guide = "none"
) +
scale_x_continuous(breaks = c(domain_vec[1], 0, domain_vec[2])) +
labs(x = "log2FC", y = NULL, title = NULL) +
theme_minimal(base_size = 6) +
theme(
panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.margin = margin(4, 8, 4, 8)
)
)
grDevices::dev.off()
# --- 4. Build per-neuron tables and save CSVs
build_one_neuron_tbl <- function(neuron_label_val) {
df_neuron_df <- global_long_df %>%
dplyr::filter(genotype == user_genotype, neuron == neuron_label_val) %>%
dplyr::select(Genes, log2FC)
df_full_df <- tibble::tibble(Genes = vatpase_gene_set) %>%
dplyr::left_join(df_neuron_df, by = "Genes") %>%
dplyr::mutate(
Color = toupper(col_fn_local(log2FC))
) %>%
dplyr::rename(Gene = Genes)
out_with_chains_df <- chain_gene_map %>%
dplyr::left_join(df_full_df, by = c("gene" = "Gene")) %>%
dplyr::mutate(
chimerax_command = paste0("color /", chain, " ", Color)
) %>%
dplyr::select(Gene = gene, Chain = chain, Color, log2FC, chimerax_command)
chain_levels_vec <- c(as.character(0:9), LETTERS)
out_with_chains_df <- out_with_chains_df %>%
dplyr::mutate(Chain = factor(Chain, levels = chain_levels_vec)) %>%
dplyr::arrange(Chain)
out_file_path <- file.path(out_dir, paste0("vATPase_colors_", neuron_label_val, "_", user_genotype, "_6wm2.csv"))
write.csv(out_with_chains_df, out_file_path, row.names = FALSE)
out_with_chains_df
}
tbl_iN <- build_one_neuron_tbl("iN")
tbl_iDA <- build_one_neuron_tbl("iDA")
list(
iN = tbl_iN,
iDA = tbl_iDA,
scale_limits = domain_vec,
colorscale_file = colorscale_file,
colorscale_pdf = colorscale_pdf
)
}
# Call function for genotype of choice:
make_vatpase_color_table("MCOLN1")
```
# SV & SNARE plot by neuron x genotype
Plotted in ascending order
```{r}
# --------------------------------------------------- #
# neuro QC plots
# SV snare plots
# --------------------------------------------------- #
svsnare_genes <- c("SYP", "VAMP2", "STXBP1", "CPLX1", "SYT1", "SNAP25","SNAP","STX1A", "STX1")
# --- 1. Filter and compute log2FC
svsnare_long_df <- clean_d50_transformed_df %>%
dplyr::filter(Genes %in% svsnare_genes, day == 50) %>%
dplyr::mutate(log2FC = log2(mean_ko / mean_ctrl)) %>%
dplyr::mutate(group = paste(genotype, neuron, sep = "_"))
# --- 2. Compute average log2FC per group
avg_svsnare_fc_summary <- svsnare_long_df %>%
dplyr::group_by(group, neuron) %>%
dplyr::summarise(avg_log2FC = mean(log2FC, na.rm = TRUE), .groups = "drop")
# --- 3. Order by increasing avg log2FC
avg_svsnare_fc_summary <- avg_svsnare_fc_summary %>%
dplyr::arrange(neuron, avg_log2FC) %>%
dplyr::group_by(neuron) %>%
dplyr::mutate(group = factor(group, levels = unique(group)))
# --- 4. Merge ordering into main data frame
svsnare_long_df <- svsnare_long_df %>%
inner_join(avg_svsnare_fc_summary %>% dplyr::select(group, neuron), by = c("group", "neuron")) %>%
dplyr::mutate(group = factor(group, levels = levels(avg_svsnare_fc_summary$group)))
# --- 5. Plot
d50_svsnare_plot <- ggplot(svsnare_long_df, aes(x = group, y = log2FC, group = Genes)) +
geom_line(color = "grey70", size = 0.3, alpha = 0.8) +
geom_hline(yintercept =0, linetype = "dashed", color = "black") +
geom_point(color = "grey70", size = 0.5, alpha = 0.8) +
geom_line(data = avg_svsnare_fc_summary,
aes(x = group, y = avg_log2FC, group = 1),
color = "red", size = 1.2, inherit.aes = FALSE) +
facet_wrap(~neuron, scales = "free_x") +
theme_bw(base_size = 6) +
theme(
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
panel.grid = element_blank()
) +
labs(
title = "SV SNARE Genes: log2FC Summary by Genotype × Neuron",
x = "Genotype × Neuron (sorted by mean log2FC)",
y = expression(log[2]*"(KO / Ctrl)")
)
d50_svsnare_plot
# --- 6. save pdf ---
ggsave(file.path(out_dir_d50_ascd_avg, "diff132_d50_Neuro_SVSNSARE.pdf"), d50_svsnare_plot, width = 4, height = 3, device = cairo_pdf)
```
# SV docking plot by neuron x genotype
Plotted in ascending order
```{r}
# --------------------------------------------------- #
# neuro QC plots
# SV docking plots
# --------------------------------------------------- #
svdocking_genes <- c("SYP","VAMP2","SYB2","STXBP1","UNC18A","CPLX1","SYT1","SNAP25","SNAP","STX1A","STX1A","UNC13B")
# --- 1. Filter and compute log2FC
svdocking_long_df <- clean_d50_transformed_df %>%
dplyr::filter(Genes %in% svdocking_genes, day == 50) %>%
dplyr::mutate(log2FC = log2(mean_ko / mean_ctrl)) %>%
dplyr::mutate(group = paste(genotype, neuron, sep = "_"))
# --- 2. Compute mean and SD per group
svdocking_summary <- svdocking_long_df %>%
dplyr::group_by(group, neuron) %>%
dplyr::summarise(
avg_log2FC = mean(log2FC, na.rm = TRUE),
sd_log2FC = sd(log2FC, na.rm = TRUE),
.groups = "drop"
)
# --- 3. Order by increasing mean per neuron
svdocking_summary <- svdocking_summary %>%
dplyr::arrange(neuron, avg_log2FC) %>%
dplyr::group_by(neuron) %>%
dplyr::mutate(group_order = dplyr::row_number()) %>%
dplyr::ungroup()
# --- 4. Plot: SD in grey70, mean in red
d50_svdocking_plot <- ggplot(svdocking_summary,
aes(x = reorder(group, group_order), y = avg_log2FC, group = 1)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
# SD band per group
geom_linerange(aes(ymin = avg_log2FC - sd_log2FC,
ymax = avg_log2FC + sd_log2FC),
color = "grey70", linewidth = 1.2, alpha = 0.8) +
# mean line and points
geom_line(color = "red", linewidth = 0.9) +
geom_point(color = "red", size = 0.8) +
facet_wrap(~ neuron, scales = "free_x") +
theme_bw(base_size = 6) +
theme(
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
panel.grid = element_blank()
) +
labs(
title = "SV Docking Genes: mean ± SD by Genotype x Neuron",
x = "Genotype x Neuron (sorted by mean log2FC)",
y = expression(log[2]*"(KO / Ctrl)")
)
# --- 5. print and save
d50_svdocking_plot
ggsave(file.path(out_dir_d50_ascd_avg, "diff132_d50_Neuro_SVdocking.pdf"), d50_svdocking_plot, width = 4, height = 3, device = cairo_pdf)
```
# =================================================
# Module 4: Time Course Analysis (d30 vs d50)
# =================================================
# Abundance of neuromarkers in Ctrl iN/iDA
```{r}
# --------------------------------------------------- #
# neuro QC plots
# Abundance of neuromarkers in Ctrl day 30 & day 50
# --------------------------------------------------- #
# --- Step 1. Define Developmental Marker Gene List ----
dev_genes <- c("BDNF", "CAMK2B", "CRTC1", "DCX", "JUN", "MAP2", "NCAM1",
"NEFH", "NEFL", "NEFM", "NES", "POU3F2", "SLC17A7", "SYN1",
"SYP", "TUBB3", "TH", "BSN", "SYNJ1", "GAP43", "PSD95", "VGLUT",
"SYN2", "SYN3", "NGN")
# --- Step 2. Prepare Data ----
dev_plot_df <- clean_d50_transformed_df %>%
dplyr::filter(Genes %in% dev_genes) %>%
dplyr::select(Genes, neuron, day, mean_ctrl) %>%
dplyr::mutate(day = as.factor(day)) # Ensure day is treated as categorical for plotting
# --- Step 3.Plot ----
Dev_Marker_Abundance <- ggplot(dev_plot_df, aes(x = neuron, y = mean_ctrl, fill = day)) +
geom_bar(stat = "summary", fun = mean, position = position_dodge(width = 0.6), width = 0.6) +
facet_wrap(~Genes, scales = "free_y") +
scale_fill_manual(values = c("30" = "grey70", "50" = "steelblue")) +
labs(
title = "Neuronal Marker Abundance (Control, day 30 vs 50)",
x = "Neuron Type",
y = "Abundnace (mean_ctrl)",
fill = "Day"
) +
theme_bw() +
theme(
legend.position = "bottom",
panel.grid = element_blank()
)
# --- Step 4. Display ----
Dev_Marker_Abundance
# --- Step 5. Save ----
ggsave(file.path(out_dir_d50_timecourse,
"diff132_NeuroMarkers_ctrl_d30_d50.pdf"),
Dev_Marker_Abundance,
width = 12, height = 12, device = cairo_pdf
)
# --------------------------------------------------- #
# Neuro QC plots
# Abundance of neuromarkers (lineplot, normalized to day 30)
# --------------------------------------------------- #
# --- Step 1. Define Marker Gene List of Interest
marker_genes_line <- c("BSN", "NCAM1", "GAP43", "SYN1", "SYN3", "NEFL", "SYP")
# --- Step 2. Define function: plots Ctrl trajectories normalized to Day 30
plot_ctrl_markers_norm <- function(data, out_dir, genes = marker_genes_line) {
# Filter the dataset for Ctrl, iN/iDA, Day 30 & 50, and selected genes
ctrl_df <- data %>%
dplyr::filter(Genes %in% genes,
genotype == "ctrl",
neuron %in% c("iN", "iDA"),
day %in% c(30, 50)) %>%
dplyr::mutate(day = factor(day, levels = c(30, 50)))
# Normalize quan to Day 30 mean within each neuron × gene group
ctrl_norm <- ctrl_df %>%
dplyr::group_by(Genes, neuron) %>%
dplyr::mutate(norm_quan = quan / mean(quan[day == 30], na.rm = TRUE)) %>%
ungroup()
# Compute mean normalized intensity for each group
ctrl_means <- ctrl_norm %>%
dplyr::group_by(Genes, neuron, day) %>%
dplyr::summarise(mean_norm = mean(norm_quan, na.rm = TRUE), .groups = "drop")
# Get Day 50 values for end-of-line labels
labels_df <- ctrl_means %>% dplyr::filter(day == 50)
# Create the line plot
p <- ggplot(ctrl_means,aes(x = day, y = mean_norm, group = interaction(Genes, neuron), color = neuron)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
ggrepel::geom_text_repel(data = labels_df,
aes(label = Genes),
nudge_x = 0.25, hjust = 0,
segment.color = NA, size = 3,
show.legend = FALSE) +
facet_wrap(~neuron, ncol = 2) +
scale_color_manual(values = c("iN" = "grey70", "iDA" = "dodgerblue3")) +
scale_x_discrete(labels = c("30" = "Day 30", "50" = "Day 50")) +
labs(title = "Ctrl NeuroDev Markers Normalized to Day 30",
x = "Day", y = "Normalized Intensity (Day 30 = 1)") +
theme_bw() +
theme(panel.grid = element_blank(),
strip.text = element_text(face = "bold"),
legend.position = "none")
# Save as PDF
ggsave(file.path(out_dir, "ctrl_d30d50_Neuro_QC_Markers_norm_facet.pdf"),
p, width = 4, height = 4, device = cairo_pdf)
# save sourcedata
readr::write_csv(ctrl_means, file.path(out_dir, "ctrl_d30d50_Neuro_QC_Markers_norm_sourcedata.csv"))
return(p)
}
# --- Step 3. Call function and Save
Neuro_QC_Lineplot_norm <- plot_ctrl_markers_norm(cleaned_df, out_dir_d50_timecourse)
```
# Pre/Post Synaptic Markers in Ctrl iN/iDA
``` {r}
# --------------------------------------------------- #dev_plot_line_df
# neuro QC plots
# Abundance of Pre/Post Synaptic Markers in Ctrl day 30 & day 50
# --------------------------------------------------- #
# --- Step 1. Define Developmental Marker Gene List ----
presyn_genes <- c("BSN", "CPLX1", "CPLX2", "ELKS1", "MUNC13A", "MUNC13B", "MUNC18", "NRXN1",
"NRXN2", "NRXN3", "PCLO", "RAB3A", "RIM1", "RIM2", "RIMBP2", "SNAP25", "STX1A",
"STX1B", "SV2A", "SV2B", "SV2C", "SYN1", "SYN2", "SYP", "SYT1", "SYT2", "SYT7",
"VAMP1", "VAMP2", "VGLUT1", "VGLUT2", "VGLUT3"
)
postsyn_genes <- c(
"DLG4", "DLG3", "DLG2", "SHANK1", "SHANK2", "SHANK3", "HOMER1", "HOMER2", "HOMER3", "NLGN1",
"NLGN2", "NLGN3", "NLGN4X", "GPHN", "SYNGAP1", "GRIN1", "GRIN2A", "GRIN2B", "GRIA1", "GRIA2",
"GRIA3", "GRIA4", "GABRA1", "GABRA2", "GABRA3", "GABRA5", "GABRB1", "GABRB3", "GABRG2",
"CAMK2A", "CAMK2B", "PPP1R9B", "NRGN", "NPTX2"
)
# --- Step 2. Define presynaptic and postsynaptic gene lists ---
presyn_df <- tibble::tibble(Genes = presyn_genes, marker_type = "Presynaptic")
postsyn_df <- tibble::tibble(Genes = postsyn_genes, marker_type = "Postsynaptic")
synapse_genes_df <- bind_rows(presyn_df, postsyn_df)
# --- Step 3. Filter and transform data for control samples at Day 30/50 ---
synapse_long <- clean_d50_transformed_df %>%
dplyr::filter(Genes %in% synapse_genes_df$Genes) %>%
dplyr::select(Genes, neuron, day, mean_ctrl) %>%
dplyr::inner_join(synapse_genes_df, by = "Genes") %>%
dplyr::mutate(
day = paste("Day", day), # Convert numeric to "Day X"
log2_abundance = log2(mean_ctrl) # Apply log2 transform
)
# --- Step 4. Summarize per gene and reshape to wide for change calculation ---
delta_df <- synapse_long %>%
dplyr::group_by(Genes, neuron, marker_type, day) %>%
dplyr::summarise(log2_abundance = mean(log2_abundance, na.rm = TRUE), .groups = "drop") %>%
tidyr::pivot_wider(names_from = day, values_from = log2_abundance) %>%
dplyr::filter(!is.na(`Day 30`) & !is.na(`Day 50`)) %>%
dplyr::mutate(direction = ifelse(`Day 50` > `Day 30`, "up", "down"))
# --- Step 5. Reshape back to long format for plotting ---
long_with_direction <- delta_df %>%
tidyr::pivot_longer(cols = c(`Day 30`, `Day 50`), names_to = "day", values_to = "log2_abundance") %>%
dplyr::mutate(day = factor(day, levels = c("Day 30", "Day 50")))
# --- Step 6. Plot function with per-gene trajectories and group average overlays ---
SynapseTrajectoryPlot <- function(df) {
ggplot(df, aes(x = day, y = log2_abundance)) +
# Individual gene lines by direction
geom_line(aes(group = Genes, color = direction), alpha = 0.3, linewidth = 0.4) +
# Overlay group mean lines
stat_summary(
aes(group = interaction(marker_type, neuron)),
fun = mean, geom = "line",
color = "grey40", linewidth = 2
) +
stat_summary(
aes(group = interaction(marker_type, neuron)),
fun = mean, geom = "point",
color = "grey30", size = 4
) +
facet_grid(marker_type ~ neuron) +
scale_color_manual(values = c("up" = "#e31a1c", "down" = "#1f78b4")) +
theme_bw() +
labs(
title = "log2 Abundance Trajectories of \n Synaptic Genes (Control Neurons)",
x = "Day", y = "log2(Abundance)",
color = "Direction"
) +
theme(
panel.grid = element_blank(),
legend.position = "bottom",
)
}
# --- Step 7. Generate and save plot ---
Synaptic_Marker_Abundance <- SynapseTrajectoryPlot(long_with_direction)
Synaptic_Marker_Abundance
ggsave(file.path(out_dir_d50_timecourse,"diff132_SynapticMarkers_ctrl_d30_d50.pdf"), Synaptic_Marker_Abundance, width = 4, height = 5, device = cairo_pdf)
```
# NeuroDev protein-list trajectory d30 to d50
```{r}
# --------------------------------------------------- #
# neuro QC plots
# NeuroDev protein trajectory for day 30 & day 50 for select genotypes
# --------------------------------------------------- #
# For Ctrl, GRN, ASAH1, GBA1 and SMPD1 have both day30 and day50 data. plot neuronal markers etc and see if they increase.
# Define your gene sets
es_drivers <- c("ATF1", "CUX1", "MKI67", "NANOG", "POU5F1", "SOX2")
neuron_drivers <- c(
"BDNF", "CAMK2B", "CRTC1", "DCX", "JUN", "MAP2", "NCAM1",
"NEFH", "NEFL", "NEFM", "NES", "POU3F2", "SLC17A7", "SYN1",
"SYP", "TUBB3", "TH", "BSN", "SYNJ1", "GAP43", "PSD95", "VGLUT",
"SYN2", "SYN3", "NGN"
)
# --- Step 1. Create marker label dataframe
marker_genes_df <- tibble(
Genes = c(es_drivers, neuron_drivers),
marker_type = c(rep("ES-driver", length(es_drivers)),
rep("Neuron-driver", length(neuron_drivers)))
)
# --- Step 2. Filter for NeuroDev == TRUE
d30_neurodev <- complete_data_annotated_day30 %>% dplyr::filter(NeuroDev == TRUE)
d50_neurodev <- complete_data_annotated_day50 %>% dplyr::filter(NeuroDev == TRUE)
# --- Step 3. Reshape to long format and add 'day'
pivot_neurodev <- function(df, day_label) {
df %>%
dplyr::select(Genes, matches("_log2_fold_change$")) %>%
tidyr::pivot_longer(
cols = -Genes,
names_to = "sample",
values_to = "log2FC"
) %>%
dplyr::mutate(
day = day_label,
sample = str_remove(sample, "_log2_fold_change$"),
genotype = str_extract(sample, "^[^_]+"),
neuron = str_extract(sample, "(?<=_)(iDA|iN)$")
)
}
d30_long <- pivot_neurodev(d30_neurodev, "Day 30")
d50_long <- pivot_neurodev(d50_neurodev, "Day 50")
# --- Step 4. Combine
combined_df <- bind_rows(d30_long, d50_long)
# Join annotation info
combined_df_annotated <- combined_df %>%
dplyr::inner_join(marker_genes_df, by = "Genes")
# Now you can summarize or plot
summary_df <- combined_df_annotated %>%
dplyr::filter(genotype %in% c("ctrl", "ASAH1", "GBA1", "SMPD1", "GRN")) %>%
dplyr::group_by(marker_type, genotype, neuron, day) %>%
dplyr::summarise(mean_log2FC = mean(log2FC, na.rm = TRUE), .groups = "drop")
# --- Step 5. Plot
NeuoDevd30d50 <- ggplot(
combined_df_annotated %>% dplyr::filter(genotype %in% c("ctrl", "ASAH1", "GBA1", "SMPD1", "GRN")),
aes(x = day, y = log2FC, group = interaction(Genes, neuron), marker_type = marker_type)
) +
# Individual traces
geom_line(aes(color = neuron), alpha = 0.2, linewidth = 0.5) +
geom_point(aes(color = neuron), alpha = 0.4, size = 1) +
# Mean lines
stat_summary(
aes(group = neuron, color = neuron),
fun = mean, geom = "line", linewidth = 1.2
) +
stat_summary(
aes(group = neuron, fill = neuron),
fun = mean, geom = "point", size = 2.5, color = "black", shape = 21
) +
facet_grid(rows = vars(marker_type), cols = vars(genotype)) +
scale_color_manual(values = c("iN" = "grey60", "iDA" = "grey80")) +
scale_fill_manual(values = c("iN" = "grey60", "iDA" = "grey80")) +
labs(
title = "log2FC of NeuroDev Genes by Marker Type and Genotype",
x = "Day", y = "log2 Fold Change"
) +
theme_bw() +
theme(
strip.text = element_text(face = "bold"),
panel.grid = element_blank()
)
NeuoDevd30d50
# --- Step 6. Save
ggsave(file.path(out_dir_d50_timecourse, "diff132_NeuoDevd30d50.pdf"), NeuoDevd30d50, width = 8, height = 6, device = cairo_pdf)
```
# Summary plots for neuronal markers Ctrl d30 / d50
```{r}
# --- Summary barplots: NeuroDev marker lists in ctrl (d30 vs d50)
# dev_genes, neuron_drivers, presyn_genes, postsyn_genes
# two versions per list: raw quan and log2(quan)
# Step 1: combine all marker lists into one long df
marker_lists_nd <- list(
dev_genes = dev_genes,
neuron_drivers = neuron_drivers,
presyn_genes = presyn_genes,
postsyn_genes = postsyn_genes
)
summary_neurodev_df <- lapply(names(marker_lists_nd), function(lst) {
clean_d50_transformed_df %>%
dplyr::filter(Genes %in% marker_lists_nd[[lst]]) %>%
dplyr::select(Genes, neuron, day, mean_ctrl) %>%
dplyr::mutate(marker_list = lst)
}) %>% dplyr::bind_rows()
readr::write_csv(summary_neurodev_df, file.path(out_dir_d50_timecourse, "neurodev_markerlist_sourcedata.csv"))
# Step 2: normalize per protein per neuron to d30 = 1, then mean ± SEM across proteins
summary_neurodev_norm <- summary_neurodev_df %>%
dplyr::group_by(Genes, neuron, day, marker_list) %>%
dplyr::summarise(mean_ctrl = mean(mean_ctrl, na.rm = TRUE), .groups = "drop") %>%
tidyr::pivot_wider(names_from = day, values_from = mean_ctrl, names_prefix = "d") %>%
dplyr::filter(!is.na(d30), !is.na(d50)) %>%
dplyr::mutate(ratio_d50 = d50 / d30) %>%
tidyr::pivot_longer(cols = c(d30, ratio_d50),
names_to = "day_norm", values_to = "norm_val") %>%
dplyr::mutate(day_norm = dplyr::recode(day_norm, "d30" = "d30", "ratio_d50" = "d50"),
norm_val = dplyr::if_else(day_norm == "d30", 1, norm_val))
summary_neurodev_stats <- summary_neurodev_norm %>%
dplyr::group_by(marker_list, neuron, day_norm) %>%
dplyr::summarise(
mean_val = mean(norm_val, na.rm = TRUE),
sem_val = sd(norm_val, na.rm = TRUE) / sqrt(dplyr::n()),
n = dplyr::n(),
.groups = "drop"
) %>%
dplyr::mutate(group = factor(paste0(neuron, "_", day_norm),
levels = c("iN_d30", "iN_d50", "iDA_d30", "iDA_d50")))
# save source data csv
readr::write_csv(summary_neurodev_stats, file.path(out_dir_d50_timecourse, "neurodev_markerlist_stats_sourcedata.csv"))
# Step 3: barplot function
plot_nd_bar <- function(lst_name, stats_df, out_dir) {
dfp <- dplyr::filter(stats_df, marker_list == lst_name)
p <- ggplot(dfp, aes(x = day_norm, y = mean_val, fill = neuron)) +
geom_hline(yintercept = 1, linetype = "dashed", linewidth = 0.3, color = "black") +
geom_bar(stat = "identity", width = 1.0, color = "black", linewidth = 0.3) +
geom_errorbar(aes(ymin = mean_val - sem_val, ymax = mean_val + sem_val),
width = 0.3, linewidth = 0.3) +
scale_fill_manual(values = c("iN" = "grey50", "iDA" = "grey80")) +
facet_wrap(~ neuron, nrow = 1) +
labs(title = paste0(lst_name, " – ctrl, normalized to d30"),
x = NULL, y = "Relative abundance (d30 = 1)") +
theme_bw(base_size = 6) +
theme(panel.grid = element_blank(),
legend.position = "none",
axis.text.x = element_text(angle = 90, hjust = 1),
strip.background = element_rect(fill = "grey90", color = NA))
ggsave(file.path(out_dir, paste0("neurodev_

diff132_d50_nDIA.Rmd at commit 6587ebc, under MIT · at the source

Overview

Authors: Felix Kraus¹, Yuchen He¹, Yizhi Jiang^1,2, Delong Li^2,3, Yohannes A. Ambaw⁴, Federico M. Gasparoli¹, Joao A. Paulo¹, Tobias C. Walther^4,5, Robert V. Farese Jr.⁴, Steven P. Gygi¹, Florian Wilfling^2,3, J. Wade Harper^1,2

ORCID iDs: Felix Kraus, Yizhi Jiang, Delong Li, Yohannes A. Ambaw, Joao A. Paulo, Steven P. Gygi, Florian Wilfling, J. Wade Harper

Department of Cell Biology, Harvard Medical School, Boston, MA 02115
Aligning Science Across Parkinson’s Collaborative Research Network, Chevy Chase, MD 20815
Mechanisms of Cellular Quality Control, Max Planck Institute of Biophysics, Frankfurt 60438, Germany
Cell Biology Program, Sloan Kettering Institute, New York, NY 10065
Howard Hughes Medical Institute, New York, NY 10065

Institutions: Harvard University (United States); Max Planck Institute of Biophysics (Germany); Howard Hughes Medical Institute (United States)

Journal: Proceedings of the National Academy of Sciences of the United States of America, volume 123, issue 27, article e2609132123

Dates: received 20 March 2026; accepted 7 June 2026; published online 1 July 2026; in print 7 July 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1073/pnas.2609132123 · PMID 42384675 · PMCID PMC13342943 · OpenAlex W7166882155

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), cellular / molecular (subfield)

Keywords: lysosome, proteomics, iNeurons, organelle, protein interactions

MeSH: Dopaminergic Neurons*, Lysosomal Storage Diseases*, Proteome*, Animals, Cerebral Cortex, Glucosylceramidase, Human Embryonic Stem Cells, Humans, Lysosomes, Mitochondria, Proteomics (* major topic)

Journal subjects: Biological Sciences, Cell Biology

Topic: Lysosomal Storage Disorders Research (Physiology, Medicine), according to OpenAlex

Funding: NIH (NS110395, GM132129); Aligning Science Across Parkinson’s (024268)

Citations: not cited yet (Europe PMC); 51 references in the paper

Abstract

Lysosomes maintain cellular homeostasis by degrading proteins delivered via endocytosis and autophagy and by recycling building blocks for organelle biogenesis. Lysosomal storage disorders (LSDs) comprise a group of diseases affecting diverse lysosomal functions. To facilitate molecular phenotyping across diverse LSD gene classes, we are developing a library of human embryonic stem cells engineered to lack individual LSD genes as a resource for the field. Here, we report our initial stem cell toolkit lacking one of 23 LSD genes, including the majority of genes associated with sphingolipidoses and neuronal ceroid lipofuscinoses, and its use in the generation of a proteomic resource for induced cortical-like and midbrain dopaminergic-like neurons. In-depth abundance and correlation profiling across organelles and suborganelle components revealed potential vulnerabilities that reflect distinct patterns of proteome alterations across both genotypes and neuronal cell types. We characterize alterations in the mitochondrial proteome associated with GBA1 and ASAH1 deficiency and identify synaptic and mitochondrial defects in ASAH1−/− induced neurons that correlate with defects in neuronal firing rates. Moreover, we developed an informatic pipeline for proteome-wide identification of individual protein-protein interactions and protein complexes that may be disrupted as a result of LSD gene deficiency. Finally, we visualized structural alterations of ASAH1-deficient endolysosomes in situ using cryoelectron tomography, revealing swollen organelles that were largely devoid of dense internal membranes characteristic of wild-type cells, but containing numerous intralumenal vesicle compartments. This toolkit and associated proteomic landscapes provide a resource for defining molecular signatures associated with LSD gene dysfunction and organelle vulnerability.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 7 matches between paragraphs and lines of code.

sauerkrausi/neuroLSD

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 6587ebc16686d3c60f7338de524fc7344c1d929f, 1 July 2026

Languages: R (11), Python (9)

Size: 41 files, 20 scripts

Software Heritage: not archived

Found in: “Data, Materials, and Software Availability”

Holds: README, license file, environment (imaging/Calcium/requirements.txt, imaging/Calcium/requirements_MPS.txt, imaging/endolyso/requirements_MPS.txt, imaging/TH/requirements_MPS.txt), 10 notebooks

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: tidyverse (10 files), cowplot (8 files), patchwork (8 files), pheatmap (8 files), ggplot2 (7 files), NumPy (7 files), tifffile (7 files), ggpubr (6 files), Matplotlib (6 files), pandas (6 files), scikit-image (6 files), data.table (5 files), reshape2 (5 files), Cellpose (4 files), Plotly (3 files), PyTorch (3 files), broom (2 files), circlize (2 files), clusterProfiler (2 files), igraph (2 files), limma (2 files), rstatix (2 files), SciPy (2 files), ComplexHeatmap (1 file), seaborn (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

22 files

imaging/Calcium/cellposeSAM_CNER_Calcium-liveCell_template_v3.Rmd, R, 1,291 lines
imaging/Calcium/neuroLSD_CalciumMaster.Rmd, R, 1,386 lines
imaging/Calcium/run_cellpose_calcium_remoteHDD.py, Python, 112 lines
imaging/Calcium/run_cellpose_calcium_remoteHDD_MPS.py, Python, 118 lines
imaging/TH/diff118_iNiDA_THeval.Rmd, R, 165 lines
imaging/TH/th_neuron_quantification_interative.py, Python, 288 lines
imaging/endolyso/HeLa_lyso_eval.Rmd, R, 646 lines
imaging/endolyso/neuro_syn_env.py, Python, 953 lines
imaging/endolyso/neuron_endo_syn_eval.Rmd, R, 507 lines
imaging/endolyso/pyEEA1TFNeval.py, Python, 200 lines
imaging/endolyso/run_cellpose_EndoLyso_HeLa.py, Python, 343 lines
imaging/endolyso/split_channels_from_composite.py, Python, 48 lines
lipidome/Lipidomics_HeLa_iN-diff133_ASAH1-WC-OrganellIeIP.Rmd, R, 423 lines, 1 match
proteome/HeLa_Ctrl-ASAH1_LysoIP.Rmd, R, 1,637 lines
proteome/diff118_iNiDA_d23.Rmd, R, 1,314 lines
proteome/diff132_d50_nDIA.Rmd, R, 4,896 lines, 6 matches
proteome/diff136_iNd35_ctrl_asah1e1_axonalproteome.Rmd, R, 2,987 lines
utils/generate_readmes.py, Python, 165 lines
utils/get_r_versions.R, R, 80 lines
utils/update_readme_csv.py, Python, 155 lines
LICENSE, License, 21 lines
README.md, Text, 110 lines

At the source: github.com/sauerkrausi/neuroLSD

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
20 scripts, each with its path and the digest of its content;
7 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

ebi.ac.uk/emdb/emd-55210, at EMBL-EBI; found in “Data, Materials, and Software Availability”
ebi.ac.uk/emdb/emd-55211, at EMBL-EBI; found in “Data, Materials, and Software Availability”
zenodo:17296003, at Zenodo; found in “Data, Materials, and Software Availability”

Data, Materials, and Software Availability

Proteomic data (.RAW files) for nDIA of iN/iDA day 30, 50, and 70 iNeurons are deposited in the MASSIVE Database (https://massive.ucsd.edu/ProteoSAFe/static/massive.jsp): accession number MSV000099237 (https://massive.ucsd.edu/ProteoSAFe/dataset.jsp?task=cafb5f85380547faae2d32c7a9f7573f). HeLa LysoIP TMT data were deposited into ProteomeXchange (50): PXD067219 (https://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD067219). Lipidomics data are deposited on Metabolomics Workbench (51): Study ID ST004217; 10.21228/M8Z556 (https://doi.org/10.21228/M8Z556). Raw Cryo-ET tomograms have been deposited at EMDB under accession numbers EMD-55210 (https://www.ebi.ac.uk/emdb/EMD-55210) (Control) and EMD-55211 (https://www.ebi.ac.uk/emdb/EMD-55211) (ASAH1−/−). A Key Resource Table containing reagents and materials, source data, segmentation models, and datasets associated with this publication are available from Zenodo.org: 10.5281/zenodo.17296003 (https://doi.org/10.5281/zenodo.17296003). The Key Resource Table also contains unique Cellosaurus (https://www.cellosaurus.org/) identifiers for all cell lines reported here, which are available upon request with requisite Material Transfer Agreements from WiCell. Scripts affiliated with this work are available on Github under the repository https://github.com/sauerkrausi/neuroLSD. Data viewer can be found at https://wren.hms.harvard.edu/ProteomeLSDNeuron/.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 12 authors, 5 keywords, 11 MeSH terms, 2 funders, 51 references.

Cite

This paper

Kraus, F., He, Y., Jiang, Y., Li, D., Ambaw, Y. A., Gasparoli, F. M., Paulo, J. A., Walther, T. C., Farese, R. V., Gygi, S. P., Wilfling, F., & Harper, J. W. (2026). A human lysosomal storage disorder toolkit for decoding proteome landscapes in cortical-like and dopaminergic-like induced neurons. Proceedings of the National Academy of Sciences of the United States of America, 123(27), e2609132123. https://doi.org/10.1073/pnas.2609132123

BibTeX

@article{kraus2026human, author = {Kraus, Felix and He, Yuchen and Jiang, Yizhi and Li, Delong and Ambaw, Yohannes A. and Gasparoli, Federico M. and Paulo, Joao A. and Walther, Tobias C. and Farese, Robert V. and Gygi, Steven P. and Wilfling, Florian and Harper, J. Wade}, title = {{A human lysosomal storage disorder toolkit for decoding proteome landscapes in cortical-like and dopaminergic-like induced neurons}}, journal = {Proceedings of the National Academy of Sciences of the United States of America}, year = {2026}, month = jul, volume = {123}, number = {27}, pages = {e2609132123}, publisher = {National Academy of Sciences}, issn = {0027-8424}, doi = {10.1073/pnas.2609132123}, url = {https://doi.org/10.1073/pnas.2609132123}, pmid = {42384675}, pmcid = {PMC13342943} }

RIS

TY - JOUR AU - Kraus, Felix AU - He, Yuchen AU - Jiang, Yizhi AU - Li, Delong AU - Ambaw, Yohannes A. AU - Gasparoli, Federico M. AU - Paulo, Joao A. AU - Walther, Tobias C. AU - Farese, Robert V. AU - Gygi, Steven P. AU - Wilfling, Florian AU - Harper, J. Wade TI - A human lysosomal storage disorder toolkit for decoding proteome landscapes in cortical-like and dopaminergic-like induced neurons T2 - Proceedings of the National Academy of Sciences of the United States of America J2 - Proc Natl Acad Sci U S A PY - 2026 DA - 2026/07/01 VL - 123 IS - 27 SP - e2609132123 SN - 0027-8424 PB - National Academy of Sciences DO - 10.1073/pnas.2609132123 UR - https://doi.org/10.1073/pnas.2609132123 LA - en ER -

CSL-JSON

{ "id": "10.1073/pnas.2609132123", "type": "article-journal", "title": "A human lysosomal storage disorder toolkit for decoding proteome landscapes in cortical-like and dopaminergic-like induced neurons", "container-title": "Proceedings of the National Academy of Sciences of the United States of America", "author": [ { "family": "Kraus", "given": "Felix" }, { "family": "He", "given": "Yuchen" }, { "family": "Jiang", "given": "Yizhi" }, { "family": "Li", "given": "Delong" }, { "family": "Ambaw", "given": "Yohannes A." }, { "family": "Gasparoli", "given": "Federico M." }, { "family": "Paulo", "given": "Joao A." }, { "family": "Walther", "given": "Tobias C." }, { "family": "Farese", "given": "Robert V." }, { "family": "Gygi", "given": "Steven P." }, { "family": "Wilfling", "given": "Florian" }, { "family": "Harper", "given": "J. Wade" } ], "container-title-short": "Proc Natl Acad Sci U S A", "volume": "123", "issue": "27", "page": "e2609132123", "DOI": "10.1073/pnas.2609132123", "PMID": "42384675", "PMCID": "PMC13342943", "ISSN": "0027-8424", "publisher": "National Academy of Sciences", "URL": "https://doi.org/10.1073/pnas.2609132123", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 1 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.