Enhancing Maturation of Human Neuromuscular Organoids via Electrical Stimulation.

Code ↔ Paper

18 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 18 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted

[1] § Experimental Section › Functional Characterization of NMOs › NMO Contraction Analysis › Signal Extraction ↔ src/utils/time_series_extraction.py, lines 136–212 · score 0.82 · minimum border length, rotated image, sub region, cropped, angle, axis
[2] § Experimental Section › Functional Characterization of NMOs › NMO Contraction Analysis › Signal Analysis ↔ src/utils/pre_processing.py, the whole file · a weak match · score 0.78 · polynomial interpolation, de trending, pre processing, noise, smoothing, conversion
[3] § Experimental Section › Immunofluorescence Analysis of NMOs › Image Analysis ↔ Image-analysis/Col4A1 or Col5A3 area analysis.ipynb, lines 180–305 · score 0.73 · Col4A1, Col5A3, NMO region, muscle regions, segmented, ilastik
[4] § Experimental Section › Bulk RNA Sequencing › Alignment and Quantification ↔ Transcriptomics-analysis/EPS_NMOs_analysis_v2.Rmd, lines 17–64 · score 0.71 · featureCounts, Homo sapiens, FR, v2, matrix, GRCh38
[5] § Experimental Section › Immunofluorescence Analysis of NMOs › Image Analysis ↔ Image-analysis/PAX7-Ki67 analysis.ipynb, lines 172–315 · score 0.71 · PAX7 masks, muscle region mask, DAPI mask, Ki67, Ilastik, quantified
[6] § Experimental Section › Functional Characterization of NMOs › NMO Contraction Analysis › Signal Analysis ↔ src/2_time_series_analysis.ipynb, lines 120–185 · score 0.70 · de trending, pre processing, polynomial interpolation, smoothing, pixels, contractile
[7] § Results › Transcriptional Analysis Reveals Chronic Electrical Pulse Stimulation‐Induced Maturation of NMOs ↔ Transcriptomics-analysis/EPS_NMOs_analysis_v2.Rmd, lines 932–1012 · score 0.69 · cellular processes, gene expression, muscle contraction, dendrite, myelination, myogenesis
[8] § Experimental Section › Bulk RNA Sequencing › Gene Set Enrichment Analysis ↔ Transcriptomics-analysis/EPS_NMOs_analysis_v2.Rmd, lines 461–484 · score 0.68 · GO Biological Process, EnrichR, Hallmark, pathways, enrichment, genes
[9] § Experimental Section › Immunofluorescence Analysis of NMOs › Image Analysis ↔ Image-analysis/PAX7-Ki67 analysis.ipynb, lines 172–315 · score 0.67 · Ki67 masks, DAPI mask, region mask, segmented, ilastik, quantified
[10] § Experimental Section › Analysis of NMO Size ↔ src/3_border_and_area_computation.ipynb, lines 29–72 · score 0.66 · segmented image, Cellpose, centroid, diameter, border, models
[11] § Results › Dual Role of Chronic Electrical Pulse Stimulation in Skeletal Muscle Maturation ↔ Image-analysis/Col4A1 or Col5A3 area analysis.ipynb, lines 180–305 · score 0.65 · Col4A1, Col5A3, neural region, NMO region, muscle
[12] § Experimental Section › Immunofluorescence Analysis of NMOs › Image Analysis ↔ Image-analysis/MYOD analysis.ipynb, lines 99–129 · score 0.64 · MyoD1, muscle region mask, DAPI mask, Ilastik, signal, filtered
[13] § Results › Transcriptional Analysis Reveals Chronic Electrical Pulse Stimulation‐Induced Maturation of NMOs ↔ Transcriptomics-analysis/EPS_NMOs_analysis_v2.Rmd, lines 932–1012 · score 0.62 · muscle contraction, fold change, chronic EPS, myogenesis, axon, enriched
[14] § Experimental Section › Functional Characterization of NMOs › NMO Contraction Analysis › Signal Extraction ↔ src/utils/time_series_extraction.py, lines 136–212 · score 0.61 · minimum border length, border points, signal, frame, mask, threshold
[15] § Experimental Section › Immunofluorescence Analysis of NMOs › Image Analysis ↔ Image-analysis/BTX & Innervation analysis.ipynb, lines 163–182 · score 0.59 · TUBB3 mask, TUBB3 images, innervation, overlaying, ilastik, NMJ
[16] § Experimental Section › Immunofluorescence Analysis of NMOs › Image Analysis ↔ Image-analysis/SOX1-Ki67 analysis.ipynb, lines 111–126 · score 0.57 · neural region mask, DAPI mask, SOX1, Ki67, ilastik, filtered
[17] § Experimental Section › Immunofluorescence Analysis of NMOs › Image Analysis ↔ Image-analysis/ChAT area analysis.ipynb, lines 140–280 · score 0.55 · ChAT, NMO region, segmented, ilastik, DAPI, quantified
[18] § Experimental Section › Bulk RNA Sequencing › Alignment and Quantification ↔ Transcriptomics-analysis/rna_seq_pipeline_v2.sh, lines 18–33 · score 0.53 · samtools, strandness, HISAT2, genome, v2, rna

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 · 1,532 lines · 71 KB · MIT · 4 matches

