Full text
Raw Data and Analysis Code for: Junction-Specific Cryptic Exon RT-qPCR Quantifies Sub-percent Nuclear TDP-43 Loss in Cells and ALS Cortex Shingo Koide Ichiko Ikegami Ryutaro Hanyu Yuka Mitsuhashi Koike Takuma Yamagishi Genri Toyama Aya Washida Mari Tada Akiyoshi Kakita Osamu Onodera Akihiro Sugai 5� 15, 2025 Contents 1 Introduction 2 1.1 RandPackageVersionsUsed.................................... 2 2 Figure 1: Workflow and primary data underlying selection of CE biomarkers 3 2.1 Figure 1C: RNA-seq Heatmap - Human Cell Lines . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Figure 1D: RNA-seq Heatmap - Human Brain Cell Types . . . . . . . . . . . . . . . . . . . . 3 2.3 Figure 1E and 1F: Baseline CE Levels in SH-SY5Y and HEK293T Cells . . . . . . . . . . . . 5 3 Figure 2: Junction-spanning primer redesign eliminates pre-mRNA co-amplification 7 3.1 Figure 2B: RT-qPCR Amplification Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.2 Figure 2C: Fold-change in CE Inclusion after TDP-43 Knockdown . . . . . . . . . . . . . . . 11 3.3 Figure 2D: Effect of Neuronal Differentiation on STMN2 Expression and Splicing . . . . . . . 14 4 Figure 3: Comparative sensitivity of CE and alternative-splicing markers 18 4.1 Figure 3C and 3D: Log2 Fold-changes of TDP-43-responsive RNA Events . . . . . . . . . . . 18 5 Figure 4: Junction-specific STMN2 CE assay detects TDP-43-deficient cells at <0.2% frequency 21 5.1 Figure 4B: Dilution-Limit Assay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 6 Figure 5: Junction-specific STMN2 assay reveals CE pathology in ALS motor cortex 27 6.1 Figure 5A: Relative Expression of TDP-43-Responsive RNA Events in ALS Motor Cortex (Junction-SpecificAssays)...................................... 27 6.2 Figure 5B: Relative Expression in ALS Motor Cortex (Conventional Assays) . . . . . . . . . . 30 1
7 Supplementary Figure S1: GAPDH/RPLP1 Shifts After TDP-43 Depletion 32 7.1 Supplementary Figure S1B and S1C: GAPDH mRNA Abundance Relative to RPLP1 . . . . 32 1 Introduction This document provides the R code used to generate the figures and perform statistical analyses presented in the manuscript “Junction-Specific Cryptic Exon RT-qPCR Quantifies Sub-percent Nuclear TDP-43 Loss in Cells and ALS Cortex”. This report is intended to accompany the raw data deposited in Zenodo, providing transparency and reproducibility of our findings. The R code is organized into sections corresponding to the figures in the main manuscript. Each section includes the necessary code to load the raw data (expected to be in a subdirectory named Rawdata/), process it, perform statistical tests, and generate the figures. To run this code, ensure that R and all required packages (listed below) are installed, and that the raw data CSV files are located in the Rawdata/ subdirectory relative to this R Markdown file. Data Availability: The raw data files (CSV format) referenced in this document are available from the Zenodo repository. Required Packages and Session Information The following R packages were used for data analysis and visualization. library(tidyverse) library(pheatmap) library(grid) library(ggbeeswarm) library(ggsignif) library(broom) library(emmeans) library(effsize) library(knitr) 1.1 R and Package Versions Used This analysis was performed using: R (base): R version 4.4.1 (2024-06-14 ucrt) knitr: 1.47 effsize: 0.8.1 emmeans: 1.10.2 broom: 1.0.6 ggsignif: 0.6.4 ggbeeswarm: 0.7.2 pheatmap: 1.0.12 lubridate: 1.9.3 forcats: 1.0.0 stringr: 1.5.1 dplyr: 1.1.4 purrr: 1.0.2 readr: 2.1.5 tidyr: 1.3.1 tibble: 3.2.1 ggplot2: 3.5.1 tidyverse: 2.0.0 2
2 Figure 1: Workflow and primary data underlying selection of CE biomarkers 2.1 Figure 1C: RNA-seq Heatmap - Human Cell Lines This section generates a heatmap showing the expression (nTPM) of candidate CE-harboring genes across various human cell lines, based on public RNA-seq data from The Human Protein Atlas. # Define markers and cell lines markers_fig1c <- c("STMN2","UNC13A","SYT7","CAMK2B","ACTL6B", "ARHGAP32","CDK7","HDGFL2","ATG4B","GPSM2") cell_lines_fig1c <- c("SH-SY5Y","SK-N-SH","HEK293","A-549","HepG2", "HeLa","U2OS","HL-60","MCF-7","HT-29") # Load and process RNA-seq data for cell lines rnaseq_cell_lines <- read_tsv("Rawdata/rna_celline.tsv",show_col_types = FALSE)%>% filter(`Gene name`%in% markers_fig1c, `Cell line`%in% cell_lines_fig1c) %>% select(`Gene name`,`Cell line`, nTPM) %>% pivot_wider(names_from = `Cell line`,values_from = nTPM) %>% column_to_rownames("Gene name") # Create log2-transformed matrix mat_cell_lines <- log2(as.matrix(rnaseq_cell_lines) +0.1) # Generate heatmap pheatmap(mat_cell_lines, cluster_rows = TRUE,cluster_cols = TRUE, fontsize_row = 10,fontsize_col = 10, main = "Gene Expression in Human Cell Lines (Fig. 1C)") grid.text("log2(nTPM + 0.1)",x = unit(0.95,"npc"), y = unit(0.9,"npc"), gp = gpar(fontsize = 8)) 2.2 Figure 1D: RNA-seq Heatmap - Human Brain Cell Types This section generates a heatmap showing the expression (nTPM) of candidate CE-harboring genes across major human brain cell types, based on single-cell RNA-seq data from The Human Protein Atlas. # Define cell types for brain scRNA-seq cell_types_fig1d <- c("Excitatory neurons","Inhibitory neurons","Astrocytes", "Oligodendrocyte precursor cells","Oligodendrocytes", "Microglial cells","Schwann cells", "Endothelial cells","Lymphatic endothelial cells", "Smooth muscle cells") # Load and process scRNA-seq data for brain cell types sc_brain_cells <- read_tsv("Rawdata/rna_single_cell_type.tsv",show_col_types = FALSE)%>% filter(`Gene name`%in% markers_fig1c, `Cell type`%in% cell_types_fig1d) %>% select(`Gene name`,`Cell type`, nTPM) %>% pivot_wider(names_from = `Cell type`,values_from = nTPM) %>% column_to_rownames("Gene name") 3
Gene Expression in Human Cell Lines (Fig. 1C) SH−SY5Y SK−N−SH HeLa HL−60 HT−29 MCF−7 HEK293 A−549 U2OS HDGFL2 CDK7 GPSM2 ATG4B CAMK2B SYT7 ARHGAP32 STMN2 ACTL6B UNC13A −2 0 2 4 6 8 log2(nTPM + 0.1) Figure 1: Figure 1C. Expression of candidate genes in human cell lines. 4
# Create log2-transformed matrix mat_brain_cells <- log2(as.matrix(sc_brain_cells) +0.1) # Generate heatmap pheatmap(mat_brain_cells, cluster_rows = TRUE,cluster_cols = TRUE, fontsize_row = 10,fontsize_col = 8, main = "Gene Expression in Brain Cell Types (Fig. 1D)") grid.text("log2(nTPM + 0.1)",x = unit(0.95,"npc"), y = unit(0.95,"npc"), gp = gpar(fontsize = 8)) 2.3 Figure 1E and 1F: Baseline CE Levels in SH-SY5Y and HEK293T Cells This section visualizes the baseline cryptic exon (CE) levels (normalized to RPLP1) in untreated differentiated SH-SY5Y neurons and HEK293T cells. # Load and process CE expression data df_fig1ef <- read.csv("Rawdata/1E-F.csv",stringsAsFactors = FALSE)%>% pivot_longer(cols = -type, names_to = "gene",values_to = "expression")%>% mutate( gene = recode(gene, "STMN2_old" ="STMN2", "ARHGAP32_old" ="ARHGAP32", "ACTL6B_old" ="ACTL6B") )%>% filter(gene %in% c("STMN2","ARHGAP32","ACTL6B", "UNC13A","CAMK2B","CDK7", "HDGFL2","SYT7","GPSM2","ATG4B")) # Determine detection limit and value for undetected points detect_min_fig1ef <- min(df_fig1ef$expression[df_fig1ef$expression >0], na.rm = TRUE) undetected_val_fig1ef <- detect_min_fig1ef /10 df_fig1ef <- df_fig1ef %>% mutate( expr_plot = if_else(expression >0, expression, NA_real_), is_undetected = (expression == 0) ) # Define gene order for plotting gene_levels_fig1ef <- c("STMN2","ARHGAP32","ACTL6B", "HDGFL2","UNC13A","CDK7", "ATG4B","GPSM2","SYT7","CAMK2B") # Function to draw the plot draw_plot_fig1ef <- function(data, cell_type_title) { ggplot(data, aes(x = factor(gene, levels = gene_levels_fig1ef), y = expr_plot)) + geom_beeswarm(size = 2,alpha = 0.7,color = "gray50")+ geom_jitter( data = data %>% filter(is_undetected), aes(x = factor(gene, levels = gene_levels_fig1ef), y = undetected_val_fig1ef), size = 2,width = 0.1,alpha = 0.7,height = 0,fill = "black",color = "gray50" 5
Gene Expression in Brain Cell Types (Fig. 1D) Schwann cells Endothelial cells Smooth muscle cells Excitatory neurons Inhibitory neurons Microglial cells Oligodendrocytes Astrocytes Oligodendrocyte precursor cells ACTL6B SYT7 UNC13A GPSM2 ATG4B CDK7 HDGFL2 ARHGAP32 CAMK2B STMN2 −2 0 2 4 6 8 log2(nTPM + 0.1) Figure 2: Figure 1D. Expression of candidate genes in human brain cell types. 6
)+ geom_hline(yintercept = undetected_val_fig1ef, linetype = "dotted",size = 0.5)+ scale_y_log10( limits = c(undetected_val_fig1ef/3,1e-2), breaks = c(undetected_val_fig1ef, 1e-7,1e-6,1e-5,1e-4,1e-3,1e-2), labels = c("Undetected",expression(10^-7), expression(10^-6), expression(10^-5), expression(10^-4), expression(10^-3), expression(10^-2)), expand = c(0,0) )+ labs( title = cell_type_title, x = NULL, y = expression(Delta*Ct~CE /RPLP1) )+ theme_bw(base_size = 10)+ theme( axis.text.y = element_text(size = rel(1)), plot.title = element_text(size = rel(1.2), hjust = 0.5,face = "bold"), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.line = element_line(color = "black",linewidth = 0.5), axis.text.x = element_text(size = rel(1), angle = 45,hjust = 1,face = "italic"), axis.title.y = element_text(size = rel(1.1)) ) } # Generate plots for SH-SY5Y and HEK293T plot_shsy5y <- draw_plot_fig1ef(df_fig1ef %>% filter(type == "SH-SY5Y"), "SH-SY5Y (Fig. 1E)") plot_hek293t <- draw_plot_fig1ef(df_fig1ef %>% filter(type == "HEK293T"), "HEK293T (Fig. 1F)") # Display plots (will be shown sequentially in PDF) print(plot_shsy5y) print(plot_hek293t) 3 Figure 2: Junction-spanning primer redesign eliminates premRNA co-amplification 3.1 Figure 2B: RT-qPCR Amplification Curves This section shows RT-qPCR amplification curves comparing pre-mRNA, conventional CE, and junctionspanning CE assays in differentiated SH-SY5Y cells under basal conditions. # Load and process data for amplification curves data_fig2b <- read.csv("Rawdata/2B.csv",stringsAsFactors = FALSE) num_cols_fig2b <- grepl("^[0-9]+$",colnames(data_fig2b)) colnames(data_fig2b)[num_cols_fig2b] <- paste0("Cycle",colnames(data_fig2b)[num_cols_fig2b]) 7
Undetected 10−7 10−6 10−5 10−4 10−3 10−2 STMN2 ARHGAP32 ACTL6B HDGFL2 UNC13A CDK7 ATG4B GPSM2 SYT7 CAMK2B ∆Ct CE RPLP1 SH−SY5Y (Fig. 1E) Figure 3: Figure 1E & 1F. Baseline CE relative expression (ΔCt to RPLP1). Plots for SH-SY5Y (top/left) and HEK293T (bottom/right) are generated sequentially. 8
Undetected 10−7 10−6 10−5 10−4 10−3 10−2 STMN2 ARHGAP32 ACTL6B HDGFL2 UNC13A CDK7 ATG4B GPSM2 SYT7 CAMK2B ∆Ct CE RPLP1 HEK293T (Fig. 1F) Figure 4: Figure 1E & 1F. Baseline CE relative expression (ΔCt to RPLP1). Plots for SH-SY5Y (top/left) and HEK293T (bottom/right) are generated sequentially. 9
print(plot_fig2d) STMN2 CE (Conventional Assay) STMN2 CE (Junction−Specific Assay) STMN2 (total mRNA) STMN2 (pre−mRNA) Undifferentiated Differentiated Undifferentiated Differentiated Undifferentiated Differentiated Undifferentiated Differentiated −2 −1 0 1 −2 −1 0 1 Differentiation Status Log2 Relative Expression (∆∆Ct, RPLP1) Figure 7: Figure 2D. Effect of neuronal differentiation on STMN2 isoforms. Plots are faceted by STMN2 isoform/assay type. # Calculate fold change statistics (Differentiated vs Undifferentiated) undiff_means_fig2d_table <- df_fig2d %>% filter(type == "Undifferentiated")%>% summarise(across(all_of(plot_vars_fig2d), mean, .names = "mean_{.col}"), .groups = 'drop') fold_stats_fig2d_table <- df_fig2d %>% filter(type == "Differentiated")%>% summarise( across( all_of(plot_vars_fig2d), list( Fold_Change_Mean = ~mean(.x /undiff_means_fig2d_table[[paste0("mean_",cur_column())]]), Fold_Change_SD = ~sd(.x /undiff_means_fig2d_table[[paste0("mean_",cur_column())]]) ), .names = "{.col}_{.fn}" ) )%>% 16
pivot_longer(everything(), names_to = c("Assay",".value"), names_pattern = "(.+)_(Fold_Change_Mean|Fold_Change_SD)")%>% mutate(Assay = recode(Assay, "STMN2" ="STMN2 (total mRNA)", "STMN2_premRNA" ="STMN2 (pre-mRNA)", "CE_conventional" ="STMN2 CE (Conventional Assay)", "CE_Refined" ="STMN2 CE (Junction-Specific Assay)")) cat("\nFold change statistics (Differentiated vs Undifferentiated):\n") ## ## Fold change statistics (Differentiated vs Undifferentiated): kable(fold_stats_fig2d_table, caption = "Fold change statistics for Figure 2D.",digits = 3, col.names = c("Assay / Isoform","Mean Fold Change","SD Fold Change")) Table 3: Fold change statistics for Figure 2D. Assay / Isoform Mean Fold Change SD Fold Change STMN2 (total mRNA) 1.571 0.037 STMN2 (pre-mRNA) 1.634 0.304 STMN2 CE (Conventional Assay) 0.806 0.071 STMN2 CE (Junction-Specific Assay) 0.212 0.007 # Perform t-tests for each gene_assay_label and report t_test_results_fig2d_table <- df_long_fig2d %>% group_by(gene_assay_label) %>% do(tidy(t.test(log2_expr_val ~type, data = .))) %>% select(Assay_Isoform = gene_assay_label, Mean_Diff_log2 = estimate, Statistic = statistic, P_Value = p.value, DF = parameter, CI_Low = conf.low, CI_High = conf.high) cat("\nT-test results (Differentiated vs Undifferentiated, on log2 expression):\n") ## ## T-test results (Differentiated vs Undifferentiated, on log2 expression): kable(t_test_results_fig2d_table, caption = "T-test results for Figure 2D.",digits = 3) Table 4: T-test results for Figure 2D. Assay_Isoform Mean_Diff_log2 Statistic P_Value DF CI_Low CI_High STMN2 (total mRNA) -0.652 -24.498 0.000 4.597 -0.722 -0.582 STMN2 (pre-mRNA) -0.701 -3.990 0.021 3.507 -1.217 -0.185 STMN2 CE (Conventional Assay) 0.312 3.562 0.026 3.812 0.064 0.560 STMN2 CE (Junction-Specific Assay) 2.239 39.413 0.000 4.491 2.088 2.390 17
4 Figure 3: Comparative sensitivity of CE and alternative-splicing markers 4.1 Figure 3C and 3D: Log2 Fold-changes of TDP-43-responsive RNA Events This section compares the log2 fold-changes of six TDP-43-responsive RNA events (three CEs and three alternative splicing events) after TDP-43 depletion in differentiated SH-SY5Y cells and HEK293T cells. Aberrant isoforms are normalized to their cognate canonical isoforms. # Load and process data for comparative sensitivity data_fig3cd <- read.csv("Rawdata/3C-D.csv",header = TRUE, check.names = FALSE,stringsAsFactors = FALSE) measurements_fig3cd <- setdiff(names(data_fig3cd), c("type","cell")) data_long_fig3cd <- data_fig3cd %>% mutate( cell = factor(cell, levels = c("SH-SY5Y","HEK293T")), type = factor(type, levels = c("Scramble","shTDP43","siCtrl","siTDP43")), condition = case_when( type %in% c("Scramble","siCtrl")~"Control", type %in% c("shTDP43","siTDP43")~"TDP-43 KD", TRUE ~as.character(type) ) )%>% pivot_longer(all_of(measurements_fig3cd), names_to = "measurement",values_to = "value")%>% mutate( measurement = factor(measurement, levels = measurements_fig3cd), log2_value = log2(value), condition = factor(condition, levels = c("Control","TDP-43 KD")) ) # Define colors for plot aesthetics group_colors_fig3cd <- c( "Control" ="gray70", "TDP-43 KD" ="#aa0000cc" ) # Function to create plot for each cell type make_plot_fig3cd <- function(df, cell_name_title) { plot_data_local <- df %>% filter(cell == cell_name_title) anova_model_fig3cd <- aov(log2_value ~measurement *condition, data = plot_data_local) cat(paste0("\nTwo-way ANOVA for ", cell_name_title, " (Log2_Fold_Change ~ RNA_Event * TDP43_Status):\n")) print(summary(anova_model_fig3cd)) anova_summary_df <- tidy(anova_model_fig3cd) interaction_p_val <- filter(anova_summary_df, term == "measurement:condition")$p.value if (!is.na(interaction_p_val) && interaction_p_val <0.05) { cat(paste0("\nPost-hoc (emmeans) for significant interaction in ", cell_name_title, " (comparing Control vs TDP-43 KD for each RNA event):\n")) 18
emm_fig3cd_interaction <- emmeans(anova_model_fig3cd, ~condition |measurement) contrast_results_interaction <- pairs(emm_fig3cd_interaction, adjust = "bonferroni") print(as.data.frame(contrast_results_interaction)) } ggplot(plot_data_local, aes(x = measurement, y = log2_value, fill = condition)) + stat_summary(fun = mean, geom = "crossbar",position = position_dodge(0.7), linewidth = 0.3,width = 0.9,color = "black")+ stat_summary(fun.data = mean_se, geom = "errorbar",position = position_dodge(0.7), linewidth = 0.3,width = 0.5,color = "black")+ geom_jitter(position = position_jitterdodge(jitter.width = 0.1,dodge.width = 0.7), shape = 21,size = 2.5,alpha = 0.6,color = "black")+ scale_y_continuous(limits = c(-7,12), expand = c(0,0), breaks = seq(-6,12,by=2)) + scale_fill_manual(values = group_colors_fig3cd) + labs(title = cell_name_title, x = "RNA Event", y = expression(Log[2]*" (Fold Change)"), fill = "TDP-43 Status")+ theme_bw(base_size = 10)+ theme( axis.text.x = element_text(size = rel(0.9), angle = 45,hjust = 1,color = "black"), axis.text.y = element_text(size = rel(1), color = "black"), plot.title = element_text(size = rel(1.2), face = "bold",hjust = 0.5), axis.title.x = element_text(size = rel(1.1)), axis.title.y = element_text(size = rel(1.1)), legend.position = "top", legend.text = element_text(size = rel(0.9)), legend.title = element_text(size = rel(1)), panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black",linewidth = 0.5) ) } # Generate plots for SH-SY5Y and HEK293T plot_shsy5y_fig3 <- make_plot_fig3cd(data_long_fig3cd, "SH-SY5Y") ## ## Two-way ANOVA for SH-SY5Y (Log2_Fold_Change ~ RNA_Event * TDP43_Status): ## Df Sum Sq Mean Sq F value Pr(>F) ## measurement 5 102.68 20.54 442 <2e-16 *** ## condition 1 266.51 266.51 5736 <2e-16 *** ## measurement:condition 5 102.68 20.54 442 <2e-16 *** ## Residuals 24 1.12 0.05 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Post-hoc (emmeans) for significant interaction in SH-SY5Y (comparing Control vs TDP-43 KD for each RNA event): ## measurement = STMN2_CE_refined: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -10.438333 0.1759985 24 -59.309 <.0001 ## ## measurement = ARHGAP32_CE: ## contrast estimate SE df t.ratio p.value 19
## Control - (TDP-43 KD) -5.573333 0.1759985 24 -31.667 <.0001 ## ## measurement = ACTL6B_CE: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -8.436667 0.1759985 24 -47.936 <.0001 ## ## measurement = POLDIP3: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -5.583333 0.1759985 24 -31.724 <.0001 ## ## measurement = STAG2: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -0.905000 0.1759985 24 -5.142 <.0001 ## ## measurement = TNIK: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -1.713333 0.1759985 24 -9.735 <.0001 plot_hek293t_fig3 <- make_plot_fig3cd(data_long_fig3cd, "HEK293T") ## ## Two-way ANOVA for HEK293T (Log2_Fold_Change ~ RNA_Event * TDP43_Status): ## Df Sum Sq Mean Sq F value Pr(>F) ## measurement 5 218.11 43.62 247.9 <2e-16 *** ## condition 1 84.99 84.99 483.0 <2e-16 *** ## measurement:condition 5 218.11 43.62 247.9 <2e-16 *** ## Residuals 36 6.33 0.18 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Post-hoc (emmeans) for significant interaction in HEK293T (comparing Control vs TDP-43 KD for each RNA event): ## measurement = STMN2_CE_refined: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -7.15500 0.2966193 36 -24.122 <.0001 ## ## measurement = ARHGAP32_CE: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -3.85000 0.2966193 36 -12.980 <.0001 ## ## measurement = ACTL6B_CE: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) 5.54875 0.2966193 36 18.707 <.0001 ## ## measurement = POLDIP3: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -6.92000 0.2966193 36 -23.330 <.0001 ## ## measurement = STAG2: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -2.24875 0.2966193 36 -7.581 <.0001 ## ## measurement = TNIK: ## contrast estimate SE df t.ratio p.value ## Control - (TDP-43 KD) -1.34250 0.2966193 36 -4.526 0.0001 20
print(plot_shsy5y_fig3) −6 −4 −2 0 2 4 6 8 10 12 STMN2_CE_refined ARHGAP32_CE ACTL6B_CE POLDIP3 STAG2 TNIK RNA Event Log2 (Fold Change) TDP−43 Status Control TDP−43 KD SH−SY5Y Figure 8: Figure 3C & 3D. Log2 fold-changes of TDP-43 responsive RNA events in SH-SY5Y (top/left) and HEK293T (bottom/right) cells after TDP-43 knockdown. Plots are generated sequentially. print(plot_hek293t_fig3) 5 Figure 4: Junction-specific STMN2 CE assay detects TDP-43deficient cells at <0.2% frequency 5.1 Figure 4B: Dilution-Limit Assay This section presents the results of a dilution-limit assay where cDNA from TDP-43-knockdown SH-SY5Y neurons was spiked into control cDNA. It compares the junction-specific STMN2 CE assay and POLDIP3 exon-skipping. 21
−6 −4 −2 0 2 4 6 8 10 12 STMN2_CE_refined ARHGAP32_CE ACTL6B_CE POLDIP3 STAG2 TNIK RNA Event Log2 (Fold Change) TDP−43 Status Control TDP−43 KD HEK293T Figure 9: Figure 3C & 3D. Log2 fold-changes of TDP-43 responsive RNA events in SH-SY5Y (top/left) and HEK293T (bottom/right) cells after TDP-43 knockdown. Plots are generated sequentially. 22
# Load and process data for dilution-limit assay df_fig4b <- read.csv("Rawdata/4B.csv") df_long_fig4b <- df_fig4b %>% pivot_longer( cols = c(STMN2_best, POLDIP3), names_to = "primer_assay", values_to = "fold_change_value" )%>% mutate( kd_fraction_label = gsub("KD_","", type), log2_fold_change = log2(fold_change_value) ) df_long_fig4b$kd_fraction_label <- factor(df_long_fig4b$kd_fraction_label, levels = c("0%","0.1%","1%","10%")) df_long_fig4b$primer_assay <- factor(df_long_fig4b$primer_assay, levels = c("STMN2_best","POLDIP3"), labels = c("STMN2 (Junction-Specific)", "POLDIP3 (Exon Skip)")) # Perform two-way ANOVA anova_model_fig4b <- aov(log2_fold_change ~primer_assay *kd_fraction_label, data = df_long_fig4b) anova_summary_fig4b <- summary(anova_model_fig4b) cat("Two-way ANOVA summary (Log2_Fold_Change ~ Assay_Type * KD_Fraction):\n") ## Two-way ANOVA summary (Log2_Fold_Change ~ Assay_Type * KD_Fraction): print(anova_summary_fig4b) ## Df Sum Sq Mean Sq F value Pr(>F) ## primer_assay 1 15.89 15.893 504.1 1.60e-13 *** ## kd_fraction_label 3 34.71 11.571 367.0 5.87e-15 *** ## primer_assay:kd_fraction_label 3 16.60 5.532 175.5 1.89e-12 *** ## Residuals 16 0.50 0.032 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova_tidy_fig4b <- tidy(anova_model_fig4b) interaction_p_fig4b <- filter(anova_tidy_fig4b, term == "primer_assay:kd_fraction_label")$p.value interaction_label_fig4b <- if (is.na(interaction_p_fig4b)) "Interaction N/A" else if (interaction_p_fig4b <0.001) { "Assay x Dilution Interaction: p < 0.001" }else { paste0("Assay x Dilution Interaction: p = ",sprintf("%.3f", interaction_p_fig4b)) } cat(paste0("\n", interaction_label_fig4b, "\n")) ## ## Assay x Dilution Interaction: p < 0.001 23
# Generate plot plot_fig4b <- ggplot(df_long_fig4b, aes(x = kd_fraction_label, y = log2_fold_change, group = primer_assay, color = primer_assay, fill = primer_assay)) + stat_summary(fun = mean, geom = "line",linewidth = 1)+ stat_summary(fun = mean, geom = "crossbar",linewidth = 0.3, width = 0.6,show.legend = FALSE)+ stat_summary(fun.data = mean_se, geom = "errorbar",width = 0.3,show.legend = FALSE)+ geom_point(size = 2,shape = 21,alpha = 0.8,show.legend = FALSE)+ scale_y_continuous(limits = c(-1,6), expand = c(0,0), breaks = seq(-1,6,1)) + scale_color_manual(values = c( "STMN2 (Junction-Specific)" ="#008080ff", "POLDIP3 (Exon Skip)" ="gray40" )) + scale_fill_manual(values = c( "STMN2 (Junction-Specific)" ="#008080ff", "POLDIP3 (Exon Skip)" ="gray40" )) + labs( x = "Knockdown cDNA Fraction", y = expression(Log[2]*" (Fold Change vs. 0% KD)"), color = "Assay Target", fill = "Assay Target" )+ theme_bw(base_size = 10)+ theme( axis.text.x = element_text(size = rel(1), color ="black"), axis.text.y = element_text(size = rel(1), color ="black"), axis.title = element_text(size = rel(1.1)), plot.title = element_text(size = rel(1.2), face = "bold",hjust = 0.5), legend.position = "top", legend.text = element_text(size = rel(0.9)), legend.title = element_text(size = rel(1)), panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black",linewidth = 0.5) )+ annotate( "text", x = Inf, y = Inf, label = interaction_label_fig4b, size = 3, hjust = 1.05, vjust = 1.5 ) print(plot_fig4b) 24
Assay x Dilution Interaction: p < 0.001 −1 0 1 2 3 4 5 6 0% 0.1% 1% 10% Knockdown cDNA Fraction Log2 (Fold Change vs. 0% KD) Assay Target STMN2 (Junction−Specific) POLDIP3 (Exon Skip) Figure 10: Figure 4B. Dilution-limit assay comparing STMN2 (Junction-Specific) and POLDIP3 (Exon Skip) assays. # Calculate linearity (slope and R-squared) for raw fold change values df_fig4b_linear_calc <- df_long_fig4b %>% mutate( dilution_fraction_numeric = as.numeric(sub("%","", kd_fraction_label)) /100 ) lm_results_fig4b_table <- df_fig4b_linear_calc %>% group_by(primer_assay) %>% do(tidy(lm(fold_change_value ~dilution_fraction_numeric, data = .))) %>% filter(term == "dilution_fraction_numeric")%>% select(Assay = primer_assay, Slope = estimate, P_Value_Slope = p.value) r2_results_fig4b_table <- df_fig4b_linear_calc %>% group_by(primer_assay) %>% do(glance(lm(fold_change_value ~dilution_fraction_numeric, data = .))) %>% select(Assay = primer_assay, R_Squared = r.squared) linearity_summary_fig4b_table <- left_join(lm_results_fig4b_table, r2_results_fig4b_table, by="Assay") cat("\nLinearity of raw fold change vs. KD fraction:\n") ## ## Linearity of raw fold change vs. KD fraction: 25
RNA_Event_Conventional Mean_Ctrl Mean_ALS Statistic P_Value DF ARHGAP32 (Conv.) 0 2.318 -2.297 0.062 5.882 UNC13A CE 0 0.147 -0.140 0.892 7.728 7 Supplementary Figure S1: GAPDH/RPLP1 Shifts After TDP43 Depletion 7.1 Supplementary Figure S1B and S1C: GAPDH mRNA Abundance Relative to RPLP1 This section analyzes the relative abundance of GAPDH mRNA (normalized to RPLP1) after TDP-43 depletion in differentiated SH-SY5Y cells and HEK293T cells, highlighting the instability of traditional housekeeping genes. # Load and process data for housekeeping gene analysis df_figs1bc <- read.csv("Rawdata/S1B-C.csv")%>% mutate( RE_GAPDH_vs_RPLP1 = 2^(-GAPDH), cell = factor(cell, levels = c("SH-SY5Y","HEK293T")), type = factor(type, levels = c("Scramble","shTDP43","siCtrl","siTDP43")) ) # Function to create plot for each cell type plot_by_cell_figs1bc <- function(cell_type_name) { df_sub_figs1bc <- df_figs1bc %>% filter(cell == cell_type_name) comp_list_figs1bc <- switch( cell_type_name, "SH-SY5Y" =list(c("Scramble","shTDP43")), "HEK293T" =list(c("siCtrl","siTDP43")) ) pval_figs1bc <- NA if (n_distinct(df_sub_figs1bc$type) >= 2) { relevant_types_figs1bc <- comp_list_figs1bc[[1]] test_data_figs1bc <- df_sub_figs1bc %>% filter(type %in% relevant_types_figs1bc) if(n_distinct(test_data_figs1bc$type) == 2&& nrow(test_data_figs1bc) >2) { test_result_figs1bc <- t.test(RE_GAPDH_vs_RPLP1 ~type, data = test_data_figs1bc, var.equal = FALSE) pval_figs1bc <- test_result_figs1bc$p.value } } pval_label_figs1bc <- if (is.na(pval_figs1bc)) "N/A" else if (pval_figs1bc <0.001)"p < 0.001" else paste0("p = ", signif(pval_figs1bc, 2)) max_y_figs1bc <- max(df_sub_figs1bc$RE_GAPDH_vs_RPLP1, na.rm = TRUE) min_y_figs1bc <- min(df_sub_figs1bc$RE_GAPDH_vs_RPLP1, na.rm = TRUE) y_pos_annot_figs1bc <- max_y_figs1bc +0.1 *(max_y_figs1bc -min_y_figs1bc) 32
STMN2 CE (Conv.) ARHGAP32 (Conv.) UNC13A CE Ctrl ALS Ctrl ALS Ctrl ALS −2 0 2 4 Group Log2 (Fold Change vs. RPLP1) Group Ctrl ALS Figure 12: Figure 5B. Relative expression of RNA events in ALS motor cortex vs. controls using conventional assays. Plots are faceted by RNA event. 33
fill_palette <- if(cell_type_name == "SH-SY5Y"){ c("Scramble"="#D3D3D3","shTDP43"="#E69F00") }else { c("siCtrl"="#D3D3D3","siTDP43"="#0072B2") } ggplot(df_sub_figs1bc, aes(x = type, y = RE_GAPDH_vs_RPLP1)) + stat_summary(fun = mean, geom = "crossbar",linewidth = 0.3,width = 0.6,color = "black")+ stat_summary(fun.data = mean_se, geom = "errorbar",width = 0.3,color = "black")+ geom_jitter(aes(fill = type), size = 2.5,stroke = 0.5,shape = 21, alpha = 0.7,width = 0.2,show.legend = FALSE)+ geom_signif( comparisons = comp_list_figs1bc, annotations = pval_label_figs1bc, y_position = y_pos_annot_figs1bc, tip_length = 0.01, textsize = 3, vjust = -0.2 )+ scale_fill_manual(values = fill_palette) + labs(x = "Condition",y = expression(GAPDH /RPLP1~(2^-(Delta*Ct))), title = cell_type_name) + scale_y_continuous(limits = c(0.25,0.95), expand = c(0,0), breaks = seq(0.3,0.9,by=0.1)) + theme_bw(base_size = 10)+ theme( axis.text.x = element_text(size = rel(1), angle=30,hjust=1,color = "black"), axis.text.y = element_text(size = rel(1), color = "black"), plot.title = element_text(size = rel(1.2), face = "bold",hjust = 0.5), axis.title.x = element_text(size = rel(1.1)), axis.title.y = element_text(size = rel(1.1)), legend.position = "none", panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black",linewidth = 0.5) ) } # Generate plots for SH-SY5Y and HEK293T plot_shsy5y_s1 <- plot_by_cell_figs1bc("SH-SY5Y") plot_hek293t_s1 <- plot_by_cell_figs1bc("HEK293T") print(plot_shsy5y_s1) print(plot_hek293t_s1) # Report t-test statistics in a table stats_s1bc_table <- df_figs1bc %>% group_by(cell) %>% filter(n_distinct(type) == 2)%>% do(tidy(t.test(RE_GAPDH_vs_RPLP1 ~type, data = ., var.equal = FALSE))) %>% 34
p < 0.001 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Scramble shTDP43 Condition GAPDH RPLP1 (2−(∆Ct)) SH−SY5Y Figure 13: Supplementary Figure S1B & S1C. GAPDH relative expression (2^-(ΔCt GAPDH - RPLP1)) after TDP-43 depletion in SH-SY5Y (top/left) and HEK293T (bottom/right) cells. Plots are generated sequentially. 35
p = 0.03 0.3 0.4 0.5 0.6 0.7 0.8 0.9 siCtrl siTDP43 Condition GAPDH RPLP1 (2−(∆Ct)) HEK293T Figure 14: Supplementary Figure S1B & S1C. GAPDH relative expression (2^-(ΔCt GAPDH - RPLP1)) after TDP-43 depletion in SH-SY5Y (top/left) and HEK293T (bottom/right) cells. Plots are generated sequentially. 36
select(Cell_Line = cell, Mean_Ctrl = estimate1, Mean_KD = estimate2, Statistic = statistic, P_Value = p.value, DF = parameter) cat("\nStatistical test results for Supplementary Figure S1B & S1C (TDP-43 KD vs Ctrl):\n") ## ## Statistical test results for Supplementary Figure S1B & S1C (TDP-43 KD vs Ctrl): kable(stats_s1bc_table, caption = "T-test results for GAPDH/RPLP1 ratio (Supplementary Figure S1B & S1C).",digits = 3) Table 9: T-test results for GAPDH/RPLP1 ratio (Supplementary Figure S1B & S1C). Cell_Line Mean_Ctrl Mean_KD Statistic P_Value DF SH-SY5Y 0.738 0.515 18.227 0.001 2.717 HEK293T 0.516 0.627 -3.826 0.030 3.064 37