---
title: "Analysis Report of EPS effect on NMO transcriptome using bulk RNA-seq"
output:
html_document:
theme: journal
toc: true
toc_float:
collapsed: false
smooth_scroll: false
date: "2025-10-16"
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE,warning = F,message = F)
```
# Data preprocessing
Qualtiy assessment of the sequenced reads was performed using *FastQC*. The data passed all important quality checks. Accordingly, reads were aligned to the homo sapiens genome assembly GRCh38 using *HISAT2* and aligned reads were then quantified using *featureCounts*. The count matrix was subsequently annotated for downstream analysis and differential gene expression was performed using *DESeq2* and gene set enrichment analysis using *EnrichR*.
```{r include = F}
# load the required libraries
library(tidyverse)
library(DESeq2)
library(DT)
library(patchwork)
# define important functions
extract_genes <- function(x,pattern) x %>% filter(Term %in% grep(x=x$Term, pattern= pattern,value = T,ignore.case=T)) %>% dplyr::select(Genes) %>% unlist %>% unname %>% str_split(';') %>% unlist %>% unique
expression_colors <- function(a,b) scale_fill_gradientn(colors = c("#2166ac", "#67a9cf", "#fbe7a2", "#f4a582", "#b2182b"),
values = scales::rescale(seq(-a,a,by = b), from = c(-a,a)),
limits = c(-a,a),oob = scales::squish)
# read in the data
rf.summ <- read.delim("/fast/AG_Gouti/BulkRNAseq/output/Anthie_NMOs/A5014/all_samples_featurecounts_RF.txt.summary")
#rf <- rowMeans(rf.summ[,-1])
## plot of assigned and unassigned counts
rf.summ %>% pivot_longer(names_to = 'sample',values_to= 'read_count', cols = 2:12) %>%
mutate(sample_id = str_extract(sample, "RNA_\\d+")) %>%
ggplot(aes(y= Status, x= read_count))+theme_bw()+geom_col()+facet_wrap(~sample_id)
## plot of % of assigned and unassigned counts
ism <- cbind(Status = rf.summ[,1], apply(rf.summ[,-1],2,function(x) 100*x/sum(x)) %>% as.data.frame)
ism %>% pivot_longer(names_to = 'sample',values_to= '% of total reads', cols = 2:12) %>%
mutate(sample_id = str_extract(sample, "RNA_\\d+")) %>%
ggplot(aes(y= Status, x= `% of total reads`))+theme_bw()+geom_col()+facet_wrap(~sample_id)
# read the FR featureCounts output
cts <- read.delim('/fast/AG_Gouti/BulkRNAseq/output/Anthie_NMOs/A5014/all_samples_featurecounts_RF.txt',header = T,skip = 1,sep = "\t" )
# remove unnecessary columns
cts <- cts[,-2:-6]
# move gene IDs to become rownames
cts <- cts %>% column_to_rownames('Geneid') %>% as.matrix
# adjust colnames to run IDs
colnames(cts) <- str_extract(colnames(cts), "RNA_\\d+")
# define sample meta data
coldata <- data.frame(sample = paste0('RNA_',c(paste0('0',1:9),10,11)),
plate_no = c('P1648',rep('P1571',5),rep('P1575',5)),
timepoint = c('d30',rep(c('d40','d40','d60','d60','d60'), times = 2)),
line = 'TTNGFP',
group = c('Control',rep(c('Control','EPS1a','Control','EPS1b','EPS2'),times = 2)))
coldata <- coldata %>% column_to_rownames('sample')
```
```{r include=FALSE}
# We choose to use RF because of the sequencing chemistry that was used
# rearrange columns of count matrix to match rows of meta data
cts <- cts[,rownames(coldata)]
ncol(cts);nrow(coldata)
identical(rownames(coldata),colnames(cts))
# exclude genes with zero expression in all samples
cts_filtered <- cts[rowSums(cts) > 0,]
dim(cts);dim(cts_filtered)
```
```{r include=FALSE}
# Converting eensembl IDs to gene symbols
library(biomaRt)
# Connect to the Ensembl BioMart database for human
mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
# Convert from ensembl IDs to gene symbols
# Retrieve mappings along with an indicator of canonical transcripts
mapping <- getBM(attributes = c("hgnc_symbol", "ensembl_gene_id", "transcript_is_canonical"),
filters = "ensembl_gene_id",
values = rownames(cts_filtered),
mart = mart)
# Filter for canonical isoforms only
canonical_mapping <- mapping %>% filter(transcript_is_canonical == 1)
length(rownames(cts_filtered)); nrow(mapping);nrow(canonical_mapping)
unique(mapping$ensembl_gene_id) %>% length
# search for duplicated Ensembl IDs and resolve them
canonical_mapping[duplicated(canonical_mapping$ensembl_gene_id),]
# search for duplicated symbols and resolve them
canonical_mapping[duplicated(canonical_mapping$hgnc_symbol),]
kar = canonical_mapping$hgnc_symbol[duplicated(canonical_mapping$hgnc_symbol)]
kar = kar[kar != ""]
cat(kar,sep = '\n') # list them
# remove rows 4234, 12582, 27361, 39893
#canonical_mapping <- canonical_mapping[-c(23619, 31588),]
# search for Ensembl IDs without a matched gene symbol
# Those IDs represent novel genes that have yet no symbols
canonical_mapping %>% filter(hgnc_symbol == "") %>% nrow()
# an easier way of mapping ensembl IDs to hgnc symbols
features_metadata <- data.frame(ensembl_gene_id = rownames(cts_filtered)) %>%
left_join(mapping[,-3], by = 'ensembl_gene_id') %>% distinct()
dim(features_metadata)
identical(features_metadata %>% arrange(ensembl_gene_id) %>% dplyr::select(hgnc_symbol,ensembl_gene_id),
canonical_mapping[,-3] %>% arrange(ensembl_gene_id))
# Since both canonical_mapping & features_metadata are identical, we can use them interchangeably
features_metadata[duplicated(features_metadata$ensembl_gene_id),]
features_metadata[duplicated(features_metadata$hgnc_symbol),]
all(features_metadata[duplicated(features_metadata$hgnc_symbol),2] != '')
# idx <- which(features_metadata[duplicated(features_metadata$hgnc_symbol),2] != '')
# ism = features_metadata[duplicated(features_metadata$hgnc_symbol),]
# ism[idx,] %>% View()
```
# Sample metadata
The following table shows information about the samples included in this dataset. The organoid line, timepoint, plate number and treatment group are specified.
```{r}
coldata$replicate <- ifelse(coldata$plate_no == 'P1575', 'Rep2','Rep1')
coldata$group_desc <- case_match(coldata$group, "Control" ~ "Control", "EPS1a" ~ "Early-stage EPS", "EPS1b" ~ "Chronic EPS, stable","EPS2" ~ "Chronic EPS, increasing")
coldata %>% datatable(options = list(pageLength = 11))
```
# Clustering samples
## Gene expression heatmap
The following heatmap shows the expression values (after undergoing variance stabilizing transformation) of the 2000 most highly variable genes in the dataset. Both samples (columns) and genes (rows) are clustered using hierarchical clustering.
```{r include = F}
library(pheatmap)
head(coldata)
# create DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = cts_filtered,
colData = coldata,
design = ~ group + timepoint)
dds <- DESeq(dds)
resultsNames(dds)
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
ntd <- normTransform(dds)
# select 2000 most highly expressed genes
#select <- order(rowMeans(counts(dds,normalized=TRUE)), decreasing=TRUE)[1:2000]
# select 2000 most highly variable genes
select <- rowSds(counts(dds,normalized= TRUE)) %>% sort(decreasing = T) %>% head(2000) %>% names()
df <- as.data.frame(colData(dds)[,c("group","timepoint","plate_no")])
pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
cluster_cols=TRUE, annotation_col=df)
```
```{r}
pheatmap(assay(vsd)[select,], cluster_rows=TRUE, show_rownames=FALSE,
cluster_cols=TRUE, annotation_col=df)
```
```{r include = F}
pheatmap(assay(rld)[select,], cluster_rows=FALSE, show_rownames=FALSE,
cluster_cols=TRUE, annotation_col=df)
```
## Heatmap of sample-to-sample distances
The following is a heatmap of the sample-to-sample distance matrix where the smaller the distance, the closer the two samples are to each other in terms of their transcriptomic profile. The clustering shows the same pattern as in the expression heatmap above.
```{r fig.dim=c(700/96,400/96),fig.dpi=300}
sampleDists <- dist(t(assay(vsd)))
library("RColorBrewer")
sampleDistMatrix <- as.matrix(sampleDists)
# rownames(sampleDistMatrix) <- colnames(sampleDistMatrix) <- paste(vsd$timepoint,vsd$plate_no,vsd$group,sep="-")
rownames(sampleDistMatrix) <- paste(vsd$timepoint,vsd$group_desc,vsd$replicate, sep = '_')
colnames(sampleDistMatrix) <- NULL
# colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
# pheatmap(sampleDistMatrix,
# clustering_distance_rows=sampleDists,
# clustering_distance_cols=sampleDists,
# col=colors,
# breaks = seq(20,60,length.out = 256))
colors <- colorRampPalette(c("#b2182b", "#f4a582","#fbe7a2"))(255)
#colors <- colorRampPalette(c("#b2182b","#fbe7a2"))(255)
#colors <- viridis::viridis(255) %>% rev
plt <- pheatmap(sampleDistMatrix,
clustering_distance_rows=sampleDists,
clustering_distance_cols=sampleDists,
col=colors,
breaks = seq(20,60,length.out = 256))
plt
# ism = as.vector(sampleDistMatrix)
# min(ism[ism!=0])
# mat_log <- log10(sampleDistMatrix + 1)
#
# # Define log-spaced breaks (using the log-transformed min and max)
# log_breaks <- seq(min(mat_log), max(mat_log), length.out = 256)
#
# # Now map these back to the original (non-log) scale to set the `breaks` param
# orig_breaks <- 10^log_breaks
```
## Principal component analysis
```{r include = F}
plotPCA(vsd, intgroup=c("timepoint", "plate_no","group"))+theme_bw()
```
The following PCA plot shows the different samples annotated by timepoint and plate number. The first PC nicely captures the variation along the timepoints from d30 through d40 to d60. This shows that the timepoint is the single most important determinant of the transcriptomic status.
```{r fig.dim=c(6.25,4.2)}
pcaData <- plotPCA(vsd, intgroup=c("timepoint", "plate_no","group"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color= timepoint, shape= plate_no)) +
geom_point(size=3) + theme_bw()+
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) #+
#coord_fixed()
```
The following PCA plot shows the different samples annotated by group and plate number.
```{r fig.dim=c(6.25,4.2)}
pcaData <- plotPCA(vsd, intgroup=c("timepoint", "plate_no","group"), returnData=TRUE)
ggplot(pcaData, aes(PC1, PC2, color= group.1, shape= plate_no)) +
geom_point(size=3) + theme_bw()+
scale_color_manual(values = c('darkolivegreen4', 'gold','tomato','dodgerblue'))+
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) #+
#coord_fixed()
```
The following PCA plot shows the different samples annotated by timepoint, group and plate number.
```{r fig.dim=c(7.8,4.2)}
pcaData <- plotPCA(vsd, intgroup=c("timepoint", "plate_no","group_desc"), returnData=TRUE)
pcaData$group_desc <- factor(pcaData$group_desc, levels = c('Control','Early-stage EPS','Chronic EPS, stable','Chronic EPS, increasing'))
#pcaData$plate_no <- factor(pcaData$plate_no, levels = c('P1648','P1571','P1575'))
plt <- ggplot(pcaData, aes(PC1, PC2, color= group_desc, shape= plate_no, label = timepoint)) +
geom_point(size=3) + theme_bw()+
ggrepel::geom_text_repel(size =5)+
scale_color_manual(values = c('darkolivegreen4', 'gold','tomato','dodgerblue'))+
labs(color = 'Group')+
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) +
theme(axis.text = element_text(size = 12, color = 'black'),
axis.title = element_text(size = 14),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14))
#coord_fixed()
plt
```
# Differential gene expression
Differential gene expression analysis was performed using DESeq2 where gene expression was modeled as a function of the group (EPS vs Control) and the plate number (P1575 vs P1571) to take into account the variability between the two replicates. A 3rd replicate is recommended to improve the statistical robustness of this analysis.
## D40 comparisons (EPS1a vs Control)
```{r include = F}
my.coldata <- coldata %>% filter(timepoint == 'd40')
my.counts <- cts_filtered[,rownames(my.coldata)]
# define the test variable and its reference level
my.coldata$group <- factor(my.coldata$group, levels = c('Control','EPS1a'))
my.coldata$plate_no <- factor(my.coldata$plate_no, levels = c('P1571','P1575'))
table(rownames(my.coldata),my.coldata$group)
table(rownames(my.coldata),my.coldata$plate_no)
# create DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = my.counts,
colData = my.coldata,
design = ~ group + plate_no)
# pre-filter low-count genes
smallestGroupSize <- 2
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep,]
# run DGE analysis
dds <- DESeq(dds)
resultsNames(dds)
res1a <- results(dds,name = "group_EPS1a_vs_Control") %>% as.data.frame
#res <- results(dds,name = "plate_no_P1575_vs_P1571") %>% as.data.frame
#res
# apply Log fold change shrinkage for visualization and ranking
#resLFC <- lfcShrink(dds, coef= "group_EPS1a_vs_Control", type="apeglm") %>% as.data.frame
#resLFC <- lfcShrink(dds, coef= "plate_no_P1575_vs_P1571", type="apeglm") %>% as.data.frame
#resLFC
# inspect results with and without LFC shrinkage
res1a %>% filter(padj < 0.05) %>% nrow
res1a %>% filter(padj < 0.05 & log2FoldChange > 0) %>% nrow
res1a %>% filter(padj < 0.05 & log2FoldChange < 0) %>% nrow
# resLFC %>% filter(padj < 0.05) %>% nrow
# resLFC %>% filter(padj < 0.05 & log2FoldChange > 0) %>% nrow
# resLFC %>% filter(padj < 0.05 & log2FoldChange < 0) %>% nrow
#View(res)
```
The following table shows the 1000 most significant DEGs arranged by P-value.
```{r}
res1a <- res1a %>% rownames_to_column('ensembl_gene_id')
res1a <- left_join(res1a,canonical_mapping[,c('ensembl_gene_id','hgnc_symbol')], by = 'ensembl_gene_id')
res1a <- res1a[,colnames(res1a)[c(1,8,2:7)]]
res1a %>% arrange(pvalue) %>% head(1000) %>% datatable(options = list(pageLength= 10))
```
```{r}
dn <- res1a %>% filter(padj < 0.05 & log2FoldChange < 0) %>% mutate(hgnc_symbol = ifelse(hgnc_symbol == "","Novel_Gene",hgnc_symbol)) %>% mutate(id_symbol = paste(ensembl_gene_id,hgnc_symbol,sep ='_')) %>% dplyr::select(id_symbol) %>% unlist %>% unname
up <- res1a %>% filter(padj < 0.05 & log2FoldChange > 0) %>% mutate(hgnc_symbol = ifelse(hgnc_symbol == "","Novel_Gene",hgnc_symbol)) %>% mutate(id_symbol = paste(ensembl_gene_id,hgnc_symbol,sep ='_')) %>% dplyr::select(id_symbol) %>% unlist %>% unname
```
- 4 genes were significantly upregulated at an FDR-adjusted P-value less than 0.05, namely: `r up`
- 4 genes were downregulated at the adjusted P-value threshold, namely: `r dn`
## D60 comparisons
```{r include=FALSE}
my.coldata <- coldata %>% filter(timepoint == 'd60')
my.counts <- cts_filtered[,rownames(my.coldata)]
# define the test variable and its reference level
my.coldata$group <- factor(my.coldata$group, levels = c('Control','EPS1b','EPS2'))
my.coldata$plate_no <- factor(my.coldata$plate_no, levels = c('P1571','P1575'))
table(rownames(my.coldata),my.coldata$group)
table(rownames(my.coldata),my.coldata$plate_no)
# create DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = my.counts,
colData = my.coldata,
design = ~ group + plate_no)
# pre-filter low-count genes
smallestGroupSize <- 2
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep,]
# run DGE analysis
dds <- DESeq(dds)
resultsNames(dds)
```
### EPS1b vs Control
```{r include = F}
res1b <- results(dds,name = "group_EPS1b_vs_Control") %>% as.data.frame()
#res2 <- results(dds,name = "group_EPS2_vs_Control") %>% as.data.frame()
#res
# apply Log fold change shrinkage for visualization and ranking
#resLFC <- lfcShrink(dds, coef= "group_EPS1b_vs_Control", type="apeglm") %>% as.data.frame()
#resLFC <- lfcShrink(dds, coef= "group_EPS2_vs_Control", type="apeglm") %>% as.data.frame()
#resLFC
# inspect results with and without LFC shrinkage
res1b %>% filter(padj < 0.05) %>% nrow
res1b %>% filter(padj < 0.05 & log2FoldChange > 0) %>% nrow
res1b %>% filter(padj < 0.05 & log2FoldChange < 0) %>% nrow
# resLFC %>% filter(padj < 0.05) %>% nrow
# resLFC %>% filter(padj < 0.05 & log2FoldChange > 0) %>% nrow
# resLFC %>% filter(padj < 0.05 & log2FoldChange < 0) %>% nrow
#View(res)
#sig.genes <- res %>% filter(padj < 0.05 & log2FoldChange > 0) %>% rownames
```
The following table shows the 1000 most significant DEGs arranged by P-value.
```{r}
res1b <- res1b %>% rownames_to_column('ensembl_gene_id')
res1b <- left_join(res1b,canonical_mapping[,c('ensembl_gene_id','hgnc_symbol')], by = 'ensembl_gene_id')
res1b <- res1b[,colnames(res1b)[c(1,8,2:7)]]
res1b %>% arrange(pvalue) %>% head(1000) %>% datatable(options = list(pageLength= 10))
```
- 46 genes were significantly upregulated and no genes significantly downregulated at an FDR-adjust P-value less than 0.05. The upregulated genes are listed in the following table.
```{r}
#gene_symbols %>% datatable(options = list(pageLength = 10))
res1b %>% filter(padj < 0.05 & log2FoldChange > 0) %>% datatable(options = list(pageLength = 10))
```
### EPS2 vs Control
```{r include=FALSE}
res2 <- results(dds,name = "group_EPS2_vs_Control") %>% as.data.frame()
# apply Log fold change shrinkage for visualization and ranking
#resLFC <- lfcShrink(dds, coef= "group_EPS2_vs_Control", type="apeglm") %>% as.data.frame()
# inspect results with and without LFC shrinkage
res2 %>% filter(padj < 0.05) %>% nrow
res2 %>% filter(padj < 0.05 & log2FoldChange > 0) %>% nrow
res2 %>% filter(padj < 0.05 & log2FoldChange < 0) %>% nrow
# resLFC %>% filter(padj < 0.05) %>% nrow
# resLFC %>% filter(padj < 0.05 & log2FoldChange > 0) %>% nrow
# resLFC %>% filter(padj < 0.05 & log2FoldChange < 0) %>% nrow
up <- res2 %>% filter(padj < 0.05 & log2FoldChange > 0) %>% rownames
dn <- res2 %>% filter(padj < 0.05 & log2FoldChange < 0) %>% rownames
```
The following table shows the 1000 most significant DEGs arranged by P-value.
```{r}
res2 <- res2 %>% rownames_to_column('ensembl_gene_id')
res2 <- left_join(res2,canonical_mapping[,c('ensembl_gene_id','hgnc_symbol')], by = 'ensembl_gene_id')
res2 <- res2[,colnames(res2)[c(1,8,2:7)]]
res2 %>% arrange(pvalue) %>% head(1000) %>% datatable(options = list(pageLength= 10))
```
- A total of 52 genes were differentially expressed at an FDR-adjusted P-value less than 0.05; 48 of which were significantly upregulated and 4 were significantly downregulated. Those are listed in the following table.
```{r}
#gene_symbols2 %>% datatable(options = list(pageLength = 10))
res2 %>% filter(padj < 0.05) %>% datatable(options = list(pageLength = 10))
```
# Functional enrichment analysis
## EnrichR
```{r include=FALSE}
# Functional enrichment
library(enrichR)
websiteLive <- getOption("enrichR.live")
if (websiteLive) {
listEnrichrSites()
setEnrichrSite("Enrichr") # Human genes
}
# if (websiteLive) {
# dbs <- listEnrichrDbs()
# head(dbs)
# }
dbs <- c('MSigDB_Hallmark_2020','GO_Biological_Process_2023','GO_Cellular_Component_2023','GO_Molecular_Function_2023',
'WikiPathways_2024_Human','Reactome_Pathways_2024','KEGG_2021_Human','Elsevier_Pathway_Collection',
'OMIM_Disease','OMIM_Expanded','Jensen_DISEASES','DisGeNET')
```
Gene set or functional enrichment analysis was performed using the statistical tests available in *EnrichR*. Gene sets of 12 different databases were tested, namely:
```{r}
cat(dbs,sep = '\n')
```
Only significant DEGs with an FDR-adjusted P-value less than 0.05 were tested for enrichment.
### D60 comparisons
#### EPS1b vs Control
```{r include = F}
g = res1b %>% filter(padj < 0.05 & log2FoldChange > 0) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
if (websiteLive) {
#enriched <- enrichr(gene_symbols$hgnc_symbol, dbs)
enriched <- enrichr(g, dbs)
}
enriched.df <- list_rbind(enriched[sapply(enriched, nrow) != 0], names_to = 'database')
```
The following table shows the 1000 most significantly enriched gene sets arranged by P-value. These represent functions that are upregulated in EPS1b vs control.
```{r}
enriched.df %>% arrange(P.value) %>% head(1000) %>%
datatable(options = list(pageLength= 10))
```
#### EPS2 vs Control
```{r include = F}
# up.symbols <- gene_symbols2 %>% filter(ensembl_gene_id %in% up)%>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
g = res2 %>% filter(padj < 0.05 & log2FoldChange > 0) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
if (websiteLive) {
#enriched <- enrichr(up.symbols, dbs)
enriched <- enrichr(g, dbs)
}
enriched.df <- list_rbind(enriched[sapply(enriched, nrow) != 0], names_to = 'database')
```
The following table shows the 1000 most significantly enriched gene sets arranged by P-value. These represent functions that are upregulated in EPS2 vs control.
```{r}
enriched.df %>% arrange(P.value) %>% head(1000) %>%
datatable(options = list(pageLength= 10))
```
## FGSEA
In contract to *EnrichR*, *FGSEA* performs functional enrichment using a ranked genome-wide gene list where genes are ranked from the most significantly upregulated to the most significantly downregulated or vice versa based on log2 fold change or some measure of statistical significance. This approach does not require the use of thresholds/cutoffs for defining statistically significant genes prior to the enrichment analysis.
```{r include = FALSE}
library(fgsea)
# load gene sets
h <- rjson::fromJSON(file = '/fast/AG_Gouti/snRNAseq/MSigDB_GeneSets/h.all.v2025.1.Hs.json.txt')
go.bp <- rjson::fromJSON(file = '/fast/AG_Gouti/snRNAseq/MSigDB_GeneSets/c5.go.bp.v2025.1.Hs.json.txt')
h.gs <- lapply(h, function(x) x$geneSymbols)
go.bp.gs <- lapply(go.bp, function(x) x$geneSymbols)
gs <- c(h.gs,go.bp.gs)
```
### D40 comparisons (EPS1a vs Control)
```{r include=FALSE}
dim(res1a); dim(res1a[complete.cases(res1a),])
df <- res1a[complete.cases(res1a),] %>% filter(hgnc_symbol != '')
sum(duplicated(df$hgnc_symbol))
df$hgnc_symbol[duplicated(df$hgnc_symbol)]
#df %>% filter(hgnc_symbol %in% df$hgnc_symbol[duplicated(df$hgnc_symbol)]) %>% View()
# df <- res1a[complete.cases(res1a),] %>% rownames_to_column("ensembl_gene_id") %>% left_join(mapping[,c("ensembl_gene_id","hgnc_symbol")], by = 'ensembl_gene_id') %>% filter(hgnc_symbol != '')
# -log10(min(df$pvalue[df$pvalue > 0]) + 1e-10)
# df <- df %>% mutate(metric = sign(log2FoldChange) * -log10(pvalue + 1e-10)) %>%
# arrange(metric) %>% dplyr::select(hgnc_symbol,metric)
# ranks <- df$metric; names(ranks) <- df$hgnc_symbol
# Apply a function to select the ID with more transcripts
ls <- split(x = df, f = df$hgnc_symbol)
df <- lapply(ls, function(x){
data.frame(stat = x$stat[which.max(x$baseMean)])
}) %>% list_rbind(names_to = 'hgnc_symbol') %>% arrange(stat)
ranks <- df$stat; names(ranks) <- df$hgnc_symbol
fgseaRes <- fgsea(pathways = gs,
stats = ranks,
eps = 0.0,
nPermSimple = 50000,
minSize = 15,
maxSize = 500)
```
The following table displays the gene sets significantly enriched at an FDR-adjusted P-value < 0.05 and arranged by P-value.
```{r}
fgseaRes %>% filter(padj < 0.05) %>% arrange(pval) %>% datatable(options = list(pageLength= 10))
```
### D60 comparisons
#### EPS1b vs Control
```{r include=FALSE}
dim(res1b); dim(res1b[complete.cases(res1b),])
df <- res1b[complete.cases(res1b),] %>% filter(hgnc_symbol != '')
sum(duplicated(df$hgnc_symbol))
df$hgnc_symbol[duplicated(df$hgnc_symbol)]
# -log10(min(df$pvalue[df$pvalue > 0]) + 1e-10)
# df <- df %>% mutate(metric = sign(log2FoldChange) * -log10(pvalue + 1e-10)) %>%
# arrange(metric) %>% dplyr::select(hgnc_symbol,metric)
# ranks <- df$metric; names(ranks) <- df$hgnc_symbol
# Apply a function to select the ID with more transcripts
ls <- split(x = df, f = df$hgnc_symbol)
df <- lapply(ls, function(x){
data.frame(stat = x$stat[which.max(x$baseMean)])
}) %>% list_rbind(names_to = 'hgnc_symbol') %>% arrange(stat)
ranks <- df$stat; names(ranks) <- df$hgnc_symbol
fgseaRes <- fgsea(pathways = gs,
stats = ranks,
eps = 0.0,
nPermSimple = 50000,
minSize = 15,
maxSize = 500)
```
The following table displays the gene sets significantly enriched at an FDR-adjusted P-value < 0.05 and arranged by P-value.
```{r}
fgseaRes %>% filter(padj < 0.05) %>% arrange(pval) %>% datatable(options = list(pageLength= 10))
```
#### EPS2 vs Control
```{r include=FALSE}
dim(res2); dim(res2[complete.cases(res2),])
df <- res2[complete.cases(res2),] %>% filter(hgnc_symbol != '')
sum(duplicated(df$hgnc_symbol))
df$hgnc_symbol[duplicated(df$hgnc_symbol)]
# -log10(min(df$pvalue[df$pvalue > 0]) + 1e-10)
# df <- df %>% mutate(metric = sign(log2FoldChange) * -log10(pvalue + 1e-10)) %>%
# arrange(metric) %>% dplyr::select(hgnc_symbol,metric)
# ranks <- df$metric; names(ranks) <- df$hgnc_symbol
# Apply a function to select the ID with more transcripts
ls <- split(x = df, f = df$hgnc_symbol)
df <- lapply(ls, function(x){
data.frame(stat = x$stat[which.max(x$baseMean)])
}) %>% list_rbind(names_to = 'hgnc_symbol') %>% arrange(stat)
ranks <- df$stat; names(ranks) <- df$hgnc_symbol
fgseaRes <- fgsea(pathways = gs,
stats = ranks,
eps = 0.0,
nPermSimple = 50000,
minSize = 15,
maxSize = 500)
```
The following table displays the gene sets significantly enriched at an FDR-adjusted P-value < 0.05 and arranged by P-value.
```{r}
fgseaRes %>% filter(padj < 0.05) %>% arrange(pval) %>% datatable(options = list(pageLength= 10))
```
# Visualization of signifcant genes
## D60 comparisons
### Bar graphs
DEGs with an adjusted P-value < 0.05 are visualized in a descending order of statistical significance.
```{r include = F}
# create DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = cts_filtered,
colData = coldata,
design = ~ group + timepoint)
#dds <- DESeq(dds)
#resultsNames(dds)
# Visualize genes of interest
ntd <- normTransform(dds)
norm.cts <- assay(ntd)
# library(vsn)
# meanSdPlot(assay(ntd))
#head(assay(ntd),3)
#head(my.counts,3)
p = 0.05 ; fc = 1
# define significant DEGs
up1b = res1b %>% filter(padj < p & log2FoldChange > log2(fc)) %>% arrange(pvalue) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
dn1b = res1b %>% filter(padj < p & log2FoldChange < -log2(fc)) %>% arrange(pvalue) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
up2 = res2 %>% filter(padj < p & log2FoldChange > log2(fc)) %>% arrange(pvalue) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
dn2 = res2 %>% filter(padj < p & log2FoldChange < -log2(fc)) %>% arrange(pvalue) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
```
#### Genes significantly upregulated by EPS1b
```{r fig.dim= c(1000/96,700/96), fig.dpi= 300}
g = up1b
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
#df <- cbind(coldata,norm.cts[rownames(norm.cts) %in% g,] %>% t)
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% filter(timepoint == 'd60' & group != 'EPS2') %>%
tidyr::pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 7:(6+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
# ism <- res1b %>% filter(padj < 0.05) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# identical(unique(df$hgnc_symbol), ism)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol))
df %>% ggplot(aes(x = group, y = exp, fill= plate_no))+geom_col(width = 0.45, position = position_dodge())+theme_bw()+facet_wrap(~hgnc_symbol)+
labs(x = '', y = 'Normalized expression')
```
#### Genes significantly upregulated by EPS2
```{r fig.dim= c(1000/96,700/96), fig.dpi= 300}
g = up2
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% filter(timepoint == 'd60' & group != 'EPS1b') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 7:(6+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
# ism <- res2 %>% filter(padj < 0.05 & log2FoldChange > 0) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# identical(unique(df$hgnc_symbol), ism)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol))
df %>% ggplot(aes(x = group, y = exp, fill= plate_no))+geom_col(width = 0.45, position = position_dodge())+theme_bw()+facet_wrap(~hgnc_symbol)+
labs(x = '', y = 'Normalized expression')
```
#### Genes significantly upregulated by both EPS1b & EPS2
```{r fig.dim=c(1100/96,700/96), fig.dpi=300}
g = intersect(up1b,up2)
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
a <- res1b %>% filter(ensembl_gene_id %in% g) %>% dplyr::select(hgnc_symbol,stat)
b <- res2 %>% filter(ensembl_gene_id %in% g) %>% dplyr::select(hgnc_symbol,stat)
g2 <- inner_join(a,b, by = "hgnc_symbol") %>% mutate(mean.stat = (stat.x + stat.y)/2) %>% arrange(desc(mean.stat)) %>%
dplyr::select(hgnc_symbol) %>% unlist %>% unname
df <- cbind(coldata,norm.cts[g,] %>% t)
# head(df)
# df %>% rownames_to_column("sample") %>%
# pivot_longer(names_to = "gene",values_to = "exp",cols = 7:48) %>%
# ggplot(aes(x = sample, y = exp, fill= group))+geom_col(width = 0.6)+theme_bw()+facet_wrap(~gene)
# df %>% rownames_to_column("sample") %>%
# pivot_longer(names_to = "gene",values_to = "exp",cols = 6:9) %>%
# ggplot(aes(x = group, y = exp,color = plate_no, group = plate_no))+geom_point(size = 1.5)+geom_line()+theme_bw()+facet_wrap(~gene)
df <- df %>% filter(timepoint == 'd60') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 7:(6+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = g2)
df %>% ggplot(aes(x = group, y = exp, fill= plate_no))+geom_col(width = 0.45, position = position_dodge())+theme_bw()+facet_wrap(~hgnc_symbol)+
labs(x = '', y = 'Normalized expression')#+
#theme(axis.text.x = element_text(angle = 30, vjust = 0.7))
```
#### Genes significantly downregulated by EPS2
```{r fig.dim=c(700/96, 240/96), fig.dpi=300}
g = dn2
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% filter(timepoint == 'd60' & group != 'EPS1b') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 7:(6+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
# ism <- res2 %>% filter(padj < 0.05 & log2FoldChange < 0) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# identical(unique(df$hgnc_symbol), ism)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol))
df %>% ggplot(aes(x = group, y = exp, fill= plate_no))+geom_col(width = 0.45, position = position_dodge())+theme_bw()+facet_wrap(~hgnc_symbol)+
labs(x = '', y = 'Normalized expression')+
theme(axis.text.x = element_text(angle = 0, vjust = 0.7))
```
### Time course plots
#### Genes significantly upregulated by EPS1b
```{r fig.dim= c(1200/96,750/96), fig.dpi= 300}
g = up1b
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
#df <- cbind(coldata,norm.cts[rownames(norm.cts) %in% g,] %>% t)
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% rownames_to_column('sample') %>%
filter(group != 'EPS2') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
# ism <- res1b %>% filter(padj < 0.05) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# identical(unique(df$hgnc_symbol), ism)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol))
#df$group <- factor(df$group, levels = c('Control','EPS1a','EPS1b','EPS2'))
#df$group2 <- case_match(df$group, c('EPS1a','EPS1b') ~ 'EPS', .default = df$group)
df$group2 <- ifelse(df$group == 'Control','Control','EPS stable')
#df$group3 <- paste(df$group2,df$plate_no, sep = '_')
df %>% ggplot(aes(x = timepoint, y = exp, color = group2, shape = plate_no))+
#geom_col(width = 0.5,position = position_dodge())+
geom_point()+
geom_line(data = df %>% filter(plate_no != 'P1575' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(plate_no != 'P1571' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_03','RNA_05')), aes(x = timepoint, y = exp, group = 1), color = 'tomato')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_08','RNA_10')), aes(x = timepoint, y = exp, group = 1), color = 'tomato')+
theme_bw()+facet_wrap(~hgnc_symbol,scales = 'free_y')+
#scale_color_manual(values = c('darkolivegreen4', 'tomato','dodgerblue','violet'))+
scale_color_manual(values = c('darkolivegreen4', 'tomato'))+
labs(x = '', y = 'Normalized expression')
```
#### Genes significantly upregulated by EPS2
```{r fig.dim= c(1300/96,750/96), fig.dpi= 300}
g = up2
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
#df <- cbind(coldata,norm.cts[rownames(norm.cts) %in% g,] %>% t)
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% rownames_to_column('sample') %>%
filter(group != 'EPS1b') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
# ism <- res1b %>% filter(padj < 0.05) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# identical(unique(df$hgnc_symbol), ism)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol))
#df$group <- factor(df$group, levels = c('Control','EPS1a','EPS1b','EPS2'))
df$group2 <- ifelse(df$group == 'Control','Control','EPS stable/increasing')
#df$group3 <- paste(df$group2,df$plate_no, sep = '_')
df %>% ggplot(aes(x = timepoint, y = exp, color = group2, shape = plate_no))+
#geom_col(width = 0.5,position = position_dodge())+
geom_point()+
geom_line(data = df %>% filter(plate_no != 'P1575' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(plate_no != 'P1571' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_03','RNA_06')) , aes(x = timepoint, y = exp, group = 1), color = 'dodgerblue')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_08','RNA_11')), aes(x = timepoint, y = exp, group = 1), color = 'dodgerblue')+
theme_bw()+facet_wrap(~hgnc_symbol,scales = 'free_y')+
#scale_color_manual(values = c('darkolivegreen4', 'tomato','dodgerblue','violet'))+
scale_color_manual(values = c('darkolivegreen4', 'dodgerblue'))+
labs(x = '', y = 'Normalized expression')
```
#### Genes significantly upregulated by both EPS1b & EPS2
```{r fig.dim=c(1150/96,700/96), fig.dpi=300}
g = intersect(up1b,up2)
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
a <- res1b %>% filter(ensembl_gene_id %in% g) %>% dplyr::select(hgnc_symbol,stat)
b <- res2 %>% filter(ensembl_gene_id %in% g) %>% dplyr::select(hgnc_symbol,stat)
g2 <- inner_join(a,b, by = "hgnc_symbol") %>% mutate(mean.stat = (stat.x + stat.y)/2) %>% arrange(desc(mean.stat)) %>%
dplyr::select(hgnc_symbol) %>% unlist %>% unname
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% rownames_to_column('sample') %>%
#filter(group != 'EPS1b') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
# ism <- res1b %>% filter(padj < 0.05) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# identical(unique(df$hgnc_symbol), ism)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = g2)
df$group <- factor(df$group, levels = c('Control','EPS1a','EPS1b','EPS2'))
#df$group2 <- case_match(df$group, c('EPS1a','EPS1b') ~ 'EPS', .default = df$group)
#df$group2 <- ifelse(df$group == 'Control','Control','EPS stable,increasing')
#df$group3 <- paste(df$group2,df$plate_no, sep = '_')
df %>% ggplot(aes(x = timepoint, y = exp, color = group, shape = plate_no))+
#geom_col(width = 0.5,position = position_dodge())+
geom_point()+
geom_line(data = df %>% filter(plate_no != 'P1575' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(plate_no != 'P1571' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_03','RNA_05')), aes(x = timepoint, y = exp, group = 1), color = 'tomato')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_08','RNA_10')), aes(x = timepoint, y = exp, group = 1), color = 'tomato')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_03','RNA_06')), aes(x = timepoint, y = exp, group = 1), color = 'dodgerblue')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_08','RNA_11')), aes(x = timepoint, y = exp, group = 1), color = 'dodgerblue')+
theme_bw()+facet_wrap(~hgnc_symbol,scales = 'free_y')+
scale_color_manual(values = c('darkolivegreen4', 'gold','tomato','dodgerblue'))+
#scale_color_manual(values = c('darkolivegreen3', 'dodgerblue'))+
labs(x = '', y = 'Normalized expression')
```
#### Genes significantly downregulated by EPS2
```{r fig.dim=c(800/96, 240/96), fig.dpi=300}
g = dn2
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
#df <- cbind(coldata,norm.cts[rownames(norm.cts) %in% g,] %>% t)
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% rownames_to_column('sample') %>%
filter(group != 'EPS1b') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
# ism <- res1b %>% filter(padj < 0.05) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# identical(unique(df$hgnc_symbol), ism)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol))
#df$group <- factor(df$group, levels = c('Control','EPS1a','EPS1b','EPS2'))
df$group2 <- ifelse(df$group == 'Control','Control','EPS stable,increasing')
#df$group3 <- paste(df$group2,df$plate_no, sep = '_')
df %>% ggplot(aes(x = timepoint, y = exp, color = group2, shape = plate_no))+
#geom_col(width = 0.5,position = position_dodge())+
geom_point()+
geom_line(data = df %>% filter(plate_no != 'P1575' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(plate_no != 'P1571' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_03','RNA_06')) , aes(x = timepoint, y = exp, group = 1), color = 'dodgerblue')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_08','RNA_11')), aes(x = timepoint, y = exp, group = 1), color = 'dodgerblue')+
theme_bw()+facet_wrap(~hgnc_symbol,scales = 'free_y')+
#scale_color_manual(values = c('darkolivegreen4', 'tomato','dodgerblue','violet'))+
scale_color_manual(values = c('darkolivegreen4', 'dodgerblue'))+
labs(x = '', y = 'Normalized expression')
```
### Heatmaps & Dot plots
DEGs with an adjusted P-value < 0.05 are mapped to the relevant biological functions and cellular processes they regulate.
#### EPS1b vs Control {.tabset .tabset-pills}
```{r include =FALSE}
#my.title = 'EPS1b vs Control comparison'
my.title = 'Chronic EPS, stable vs Control comparison'
g1 <- res1b %>% filter(padj < 0.05 & log2FoldChange > 0) %>% arrange(pvalue) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
g2 <- res1b %>% filter(padj < 0.05 & log2FoldChange > 0) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# perform enrichment analysis for significant genes
enriched <- enrichr(g2, dbs)
enriched.df <- list_rbind(enriched[sapply(enriched, nrow) != 0], names_to = 'database')
# Biological functions in Neurons & Muscles
gene.functions = list(`Neuronal functions` = extract_genes(enriched.df,'neuron'),
`Axon functions` = extract_genes(enriched.df,'axon'),
`Dendrite functions` = extract_genes(enriched.df,'dendrite'),
`Myelination`= extract_genes(enriched.df,'myelin'),
`Synaptic functions` = extract_genes(enriched.df,'synaps|synaptic transmission'),
`Myogenesis/Myotube formation`= extract_genes(enriched.df,'myogenesis|myotube|myoblast fusion'),
`Muscle contraction`= extract_genes(enriched.df,'muscle contract|muscular contract'),
`Angiogenesis`= extract_genes(enriched.df,'angiogenesis|vascular endothelial growth|blood vessel remodeling|blood vessel morphogenesis|blood vessel sprouting'),
`Endothelial cell functions`= extract_genes(enriched.df,'\\b(endothelial cell|endothelial cells)\\b'),
`Apoptosis/Cell death` = extract_genes(enriched.df,'apopto|cell death'),
`Mitochondrial functions` = extract_genes(enriched.df,'mitochondr|atp synthesis|oxidative phosphory'),
`Oxidative stress` = extract_genes(enriched.df,'oxidative stress')
)
# create gene-function data frame
gene.functions.df <- lapply(gene.functions, function(x) data.frame(gene = x)) %>% list_rbind(names_to = 'annotation')
my.genes <- unique(gene.functions.df$gene)
orf.genes <- grep(x = my.genes,pattern = 'orf',ignore.case = T,value = T)
my.genes[my.genes %in% orf.genes] <- str_to_title(orf.genes)
gene.functions.df$gene[gene.functions.df$gene %in% orf.genes] <- str_to_title(orf.genes)
# Gene expression heatmaps
g1sub <- canonical_mapping %>% filter(hgnc_symbol %in% unique(gene.functions.df$gene)) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
#mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g1)
#df <- cbind(coldata,norm.cts[g1,] %>% t %>% scale)
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g1sub)
df <- cbind(coldata,norm.cts[g1sub,] %>% t %>% scale)
df <- df %>% rownames_to_column('sample') %>% filter(timepoint == 'd60' & group != 'EPS2') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g1sub))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
df$group_plate <- paste(df$group,df$plate_no, sep = '_')
df$group_plate <- factor(df$group_plate, levels = c("Control_P1571","Control_P1575","EPS1b_P1571","EPS1b_P1575"))
df$group_replicate <- paste(df$group_desc,df$replicate, sep = ' ')
df$group_replicate <- factor(df$group_replicate, levels = c("Control Rep1","Control Rep2","Chronic EPS, stable Rep1","Chronic EPS, stable Rep2"))
#table(df$group_plate,df$group_replicate)
#df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol) %>% rev)
plt1v <- df %>% ggplot(aes(x= group_replicate, y= hgnc_symbol, fill= exp)) + # instead of x = group_plate
geom_tile() +
theme_bw()+
theme(panel.grid.major = element_blank())+
scale_x_discrete(position='bottom')+
#scale_fill_gradient2(low = 'dodgerblue', mid = 'white', high = 'tomato', midpoint = 0, limits = c(-1.2,1.2))+
expression_colors(1.5,0.75)+
labs(x= '',y= '',fill= 'Scaled expression', title = my.title)+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30,hjust = 1,size = 14),
axis.text.y = element_text(size = 12),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
plt1h <- df %>% ggplot(aes(y= group_replicate, x= hgnc_symbol, fill= exp)) + # instead of y = group_plate
geom_tile() +
theme_bw()+
theme(panel.grid.major = element_blank())+
scale_x_discrete(position='bottom')+
#scale_fill_gradient2(low = 'dodgerblue', mid = 'white', high = 'tomato', midpoint = 0, limits = c(-1.2,1.2))+
expression_colors(1.5,0.75)+
labs(x= '',y= '',fill= 'Scaled expression', title = my.title)+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust = 1,size = 12),
axis.text.y = element_text(size = 14),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
```
```{r include = F}
# Gene-Functions Dot plots
global_size_range <- range(c(-log10(res1b$pvalue), -log10(res2$pvalue)))
dt <- res1b %>% filter(hgnc_symbol %in% g2) %>% dplyr::select(hgnc_symbol,log2FoldChange,pvalue) %>% mutate(metric = -log10(pvalue))
# unique(gene.functions.df$gene)
# unique(dt$hgnc_symbol)
colnames(gene.functions.df)[2] <- 'hgnc_symbol'
dt <- left_join(gene.functions.df, dt, by = 'hgnc_symbol')
#summary(dt$log2FoldChange)
dt$annotation <- factor(dt$annotation, levels = names(gene.functions))
plt2v <- dt %>% ggplot(aes(x = annotation, y = hgnc_symbol, fill = log2FoldChange, size = metric))+geom_point(shape=21)+theme_bw()+expression_colors(1.5,0.75)+theme_bw()+
scale_size_continuous(limits = global_size_range, range = c(1,10))+
labs(x ='',y = '', size = 'statistical significance',title = my.title)+
guides(fill = guide_colourbar(order = 1),size = guide_legend(order = 2))+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust= 1,size = 14),
axis.text.y = element_text(size = 12),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
dt$annotation <- factor(dt$annotation, levels = names(gene.functions) %>% rev)
plt2h <- dt %>% ggplot(aes(y = annotation, x = hgnc_symbol, fill = log2FoldChange, size = metric))+geom_point(shape=21)+theme_bw()+expression_colors(1.5,0.75)+theme_bw()+
scale_size_continuous(limits = global_size_range, range = c(1,10))+
labs(x ='',y = '', size = 'statistical significance',title = my.title)+
guides(fill = guide_colourbar(order = 1),size = guide_legend(order = 2))+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust= 1,size = 12),
axis.text.y = element_text(size = 14),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
```
##### Vertical heatmap & dot plot
```{r fig.dim=c(1300/96,940/96), fig.dpi=300}
plt1v + plt2v + plot_layout(widths = c(4,12))
```
##### Horizontal heatmap & dot plot
```{r fig.dim=c(1350/96,700/96), fig.dpi=300}
plt1h / plt2h + plot_layout(heights = c(4,12))
```
#### {-}
#### EPS2 vs Control {.tabset .tabset-pills}
```{r include =FALSE}
#my.title = 'EPS2 vs Control comparison'
my.title = 'Chronic EPS, increasing vs Control comparison'
g1 <- res2 %>% filter(padj < 0.05 & log2FoldChange > 0) %>% arrange(pvalue) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
g2 <- res2 %>% filter(padj < 0.05 & log2FoldChange > 0) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# perform enrichment analysis for significant genes
enriched <- enrichr(g2, dbs)
enriched.df <- list_rbind(enriched[sapply(enriched, nrow) != 0], names_to = 'database')
# Biological functions in Neurons & Muscles
gene.functions = list(`Neuronal functions` = extract_genes(enriched.df,'neuron'),
`Axon functions` = extract_genes(enriched.df,'axon'),
`Dendrite functions` = extract_genes(enriched.df,'dendrite'),
`Myelination`= extract_genes(enriched.df,'myelin'),
`Synaptic functions` = extract_genes(enriched.df,'synaps|synaptic transmission'),
`Myogenesis/Myotube formation`= extract_genes(enriched.df,'myogenesis|myotube|myoblast fusion'),
`Muscle contraction`= extract_genes(enriched.df,'muscle contract|muscular contract'),
`Angiogenesis`= extract_genes(enriched.df,'angiogenesis|vascular endothelial growth|blood vessel remodeling|blood vessel morphogenesis|blood vessel sprouting'),
`Endothelial cell functions`= extract_genes(enriched.df,'\\b(endothelial cell|endothelial cells)\\b'),
`Apoptosis/Cell death` = extract_genes(enriched.df,'apopto|cell death'),
`Mitochondrial functions` = extract_genes(enriched.df,'mitochondr|atp synthesis|oxidative phosphory'),
`Oxidative stress` = extract_genes(enriched.df,'oxidative stress')
)
# create gene-function data frame
gene.functions.df <- lapply(gene.functions, function(x) data.frame(gene = x)) %>% list_rbind(names_to = 'annotation')
my.genes <- unique(gene.functions.df$gene)
orf.genes <- grep(x = my.genes,pattern = 'orf',ignore.case = T,value = T)
my.genes[my.genes %in% orf.genes] <- str_to_title(orf.genes)
gene.functions.df$gene[gene.functions.df$gene %in% orf.genes] <- str_to_title(orf.genes)
# Gene expression heatmaps
g1sub <- canonical_mapping %>% filter(hgnc_symbol %in% unique(gene.functions.df$gene)) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
#mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g1)
#df <- cbind(coldata,norm.cts[g1,] %>% t %>% scale)
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g1sub)
df <- cbind(coldata,norm.cts[g1sub,] %>% t %>% scale)
df <- df %>% rownames_to_column('sample') %>% filter(timepoint == 'd60' & group != 'EPS1b') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g1sub))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
df$group_plate <- paste(df$group,df$plate_no, sep = '_')
df$group_plate <- factor(df$group_plate, levels = c("Control_P1571","Control_P1575","EPS2_P1571","EPS2_P1575"))
df$group_replicate <- paste(df$group_desc,df$replicate, sep = ' ')
df$group_replicate <- factor(df$group_replicate, levels = c("Control Rep1","Control Rep2","Chronic EPS, increasing Rep1","Chronic EPS, increasing Rep2"))
#table(df$group_plate,df$group_replicate)
#df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol) %>% rev)
plt1v <- df %>% ggplot(aes(x= group_replicate, y= hgnc_symbol, fill= exp)) + # instead of x= group_plate
geom_tile() +
theme_bw()+
theme(panel.grid.major = element_blank())+
scale_x_discrete(position='bottom')+
#scale_fill_gradient2(low = 'dodgerblue', mid = 'white', high = 'tomato', midpoint = 0, limits = c(-1.2,1.2))+
expression_colors(1.5,0.75)+
labs(x= '',y= '',fill= 'Scaled expression', title = my.title)+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30,hjust = 1,size = 14),
axis.text.y = element_text(size = 12),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
plt1h <- df %>% ggplot(aes(y= group_replicate, x= hgnc_symbol, fill= exp)) + # instead of y= group_plate
geom_tile() +
theme_bw()+
theme(panel.grid.major = element_blank())+
scale_x_discrete(position='bottom')+
#scale_fill_gradient2(low = 'dodgerblue', mid = 'white', high = 'tomato', midpoint = 0, limits = c(-1.2,1.2))+
expression_colors(1.5,0.75)+
labs(x= '',y= '',fill= 'Scaled expression', title = my.title)+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust = 1,size = 12),
axis.text.y = element_text(size = 14),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
```
```{r include = F}
# Gene-Functions Dot plots
global_size_range <- range(c(-log10(res1b$pvalue), -log10(res2$pvalue)))
dt <- res2 %>% filter(hgnc_symbol %in% g2) %>% dplyr::select(hgnc_symbol,log2FoldChange,pvalue) %>% mutate(metric = -log10(pvalue))
# unique(gene.functions.df$gene)
# unique(dt$hgnc_symbol)
colnames(gene.functions.df)[2] <- 'hgnc_symbol'
dt <- left_join(gene.functions.df, dt, by = 'hgnc_symbol')
#summary(dt$log2FoldChange)
dt$annotation <- factor(dt$annotation, levels = names(gene.functions))
plt2v <- dt %>% ggplot(aes(x = annotation, y = hgnc_symbol, fill = log2FoldChange, size = metric))+geom_point(shape=21)+theme_bw()+expression_colors(1.5,0.75)+theme_bw()+
scale_size_continuous(limits = global_size_range, range = c(1,10))+
labs(x ='',y = '', size = 'statistical significance',title = my.title)+
guides(fill = guide_colourbar(order = 1),size = guide_legend(order = 2))+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust= 1,size = 14),
axis.text.y = element_text(size = 12),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
dt$annotation <- factor(dt$annotation, levels = names(gene.functions) %>% rev)
plt2h <- dt %>% ggplot(aes(y = annotation, x = hgnc_symbol, fill = log2FoldChange, size = metric))+geom_point(shape=21)+theme_bw()+expression_colors(1.5,0.75)+theme_bw()+
scale_size_continuous(limits = global_size_range, range = c(1,10))+
labs(x ='',y = '', size = 'statistical significance',title = my.title)+
guides(fill = guide_colourbar(order = 1),size = guide_legend(order = 2))+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust= 1,size = 12),
axis.text.y = element_text(size = 14),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
```
##### Vertical heatmap & dot plot
```{r fig.dim=c(1420/96,920/96), fig.dpi=300}
plt1v + plt2v + plot_layout(widths = c(4,11))
```
##### Horizontal heatmap & dot plot
```{r fig.dim=c(1300/96,700/96), fig.dpi=300}
plt1h / plt2h + plot_layout(heights = c(4,11))
```
#### {-}
# Dynamic genes
Principal component analysis shows that PC1 neatly separates samples based on the developmental time point. Therefore, genes whose expression highly correlates with PC1 are those which significantly change over the developmental timeline. To identify those dynamic genes, a correlation analysis between gene expression and PC1 was performed and genes with the most significant correlations were tested for functional enrichment using *EnrichR*.
```{r include=FALSE}
norm.cts <- assay(ntd)
identical(rownames(t(norm.cts)), rownames(pcaData))
# Calculate correlation between gene expression and first principal component (PC1)
# cor.res <- cor(pcaData$PC1,t(norm.cts))
# hist(cor.res[1,])
# cor.res[1,] %>% sort() %>% tail(50)
# sum(abs(cor.res[1,]) > 0.9 )
cor.res <- apply(t(norm.cts),2,cor.test, y= pcaData$PC1)
cor.coef <- sapply(cor.res, function(x) unname(x$estimate))
pval <- sapply(cor.res, function(x) x$p.value)
padj <- p.adjust(pval,method = 'fdr')
#data.frame(cor.coef = cor.coef, pval = pval) %>% ggplot(aes(x = cor.coef, y = pval))+geom_point()
sum(padj < 5/100)
g.id <- names(padj[padj < 5/100])
g.sym <- canonical_mapping %>% filter(ensembl_gene_id %in% g.id) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname %>% unique
if (websiteLive) {
enriched <- enrichr(g.sym, dbs)
}
enriched.df <- list_rbind(enriched[sapply(enriched, nrow) != 0], names_to = 'database')
top.g.id <- sort(pval, decreasing = F) %>% head(1000) %>% names()
top.g.sym <- canonical_mapping %>% filter(ensembl_gene_id %in% top.g.id) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname %>% unique
length(top.g.id);length(unique(top.g.sym))
top.g.pval <- pval[names(pval) %in% top.g.id]
top.g.coef <- cor.coef[names(cor.coef) %in% top.g.id]
identical(rownames(as.data.frame(top.g.pval)), rownames(as.data.frame(top.g.coef)))
top.g.df <- inner_join(as.data.frame(top.g.coef) %>% rownames_to_column('ensembl_gene_id') ,
as.data.frame(top.g.pval) %>% rownames_to_column('ensembl_gene_id'), by = 'ensembl_gene_id')
colnames(top.g.df)[c(2,3)] <- c('cor.coef','pval')
top.g.df <- left_join(top.g.df, canonical_mapping %>% dplyr::select(ensembl_gene_id,hgnc_symbol), by = 'ensembl_gene_id')
```
The top dynamic genes are shown in the following table together with their correlation coefficient with PC1, the P-value for the correlation and the FDR-adjusted P-value. Genes are ranked based on the statistical significance of their correlation. A negative correlation indicates that the gene decreases in expression and a positive correlation indicates that it increases in expression with developmental time.
```{r}
#colnames(top.g.df)[c(1,4,2,3)]
top.g.df[,c(1,4,2,3)] %>% mutate(padj = p.adjust(pval, method= 'fdr')) %>% arrange(pval) %>%
datatable(options = list(pageLength= 10))
```
The results of functional enrichment analysis of those genes are shown in the following table. Gene sets are arranged by P-value.
```{r}
enriched.df %>% arrange(P.value) %>% head(1000) %>%
datatable(options = list(pageLength= 10))
```
The developmental time profiles of the expression of the top 50 genes are shown here.
```{r include=FALSE}
top.g.id <- sort(pval, decreasing = F) %>% head(50) %>% names()
top.g.sym <- canonical_mapping %>% filter(ensembl_gene_id %in% top.g.id) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname %>% unique
length(top.g.id);length(unique(top.g.sym))
top.g.pval <- pval[names(pval) %in% top.g.id]
top.g.coef <- cor.coef[names(cor.coef) %in% top.g.id]
identical(rownames(as.data.frame(top.g.pval)), rownames(as.data.frame(top.g.coef)))
top.g.df <- inner_join(as.data.frame(top.g.coef) %>% rownames_to_column('ensembl_gene_id') ,
as.data.frame(top.g.pval) %>% rownames_to_column('ensembl_gene_id'), by = 'ensembl_gene_id')
colnames(top.g.df)[c(2,3)] <- c('cor.coef','pval')
top.g.df <- left_join(top.g.df, canonical_mapping %>% dplyr::select(ensembl_gene_id,hgnc_symbol), by = 'ensembl_gene_id')
sorted.sym <- top.g.df %>% arrange(pval) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
```
```{r fig.dim= c(1200/96,900/96), fig.dpi= 300}
df <- cbind(coldata,norm.cts[top.g.id,] %>% t)
df <- df %>% pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 7:(6+length(top.g.id))) %>%
left_join(canonical_mapping[,-3], by = 'ensembl_gene_id')
#identical(unique(df$hgnc_symbol), sorted.sym)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = sorted.sym) # or levels = unique(df$hgnc_symbol)
plt <- df %>% ggplot(aes(x = timepoint, y = exp))+theme_bw()+
geom_jitter(color = 'gray', width = 0.2, height = 0)+
geom_point(aes(x = timepoint, y= mean.exp),size = 1, data = df %>% group_by(hgnc_symbol,timepoint) %>% summarise(n_samples = n(), mean.exp = mean(exp)))+
geom_line(aes(x = timepoint, y= mean.exp, group = 1), data = df %>% group_by(hgnc_symbol,timepoint) %>% summarise(n_samples = n(), mean.exp = mean(exp)))+
facet_wrap(~ hgnc_symbol, scales = 'free')+
labs(x = '', y = 'Normalized expression')
plt
```
## Dynamic genes significantly upregulated by EPS1b
```{r fig.dim=c(850/96, 500/96), fig.dpi=300}
g = intersect(up1b,g.id)
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% rownames_to_column('sample') %>%
filter(group != 'EPS2') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol))
#df$group <- factor(df$group, levels = c('Control','EPS1a','EPS1b','EPS2'))
#df$group2 <- case_match(df$group, c('EPS1a','EPS1b') ~ 'EPS', .default = df$group)
df$group2 <- ifelse(df$group == 'Control','Control','EPS stable')
#df$group3 <- paste(df$group2,df$plate_no, sep = '_')
df %>% ggplot(aes(x = timepoint, y = exp, color = group2, shape = plate_no))+
#geom_col(width = 0.5,position = position_dodge())+
geom_point()+
geom_line(data = df %>% filter(plate_no != 'P1575' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(plate_no != 'P1571' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_03','RNA_05')), aes(x = timepoint, y = exp, group = 1), color = 'tomato')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_08','RNA_10')), aes(x = timepoint, y = exp, group = 1), color = 'tomato')+
theme_bw()+facet_wrap(~hgnc_symbol,scales = 'free_y')+
#scale_color_manual(values = c('darkolivegreen4', 'tomato','dodgerblue','violet'))+
scale_color_manual(values = c('darkolivegreen4', 'tomato'))+
labs(x = '', y = 'Normalized expression')
```
### Heatmap & Dot plot
```{r include=FALSE}
my.title = "Chronic EPS, stable vs Control comparison"
# perform enrichment analysis for significant genes
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
my.genes <- unique(mapping.df$hgnc_symbol)
# print numbers of ensembl IDs and corresponding gene symbols
#length(g);length(my.genes)
enriched <- enrichr(my.genes, dbs)
enriched.df <- list_rbind(enriched[sapply(enriched, nrow) != 0], names_to = 'database')
# Biological functions in Neurons & Muscles
gene.functions = list(`Neuronal functions` = extract_genes(enriched.df,'neuron'),
`Axon functions` = extract_genes(enriched.df,'axon'),
`Dendrite functions` = extract_genes(enriched.df,'dendrite'),
`Myelination`= extract_genes(enriched.df,'myelin'),
`Synaptic functions` = extract_genes(enriched.df,'synaps|synaptic transmission'),
`Myogenesis/Myotube formation`= extract_genes(enriched.df,'myogenesis|myotube|myoblast fusion'),
`Muscle contraction`= extract_genes(enriched.df,'muscle contract|muscular contract'),
`Angiogenesis`= extract_genes(enriched.df,'angiogenesis|vascular endothelial growth|blood vessel remodeling|blood vessel morphogenesis|blood vessel sprouting'),
`Endothelial cell functions`= extract_genes(enriched.df,'\\b(endothelial cell|endothelial cells)\\b'),
`Apoptosis/Cell death` = extract_genes(enriched.df,'apopto|cell death'),
`Mitochondrial functions` = extract_genes(enriched.df,'mitochondr|atp synthesis|oxidative phosphory'),
`Oxidative stress` = extract_genes(enriched.df,'oxidative stress')
)
# create gene-function data frame
gene.functions.df <- lapply(gene.functions, function(x) data.frame(gene = x)) %>% list_rbind(names_to = 'annotation')
my.genes <- unique(gene.functions.df$gene)
orf.genes <- grep(x = my.genes,pattern = 'orf',ignore.case = T,value = T)
my.genes[my.genes %in% orf.genes] <- str_to_title(orf.genes)
gene.functions.df$gene[gene.functions.df$gene %in% orf.genes] <- str_to_title(orf.genes)
# Heatmap
g.sub <- canonical_mapping %>% filter(hgnc_symbol %in% unique(gene.functions.df$gene)) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
df <- cbind(coldata,norm.cts[g.sub,] %>% t %>% scale)
df <- df %>% rownames_to_column('sample') %>% filter(timepoint == 'd60' & group != 'EPS2') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g.sub))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
df$group_plate <- paste(df$group,df$plate_no, sep = '_')
df$group_plate <- factor(df$group_plate, levels = c("Control_P1571","Control_P1575","EPS1b_P1571","EPS1b_P1575"))
df$group_replicate <- paste(df$group_desc,df$replicate, sep = ' ')
df$group_replicate <- factor(df$group_replicate, levels = c("Control Rep1","Control Rep2","Chronic EPS, stable Rep1","Chronic EPS, stable Rep2"))
#table(df$group_plate,df$group_replicate)
plt1h <- df %>% ggplot(aes(y= group_replicate, x= hgnc_symbol, fill= exp)) + # instead of y = group_plate
geom_tile() +
theme_bw()+
theme(panel.grid.major = element_blank())+
scale_x_discrete(position='bottom')+
#scale_fill_gradient2(low = 'dodgerblue', mid = 'white', high = 'tomato', midpoint = 0, limits = c(-1.2,1.2))+
expression_colors(1.5,0.75)+
labs(x= '',y= '',fill= 'Scaled expression', title = my.title)+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust = 1,size = 12),
axis.text.y = element_text(size = 14),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
# Gene-Function dot plot
dt <- res1b %>% filter(hgnc_symbol %in% my.genes) %>% dplyr::select(hgnc_symbol,log2FoldChange,pvalue) %>% mutate(metric = -log10(pvalue))
# unique(gene.functions.df$gene)
# unique(dt$hgnc_symbol)
colnames(gene.functions.df)[2] <- 'hgnc_symbol'
dt <- left_join(gene.functions.df, dt, by = 'hgnc_symbol')
summary(dt$log2FoldChange)
dt$annotation <- factor(dt$annotation, levels = names(gene.functions) %>% rev)
plt2h <- dt %>% ggplot(aes(y = annotation, x = hgnc_symbol, fill = log2FoldChange, size = metric))+geom_point(shape=21)+theme_bw()+expression_colors(1.5,0.75)+theme_bw()+
labs(x ='',y = '', size = 'statistical significance',title = my.title)+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust= 1,size = 12),
axis.text.y = element_text(size = 14),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
```
```{r fig.dim=c(800/96,710/96), fig.dpi=300}
plt1h / plt2h + plot_layout(heights = c(4,12))
```
## Dynamic genes significantly upregulated by EPS2
```{r fig.dim=c(850/96, 550/96), fig.dpi=300}
g = intersect(up2,g.id)
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
#df <- cbind(coldata,norm.cts[rownames(norm.cts) %in% g,] %>% t)
df <- cbind(coldata,norm.cts[g,] %>% t)
df <- df %>% rownames_to_column('sample') %>%
filter(group != 'EPS1b') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
# ism <- res1b %>% filter(padj < 0.05) %>% arrange(pvalue) %>% dplyr::select(hgnc_symbol) %>% unlist %>% unname
# identical(unique(df$hgnc_symbol), ism)
df$hgnc_symbol <- factor(df$hgnc_symbol, levels = unique(df$hgnc_symbol))
#df$group <- factor(df$group, levels = c('Control','EPS1a','EPS1b','EPS2'))
df$group2 <- ifelse(df$group == 'Control','Control','EPS stable/increasing')
#df$group3 <- paste(df$group2,df$plate_no, sep = '_')
df %>% ggplot(aes(x = timepoint, y = exp, color = group2, shape = plate_no))+
#geom_col(width = 0.5,position = position_dodge())+
geom_point()+
geom_line(data = df %>% filter(plate_no != 'P1575' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(plate_no != 'P1571' & group == 'Control'), aes(x = timepoint, y = exp, group = 1), color = 'darkolivegreen4')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_03','RNA_06')) , aes(x = timepoint, y = exp, group = 1), color = 'dodgerblue')+
geom_line(data = df %>% filter(sample %in% c('RNA_01','RNA_08','RNA_11')), aes(x = timepoint, y = exp, group = 1), color = 'dodgerblue')+
theme_bw()+facet_wrap(~hgnc_symbol,scales = 'free_y')+
#scale_color_manual(values = c('darkolivegreen4', 'tomato','dodgerblue','violet'))+
scale_color_manual(values = c('darkolivegreen4', 'dodgerblue'))+
labs(x = '', y = 'Normalized expression')
```
### Heatmap & Dot plot
```{r include=FALSE}
my.title = "Chronic EPS, increasing vs Control comparison"
# perform enrichment analysis for significant genes
mapping.df <- canonical_mapping %>% filter(ensembl_gene_id %in% g)
my.genes <- unique(mapping.df$hgnc_symbol)
# print numbers of ensembl IDs and corresponding gene symbols
#length(g);length(my.genes)
enriched <- enrichr(my.genes, dbs)
enriched.df <- list_rbind(enriched[sapply(enriched, nrow) != 0], names_to = 'database')
# Biological functions in Neurons & Muscles
gene.functions = list(`Neuronal functions` = extract_genes(enriched.df,'neuron'),
`Axon functions` = extract_genes(enriched.df,'axon'),
`Dendrite functions` = extract_genes(enriched.df,'dendrite'),
`Myelination`= extract_genes(enriched.df,'myelin'),
`Synaptic functions` = extract_genes(enriched.df,'synaps|synaptic transmission'),
`Myogenesis/Myotube formation`= extract_genes(enriched.df,'myogenesis|myotube|myoblast fusion'),
`Muscle contraction`= extract_genes(enriched.df,'muscle contract|muscular contract'),
`Angiogenesis`= extract_genes(enriched.df,'angiogenesis|vascular endothelial growth|blood vessel remodeling|blood vessel morphogenesis|blood vessel sprouting'),
`Endothelial cell functions`= extract_genes(enriched.df,'\\b(endothelial cell|endothelial cells)\\b'),
`Apoptosis/Cell death` = extract_genes(enriched.df,'apopto|cell death'),
`Mitochondrial functions` = extract_genes(enriched.df,'mitochondr|atp synthesis|oxidative phosphory'),
`Oxidative stress` = extract_genes(enriched.df,'oxidative stress')
)
# create gene-function data frame
gene.functions.df <- lapply(gene.functions, function(x) data.frame(gene = x)) %>% list_rbind(names_to = 'annotation')
my.genes <- unique(gene.functions.df$gene)
orf.genes <- grep(x = my.genes,pattern = 'orf',ignore.case = T,value = T)
my.genes[my.genes %in% orf.genes] <- str_to_title(orf.genes)
gene.functions.df$gene[gene.functions.df$gene %in% orf.genes] <- str_to_title(orf.genes)
# Heatmap
g.sub <- canonical_mapping %>% filter(hgnc_symbol %in% unique(gene.functions.df$gene)) %>% dplyr::select(ensembl_gene_id) %>% unlist %>% unname
df <- cbind(coldata,norm.cts[g.sub,] %>% t %>% scale)
df <- df %>% rownames_to_column('sample') %>% filter(timepoint == 'd60' & group != 'EPS1b') %>%
pivot_longer(names_to = "ensembl_gene_id",values_to = "exp",cols = 8:(7+length(g.sub))) %>%
left_join(mapping.df[,-3], by = 'ensembl_gene_id')
df$group_plate <- paste(df$group,df$plate_no, sep = '_')
df$group_plate <- factor(df$group_plate, levels = c("Control_P1571","Control_P1575","EPS2_P1571","EPS2_P1575"))
df$group_replicate <- paste(df$group_desc,df$replicate, sep = ' ')
df$group_replicate <- factor(df$group_replicate, levels = c("Control Rep1","Control Rep2","Chronic EPS, increasing Rep1","Chronic EPS, increasing Rep2"))
#table(df$group_plate,df$group_replicate)
plt1h <- df %>% ggplot(aes(y= group_replicate, x= hgnc_symbol, fill= exp)) + # instead of y = group_plate
geom_tile() +
theme_bw()+
theme(panel.grid.major = element_blank())+
scale_x_discrete(position='bottom')+
#scale_fill_gradient2(low = 'dodgerblue', mid = 'white', high = 'tomato', midpoint = 0, limits = c(-1.2,1.2))+
expression_colors(1.5,0.75)+
labs(x= '',y= '',fill= 'Scaled expression', title = my.title)+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust = 1,size = 12),
axis.text.y = element_text(size = 14),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
# Gene-Function dot plot
dt <- res1b %>% filter(hgnc_symbol %in% my.genes) %>% dplyr::select(hgnc_symbol,log2FoldChange,pvalue) %>% mutate(metric = -log10(pvalue))
# unique(gene.functions.df$gene)
# unique(dt$hgnc_symbol)
colnames(gene.functions.df)[2] <- 'hgnc_symbol'
dt <- left_join(gene.functions.df, dt, by = 'hgnc_symbol')
summary(dt$log2FoldChange)
dt$annotation <- factor(dt$annotation, levels = names(gene.functions) %>% rev)
plt2h <- dt %>% ggplot(aes(y = annotation, x = hgnc_symbol, fill = log2FoldChange, size = metric))+geom_point(shape=21)+theme_bw()+expression_colors(1.5,0.75)+theme_bw()+
labs(x ='',y = '', size = 'statistical significance',title = my.title)+
theme(plot.title = element_text(hjust = 0.5, size = 16),
axis.text.x = element_text(angle = 30, hjust= 1,size = 12),
axis.text.y = element_text(size = 14),
legend.position = 'right',legend.title.position = 'left', legend.title = element_text(angle = 90))
```
```{r fig.dim=c(800/96,710/96), fig.dpi=300}
plt1h / plt2h + plot_layout(heights = c(4,12))
```
# R session information
```{r}
utils::sessionInfo()
```

EPS_NMOs_analysis_v2.Rmd at commit 8eb2cca, under MIT · at the source

Overview

Authors: Chrysanthi‐Maria Moysidou¹, Inês Afonso Martins¹, Ismail Amr El‐Shimy¹, Iacopo Bicci¹, Donatella Cea², Christina Bukas², Isra Mekki², Mara‐Camelia Rusu¹, Aylin Nebol¹, Ines Lahmann¹, Marie Piraud², Enrico Klotzsch^3,4,5, Mina Gouti^1,6

ORCID iDs: Chrysanthi‐Maria Moysidou, Inês Afonso Martins, Ismail Amr El‐Shimy, Iacopo Bicci, Donatella Cea, Enrico Klotzsch, Mina Gouti

Max Delbrück Center for Molecular Medicine (MDC) Berlin Germany
Helmholtz AI, Computational Health Center (CHC) Helmholtz Munich Neuherberg Germany
Berlin Institute of Health Center for Regenerative Therapies (BCRT) Berlin Germany
Berlin Institute of Health Experimental and Clinical Research Center (ECRC) Berlin Germany
Charité Universitätsmedizin Berlin Germany
Department of Pediatric Neurology Charité Universitätsmedizin Berlin Berlin Germany

Institutions: Max Delbrück Center (Germany); Helmholtz Munich (Germany); Berlin-Brandenburger Centrum für Regenerative Therapien (Germany); Charité - Universitätsmedizin Berlin (Germany)

Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany), volume 13, issue 50, article e22762

Dates: received 13 November 2025; accepted 5 June 2026; published online 22 June 2026; in print September 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1002/advs.202522762 · PMID 42325128 · PMCID PMC13337066 · OpenAlex W7165545725

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), developmental (subfield)

Keywords: bioengineering, biophysical cues, electrical stimulation, EPS‐NMOs, neuromuscular organoids, neuromuscular systems, organoid maturation

MeSH: Electric Stimulation*, Muscle, Skeletal*, Neuromuscular Junction*, Organoids*, Pluripotent Stem Cells*, Cell Differentiation, Humans (* major topic)

Topic: Planarian Biology and Electrostimulation (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: European Molecular Biology Organization (EMBO); Einstein Stiftung Berlin (EZ‐2020‐597‐2); Deutsche Forschungsgemeinschaft (KL 3278/2‐1); Helmholtz Association; HORIZON EUROPE European Research Council (101002689); HORIZON EUROPE Marie Sklodowska-Curie Actions (101149182)

Citations: cited by 1 paper (Europe PMC); 93 references in the paper

Abstract

Organoids derived from human pluripotent stem cells (hPSCs) are emerging as powerful models for studying development and disease. Despite their physiological relevance, the predictive power of organoids remains limited by the immature state of the constituent cells, posing a major challenge for mechanistic studies of adult physiology and late‐onset disorders. Here, we establish a strategy for enhancing the maturation status of human neuromuscular organoids (NMOs) through chronic Electrical Pulse Stimulation (EPS). We demonstrate that low‐frequency EPS, applied during the early stages of NMO development and maintained over several weeks, enhances neuromuscular maturation and functional output. Independent of stimulation waveform dynamics, EPS‐trained NMOs (EPS‐NMOs) display stronger and more frequent spontaneous contractions that persist long after stimulation has ceased. Quantitative imaging and transcriptomic analyses reveal a robust improvement in EPS‐NMO skeletal muscle and neural tissue morphology, coordinated regulation of lineage‐specific biomarkers, and upregulation of gene programs associated with neuromuscular maturation. Mechanobiological measurements further demonstrate increased tissue stiffness and faster relaxation dynamics in EPS‐NMOs, consistent with enhanced excitation‐contraction coupling (ECC) and force generation. Collectively, these findings establish EPS as a powerful, non‐invasive, and on‐demand modality for promoting the morphological and functional maturation of complex organoid systems.

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

Repositories

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

Gouti-Lab/EPS-Project

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 8eb2ccab44fac1b8c61fbd493b55a05d963c75bb, 23 March 2026

Languages: Jupyter (8), R (1), Shell (1)

Size: 15 files, 10 scripts

Software Heritage: not archived

Found in: “Data Availability Statement”

Holds: README, license file, 9 notebooks

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

Tools: Matplotlib (8 files), NumPy (8 files), pandas (8 files), scikit-image (8 files), imageio (7 files), h5py (4 files), tifffile (4 files), DESeq2 (1 file), patchwork (1 file), pheatmap (1 file), SAMtools (1 file), SciPy (1 file), Subread (featureCounts) (1 file), tidyverse (1 file)

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

27 September 2026: the link answers

12 files

Image-analysis/Apoptosis levels.ipynb, Jupyter, 128 lines
Image-analysis/BTX & Innervation analysis.ipynb, Jupyter, 246 lines, 1 match
Image-analysis/ChAT area analysis.ipynb, Jupyter, 280 lines, 1 match
Image-analysis/Col4A1 or Col5A3 area analysis.ipynb, Jupyter, 319 lines, 2 matches
Image-analysis/GFAP analysis.ipynb, Jupyter, 303 lines
Image-analysis/MYOD analysis.ipynb, Jupyter, 129 lines, 1 match
Image-analysis/PAX7-Ki67 analysis.ipynb, Jupyter, 332 lines, 2 matches
Image-analysis/SOX1-Ki67 analysis.ipynb, Jupyter, 154 lines, 1 match
Transcriptomics-analysis/EPS_NMOs_analysis_v2.Rmd, R, 1,532 lines, 4 matches
Transcriptomics-analysis/rna_seq_pipeline_v2.sh, Shell, 68 lines, 1 match
LICENSE, License, 21 lines
README.md, Text, 12 lines

At the source: github.com/Gouti-Lab/EPS-Project

HelmholtzAI-Consultants-Munich/NMOs-Contraction

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 58bf0ae2227d6ef5a657d2191f2ed4ae94baf63d, 12 June 2026

Languages: Python (5), Jupyter (3)

Size: 16 files, 8 scripts

Software Heritage: not archived

Found in: “Data Availability Statement”

Holds: README, environment (requirements.txt), 3 notebooks

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

Tools: NumPy (6 files), pandas (5 files), Matplotlib (4 files), scikit-image (3 files), Plotly (2 files), Cellpose (1 file), Pillow (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file)

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

27 September 2026: the link answers

9 files

src/1_time_series_extraction.ipynb, Jupyter, 447 lines
src/2_time_series_analysis.ipynb, Jupyter, 969 lines, 1 match
src/3_border_and_area_computation.ipynb, Jupyter, 278 lines, 1 match
src/utils/feature_extraction.py, Python, 74 lines
src/utils/imaging_viewers.py, Python, 152 lines
src/utils/plot_boxplot.py, Python, 18 lines
src/utils/pre_processing.py, Python, 49 lines, 1 match
src/utils/time_series_extraction.py, Python, 216 lines, 2 matches
README.md, Text, 111 lines

At the source: github.com/HelmholtzAI-Consultants-Munich/NMOs-Contraction

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:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
18 scripts, each with its path and the digest of its content;
18 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

No dataset and no data link were found in the paper.

Data Availability Statement

The data that support the findings of this study are available from the corresponding authors upon reasonable request. All analyses of sequencing data were executed on a SLURM‐managed HPC cluster (Max Cluster) at MDC‐Berlin using GNU parallel to parallelize QC, alignment, and counting across samples. Software versions and links to reference files (genome build and GTF) are provided in Table S3 and command‐line arguments and analysis scripts are available on GitHub (https://github.com/Gouti‐Lab/EPS‐Project (https://github.com/Gouti-Lab/EPS-Project)). Image analysis was carried out using FiJi, iLastik and python scripts. Code is available on GitHub (https://github.com/Gouti‐Lab/EPS‐Project (https://github.com/Gouti-Lab/EPS-Project)). The code used for size analysis and spontaneous contraction characterization in Python is available at https://github.com/HelmholtzAI‐Consultants‐Munich/NMOs‐Contraction (https://github.com/HelmholtzAI-Consultants-Munich/NMOs-Contraction). The programme for quantification of mechanobiological properties is available at https://github.com/enricoklotzsch/organoid_stress_relax/.

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, 13 authors, 7 keywords, 7 MeSH terms, 6 funders, 89 references.

Cite

This paper

Moysidou, C., Afonso Martins, I., El‐Shimy, I. A., Bicci, I., Cea, D., Bukas, C., Mekki, I., Rusu, M., Nebol, A., Lahmann, I., Piraud, M., Klotzsch, E., & Gouti, M. (2026). Enhancing Maturation of Human Neuromuscular Organoids via Electrical Stimulation. Advanced science (Weinheim, Baden-Wurttemberg, Germany), 13(50), e22762. https://doi.org/10.1002/advs.202522762

BibTeX

@article{moysidou2026enhancing, author = {Moysidou, Chrysanthi‐Maria and Afonso Martins, Inês and El‐Shimy, Ismail Amr and Bicci, Iacopo and Cea, Donatella and Bukas, Christina and Mekki, Isra and Rusu, Mara‐Camelia and Nebol, Aylin and Lahmann, Ines and Piraud, Marie and Klotzsch, Enrico and Gouti, Mina}, title = {{Enhancing Maturation of Human Neuromuscular Organoids via Electrical Stimulation}}, journal = {Advanced science (Weinheim, Baden-Wurttemberg, Germany)}, year = {2026}, month = jun, volume = {13}, number = {50}, pages = {e22762}, publisher = {Wiley}, issn = {2198-3844}, doi = {10.1002/advs.202522762}, url = {https://doi.org/10.1002/advs.202522762}, pmid = {42325128}, pmcid = {PMC13337066} }

RIS

TY - JOUR AU - Moysidou, Chrysanthi‐Maria AU - Afonso Martins, Inês AU - El‐Shimy, Ismail Amr AU - Bicci, Iacopo AU - Cea, Donatella AU - Bukas, Christina AU - Mekki, Isra AU - Rusu, Mara‐Camelia AU - Nebol, Aylin AU - Lahmann, Ines AU - Piraud, Marie AU - Klotzsch, Enrico AU - Gouti, Mina TI - Enhancing Maturation of Human Neuromuscular Organoids via Electrical Stimulation T2 - Advanced science (Weinheim, Baden-Wurttemberg, Germany) J2 - Adv Sci (Weinh) PY - 2026 DA - 2026/06/22 VL - 13 IS - 50 SP - e22762 SN - 2198-3844 PB - Wiley DO - 10.1002/advs.202522762 UR - https://doi.org/10.1002/advs.202522762 LA - en ER -

CSL-JSON

{ "id": "10.1002/advs.202522762", "type": "article-journal", "title": "Enhancing Maturation of Human Neuromuscular Organoids via Electrical Stimulation", "container-title": "Advanced science (Weinheim, Baden-Wurttemberg, Germany)", "author": [ { "family": "Moysidou", "given": "Chrysanthi‐Maria" }, { "family": "Afonso Martins", "given": "Inês" }, { "family": "El‐Shimy", "given": "Ismail Amr" }, { "family": "Bicci", "given": "Iacopo" }, { "family": "Cea", "given": "Donatella" }, { "family": "Bukas", "given": "Christina" }, { "family": "Mekki", "given": "Isra" }, { "family": "Rusu", "given": "Mara‐Camelia" }, { "family": "Nebol", "given": "Aylin" }, { "family": "Lahmann", "given": "Ines" }, { "family": "Piraud", "given": "Marie" }, { "family": "Klotzsch", "given": "Enrico" }, { "family": "Gouti", "given": "Mina" } ], "container-title-short": "Adv Sci (Weinh)", "volume": "13", "issue": "50", "page": "e22762", "DOI": "10.1002/advs.202522762", "PMID": "42325128", "PMCID": "PMC13337066", "ISSN": "2198-3844", "publisher": "Wiley", "URL": "https://doi.org/10.1002/advs.202522762", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 22 ] ] } }

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

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Subread (featureCounts), SAMtools, DESeq2, 12 other tools
[2] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: imageio, tifffile, DESeq2, 12 other tools
[3] doi:10.1371/journal.pcbi.1014571 [code]: SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: Cellpose, imageio, tifffile, 8 other tools, 2 references
[4] doi:10.1038/s41467-026-71458-0 [code]: Early differential impact of MeCP2 mutations on functional networks in Rett syndrome patient-derived human cortical organoids.
Journal: Nature communications
In common: Cellpose, tifffile, pheatmap, 7 other tools, 2 references
[5] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: Subread (featureCounts), SAMtools, DESeq2, 9 other tools
[6] doi:10.1093/neuonc/noag128 [code]: Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response.
Journal: Neuro-oncology
In common: Cellpose, imageio, tifffile, 8 other tools
[7] doi:10.1073/pnas.2609132123 [code]: 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
In common: Cellpose, tifffile, pheatmap, 9 other tools
[8] doi:10.1038/s41467-026-71803-3 [code]: Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: SAMtools, DESeq2, pheatmap, 9 other tools
[9] doi:10.1186/s11689-026-09713-0 [code]: DRP1 mutations associated with EMPF1 encephalopathy perturb the transcriptional profile and maturation of cortical neurons.
Journal: Journal of neurodevelopmental disorders
In common: tifffile, DESeq2, pheatmap, 9 other tools
[10] doi:10.1038/s42003-026-10063-9 [code]: Conserved Kir channel mechanisms governing intrinsic excitability in human and rodent parvalbumin neurons.
Journal: Communications biology
In common: Cellpose, tifffile, Plotly, 7 other tools, 1 reference

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.