Spaziotemporal Cortical Transcriptome: Modelling Trajectories Across the Lifespan
An Attempt, a Failure, and the Lessons Learned
Luca Sclisizzo
2026-07-18
1 Field Standards and Motivation
1.1 Field Standards
To date, the Brainspan / PsychENCODE (Li et al. 2018) database is widely used in Behavioral Genetic literature. It is a high throughput database that contains postmortem brain samples collected from the 7th Post Conception Week (PCW) to adulthood. These brains are classified in 9 developmental windows and each brain is sampled from 16 different Regions Of Interest (ROI), leading to a hierarchical structure of the data.
In the Behavioral Genetic field, there are two standard ways to analyse this database:
- Differential Expression (DE) tests between different developmental windows using gene-wise models
- Non-parametric transcriptomic trajectories across lifespan
In this project I will focus on the second way to tackle this problem, that is quite unique since RNA-seq expression trajectories are usually studied within the cell life, but they are rarely evaluated across the organism lifespan. Moreover, the DE tradition (the most widely used RNA-seq analyses) does not focus on lifelong expression trajectories, leaving it an under explored research area, with new and open analytical challenges.
To date, the pipeline for the lifetime trajectory analysis consists on 3 steps:
- Selecting a list of genes of interest. Behavioral geneticists usually take the genes that are significant in a Genome Wide Association Study (GWAS).
- Get the mean across gene expression for each biological sample - that is, a ROI in a specific developmental stage. Usually, the selected genes are taken from GWAS based on European ancestry, thus only ROIs from European samples are included, limiting the sample size of the dataset to (\(\mathrm{n}\sim22\)).
- Smooth the ROI means across stages, usually using a local regression (LOESS).
1.2 The Grey Area
RNA-seq data present several intrinsic analytical challenges. In particular, normalization is an important step in the transcriptome analyses. It is well known that raw RNA counts contain both biological and technical source of variations: chips, batch, sample composition and preparation do matter in the RNA raw counts. Therefore, it is crucial to be able to separate the two sources of variability through the normalization process.
During the current “genomic era”, several normalization methods and packages have been developed. However, there is still no consensus on the best normalization method. Moreover, due to the uniqueness of this dataset, no study has yet explored the influence of the normalization methods in the trajectory shapes. In addition, there is no standard strategy across studies: scholars are using log-raw counts (Stauffer et al. 2023), upperquartile (Gandal et al. 2018), RPKM (Grotzinger et al. 2026; Dear et al. 2024) and more. This heterogeneity hinders the development of a coherent analytical framework across studies and limits the reader’s ability to directly compare findings across studies and phenotypes.
2 The Project
2.1 General Objectives
The project aimed to answer two main questions:
- Is it possible to model transcriptome trajectories with a parametric model within a bayesian framework?
- Does the normalization procedure influence the shape of the inferred trajectories?
2.2 Schizophrenia as Case Study
This study integrates data from the BrainSpan and PsychENCODE (Li et al. 2018) resources with the schizophrenia GWAS from the Psychiatric Genomics Consortium (PGC) (Trubetskoy et al. 2022). The GWAS summary statistics are firstly analyzed using hMAGMA (Leeuw et al. 2016), which maps Single Nucleotide Polymorphisms (SNPs) to genes. This mapping incorporates cortical annotations and allows the integration of:
- genomic position of SNPs
- functional effects through expression quantitative trait loci (eQTLs)
- 3D genome architecture via chromatin interaction (Hi-C) data
This analytical tool enables the identification of genes involved in the polygenic architecture of schizophrenia. The hMAGMA results were previously performed and the list of the Bonferroni significant genes is provided in SCZ_genes.xlsx file. The hMAGMA analysis were submitted with default setting of the algorithm:
FILE="SCZ_sumstats.tsv" # schizophrenia GWAS summary statistics
annot="./Adult_brain.genes.annot" # Adult brain annotations
echo ">> H-MAGMA analysis..."
/Magma/magma \
--bfile /g1000_eur \ # EUR LD files from 1000Genome reference panel
--pval $FILE snp-id=SNP pval=PVAL ncol=NEFF \
--gene-annot $annot \
--out ./hMAGMA_results2.3 Data Sources
The datasets used in this study are publicly available:
- Psychiatric Genomic Consortium figshare for schizophrenia summary statistic
- PsychENCODE resource including sample metadata and raw counts
3 Exploratory Analyses
I started my analyses comparing different normalization methods using the standard non parametric approach. In particular, I took the Brainspan raw counts and used the package edgeR to normalize the data. I mostly constrained the normalization methods among the between-sample methods, as they allow comparisons between biological samples (Evans et al. 2018). Thus, the selected algorithm were: Upperquartile, RLE and TMM. However, since RPKM normalization is widely used for these analyses (Grotzinger et al. 2026; Dear et al. 2024), I decided to include it. It should be noted that this is a within-sample normalization and should not be used to compare different biological samples (Harvard Chan Bioinformatics Core 2026). With these normalizing strategies, I computed the trajectories using LOESS as a smoothing process (see Field Standards for pipeline details).
Figure 1: The figure plots the SCZ trajectories with different normalization methods. RLE is the standard normalization in DESeq2 package, TMM is a native edgeR method, UpperQuartile is a standard between-sample normalization and RPKM is a normalization that does not account for library composition (a within-sample normalization). All the trajectories are squeezed between (-1,1) to better compare the shapes.
The differences between normalizations are non trival, further motivating the aims of this project. Compared with the LOESS strategy, a bayesian model has the advantage to fully capture the uncertainity of the parameters, enabling a more through assessment of the impact of the normalization step on the transcriptomic trajectories.
3.1 PCA on Samples
The PCA exploratory step aimed to better characterize expression variability across developmental stages and to assess the similarity among ROIs within each stage of expression (called window in this context). To achieve this, a PCA was performed after a logarithmic transformation to account for the skeweness of the data. This analysis can help evaluate the extent to which different normalization methods capture expression variability, as well as the importance of hierarchical structure. In particular, it could addresses whether a complete pooling approach is appropriate or whether the hierarchical structure of the data accounts for a substantial proportion of the observed variability.
The PC1 and PC2 were very similar between the normalization methods. In particular, the PCA coloring highlights a time-dependent pattern where ROI that belong to the same developmental stage are closer (eg. within window 1) than ROIs that belong to different life stages (eg. window 1 and 8). It is important to note that the biggest source of variability in the PC1-PC2 scatter is mostly accounted by windows differences: within the same window (scatters with the same color) the points appeared tightly packed compared to the variation across windows (colored blobs). This is important to note for the specification of the hierarchical structure (more on this in the Data Modelling section).
The PC3 and PC4 did not highlight normalization specificities: all the techniques gave similar results. Moreover, the 3rd and 4th components were not capturing any timing or spatial effect, since the samples tended to be randomly scrambled around the axes origin.
Note that also a GLM-PCA could have been appropriate, since it has shown good performance in RNA-seq analyses where both the Poisson and Negative Binomial families appear to be appropriate distributional choices, as discussed in the GLM-PCA vignette. GLM-PCA was also tested, but results were sensitive to initialization, therefore I decided to use a more standard approach. As a reader’s note, the PCA_and_exploratory_analyses.R contains also the code for the GLM-PCA.
4 Data Modelling
4.1 Data Description & Modelling
The PsychENCODE dataset presents three types of hierarchical structure. The first and more intuitive is the brain sample hierarchy: different portion of the brain (ROI) are sampled from the same individual, leading to different within-subject clusters that can be distinguished by their age. The second hierarchy is the ROI, the regional brain area: different portions of the brain could present region-specific developmental trajectory; therefore, we could create a ROI-specific hierarchy, that encompasses all the developmental stages. The third, and more subtle, nested structures are the genes: each gene could present a specific and unique expression trajectory, however, since genes are often organised in clusters with expression similarities (Montenegro 2022), we could suppose that exist a trade-off between complete and partial pooling also among genes.
For the subsequent analyses and modelling, I used the R package rstanarm (Goodrich et al. 2025) and CmdStanR (Gabry et al. 2025) to fit a hierarchical Bayesian model using a negative binomial (\(\mathrm{NegBin}\)) distribution to model RNA expression. This distributional choice is standard in transcriptomic analyses (Robinson et al. 2010; Love et al. 2014). The offset in the \(\mathrm{NegBin}\) model can be extracted using the function edgeR::getOffset() that computes: \[\text{offset}_s = \log(\text{norm.factor}_s \times \text{lib.size}_s)\] where \(s\) refers to each biological sample, \(\mathrm{norm.factor}\) adjust for composition bias of the library and \(\mathrm{lib.size}\) takes into account the sequencing depth, also known as library size.
The first step of this project was focused on evaluating a proper modelling strategy for this dataset. For this aim, I used the TMM normalization method to evaluate model building. Finding a proper model is crucial to assess whether the choice of normalization could influence trajectory estimates (see the General Objectives section).
4.2 Complete Pooling or No Pooling?
A dataset that present a hierarchical structure often require to find a trade-off between complete and no pooling for each hierarchy. Since the dataset could present three types of hierarchy (see Data Description & Modelling), I started the modelling process taking the simplest approach. Therefore, the first model accounted both the subject-dependent and ROI-dependent structures, but it assumed that all the genes were equally expressed along the development, determining a unique and homogeneous population across genes. This is usually called complete pooling at the gene level. With these modelling choices in mind, two random components were included: a random intercept for the subject and a random component for the ROI. However, due to the high dependence between ROI variability and subject variability, the model was hardly identifiable, leading to a non convergence of the MCMC algorithm. This MCMC behavior is in line with the PCA analysis that revealed a higher variability across ages compared to the variability across ROIs.
Thus, I removed the ROI random effect and fitted the model as following: \[y_{ij} \sim \mathrm{NegBin}(\mu_{ij}, \phi) \\ \mu_{ij} = \exp(\eta_{ij}) \\ \eta_{ij} = \alpha + \sum_{k=1}^{4} B_k \, \text{window}_{ijk} + \boldsymbol{\beta}^\top \mathbf{Z}_j + u_j\] The model was specified with weakly informative priors automatically implemented by the rstanarm package. In particular, the package specifies: \[B_k \sim \mathcal{N}(0, 2.5) \quad \forall k \\ \boldsymbol{\beta} \sim \mathcal{N}(\mathbf{0}, 2.5) \\ u_j \sim \mathcal{N}(\mathbf{0}, \mathbf{\Sigma})\\ \phi \sim ~ \exp(1)\] The \(\mathrm{u_j}\) term was a random intercept for the \(\mathrm{j^{th}}\) sampled subject/brain and the term \(\sum_{k=1}^{4} \beta_k B_k(window)\) is a spline with 4 degree of freedom, where each \(\mathrm{B_k}\) is a basis spline component. A spline term was fitted since I did not expect an apriori linear trend. The \(\mathbf{Z}_j\) included control covariates, in particular \(\mathrm{sex}\), \(\mathrm{Seq.site}\) and the \(\mathrm{offset}\) computed as discussed above. As mentioned above, this model did not include any gene-wise term, resulting in complete pooling among genes while accounting for the subject-specific hierarchical structure.
The model converged properly with good \(\mathrm{R_{hat}}\) and \(\mathrm{n_{eff}}\):
head(summary(model_bayes_raw_counts_TMM)[, c("Rhat", "n_eff")], n = 10)## Rhat n_eff
## (Intercept) 0.9997845 3755
## ns(Window, df = 4)1 1.0001443 3608
## ns(Window, df = 4)2 1.0009187 2831
## ns(Window, df = 4)3 1.0006084 3928
## ns(Window, df = 4)4 0.9993619 3614
## SexM 0.9996655 3256
## Sequencing.SiteYALE 0.9995540 3572
## areaDFC 0.9999103 4531
## areaIPC 0.9999911 4382
## areaITC 0.9999271 4008However, the Posterior Predictive Checks (PPC) presented some issues with the model that could not replicate both the median, the sd and the Sum_0s of the original data:
wrap_plots(list(median_PPC, sd_PPC, sum_0s)) +
plot_annotation(title = 'PPC in complete pooling across genes',
theme = theme(plot.title = element_text(hjust = 0.5)))
Figure 2: PPC of the complete pooling across genes. This model is not fitting well.
The PPC analyses suggested that the complete pooling model at the gene-level was not adequate in capturing the distribution of the data. In fact, most of the RNA-seq models available in literature are meant to be gene-wise and the distributional properties of the data had been assessed and discussed in a gene-wise perspective. Thus, a complete pooling model could violate the \(\mathrm{NegBin}\) distribution, leading to a poor fit of the model. I assumed that the violation of the distribution was greater as the gene heterogeneity increases.
To assess if the poor fit was due to the complete pooling approach, I followed the opposite direction fitting a single gene model. This approach is called no pooling as it assumes a different population for each element of the hierarchy.
Before fitting the model, I wanted to assess the quality of the gene signal and then fit different models using genes with different quality scores. To achieved this, I performed the TMM scaling calculation and then normalized the raw counts using the edgeR::rpkm(). This step accounted for sequencing depth, compositional bias of the library and the transcript length of the gene:
edgeR::rpkm(
dge, # dge object with raw counts matrix and TMM scaling factors
gene.length = scz_genes$transcript_length, # retrieved with biomaRt::useEnsembl()
normalized.lib.size = TRUE,
log = FALSE,
prior.count = 0
)With these normalized counts, I created a quality score for each gene: since I was interested in modelling the trajectory across developmental stages, I defined the score as \[\mathrm{score}_i = \frac{\mathrm{V}_{\mathrm{between}, i}}{\mathrm{V}_{\mathrm{within}, i} + \epsilon}\] where \(\mathrm{V}_{between}\) was the variance between the windows, and \(\mathrm{V}_{within}\) denotes the average within-window variance across all windows. Thus, the \(\mathrm{score}\) defined the proportion of signal over the replication noise seen as the cortex variability within each stage.
gene_score <- read.csv('gene_variability_score.csv')
head(gene_score, n= 5)## gene var_between var_within score
## 1 ENSG00000176371 323.77246 16.777380 19.29815
## 2 ENSG00000066135 83.90236 6.372103 13.16714
## 3 ENSG00000188070 670.57362 55.804506 12.01648
## 4 ENSG00000163939 99.94313 8.582936 11.64440
## 5 ENSG00000147535 57.99455 5.222241 11.10530I then fitted the first model using the ENSG00000176371 gene, as it was the gene with higher score. The model specification followed the logic explained above with the gene-wise no pooling perspective. The model had a good degree of convergence:
head(summary(model_single_gene_first_score)[, c("Rhat", "n_eff")], n= 10)## Rhat n_eff
## (Intercept) 1.0001546 2971
## ns(Window, df = 4)1 1.0011323 3068
## ns(Window, df = 4)2 1.0006675 3206
## ns(Window, df = 4)3 0.9998096 3349
## ns(Window, df = 4)4 1.0005417 2770
## SexM 0.9991130 3418
## Sequencing.SiteYALE 0.9993962 3274
## areaDFC 1.0008459 3317
## areaIPC 1.0020416 3072
## areaITC 1.0005402 3451Moreover, the PPC showed a good fit of the model, as the model was able to replicate from the posterior all the original data distributional features:
wrap_plots(list(median_single_first_gene, sd_single_first_gene, sum_0s_single_first_gene)) +
plot_annotation(title = 'PPC in no pooling across genes',
theme = theme(plot.title = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0)
))
Figure 3: Gene with the highest score: highest ratio between developmental variance and within-stage cortical variance
Therefore, I then selected another gene to check whether the \(\mathrm{score}\) could hinder the ability of the model to properly fit the data. Thus, I randomly picked the 50th gene with a ratio score of \(\sim4\) and fit the gene-wise model:
gene_50_score <- gene_score[50,]
print(gene_50_score)## gene var_between var_within score
## 50 ENSG00000120733 61.73084 14.87637 4.14959Again, the mixing was good and the \(\mathrm{n_{eff}}\) appropriate. The PPC revealed that the model was able to replicate the data with a good fit. However, compared with the previous single gene model, the \(\mathrm{y_{rep}}\) of the median were more disperse suggesting a partial decrease of model’s ability to replicate this distributional point estimate. The sd was as good as before, however, the sd scale was greater for this model than in the previous.
wrap_plots(list(median_single_50_gene, sd_single_50_gene, sum_0s_single_50_gene)) +
plot_annotation(title = 'PPC in no pooling across genes',
theme = theme(plot.title = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0)
))
Figure 4: Gene with the 50th score: 50th in terms of ratio between developmental variance and within-stage cortical variance. For this gene, the ration was ~4
With an achievement of a reliable fit, I had a look at the spline trajectories of both the single gene models. I created a new dataset to predict on; and then I used the rstanarm::linpred() function to predict the values on the link scale for each single-gene model:
rstanarm::posterior_linpred(
object = model, # the rstanarm model object
newdata = newdata_bayes, # the newdata df
re.form = NA # ignore random groups
)Then I decided to plot the trajectories for Sex = F and Seq.Site = YALE and area = DFC:
spline_traj_single_gene
Figure 5: Genewise trajectories predicted in the linear space (Negative Binomial link). These trajectories contain predictions from the former two gene-wise models. Light shaded areas represent 50% credible intervals (25–75% posterior quantiles)
It is evident that the two randomly selected genes show different trajectories, leading to a variability that was not captured by the complete pooling approach and that was diagnosed with a poor fit of the first model. In this plot, however, the variability of the prediction is not evident but would become detectable if I had used the rstanarm::posterior_epred() function.
In summary, the complete pooling strategy oversimplified the data structure overlooking the gene-wise heterogeneity and resulting in poor model fit. On the other hand, the single gene model (no pooling) achieved a better fit, especially with genes with higher scores (see Figure 2). However, assessing the trajectory landscape of schizophrenia using single gene models would be very hard, due to the high number of genes associated with this clinical trait (\(\sim500\)).
4.3 Partial Pooling
To overcome the limitation of both no pooling and complete pooling strategies, I followed the partial pooling path. However, including an interaction term to capture the gene-wise trajectory specificities (\(\sum_{k=1}^{4} B_k*\mathrm{gene}\)) would have been computationally intensive, as it would have resulted in a highly dimensional parameter space (\(\mathrm{d}>2000\)). Thus, I decided to cluster the genes before fitting the bayesian model using a coexpression network analysis (WCGNA). This approach has the advantage of clustering the genes based on their expression similarities. I saw this as an opportunity to create a model that was a trade-off between the complete and no pooling approach.
4.3.1 Step 1: Gene Co-Expression Definition
The package BioNERO (Almeida-Silva and Venancio 2022) was used to perform the network analysis and identify network of CoExpression. These networks will be be included in the hierarchical model, leading to a biologically informed partial pooling approach. I followed the pipeline specified in the Bioconductor vignette which performs a Quality Check (QC) and, subsequently, extracts the modules.
The analysis revealed 9 identifiable modules, with a number of genes spanning from \(\sim30\) to \(>100\)
source("Network_analysis.R")
plot_genes_module + scale_fill_viridis_d(option = "cividis") + theme(plot.title = element_text(hjust = 0.5))
Figure 6: Gene co-expression modules ordered by the number of genes within each module.
and a variable, but acceptable, stability across resampling
module_stability(counts_preprocessed, net, nRuns = 20) + # Assessing modules stability
scale_fill_viridis_d(option = "cividis") +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "none"
)
Figure 7: Module stability across 20 resampling cycles. Within each cycle a new model estimation is computed with a resampling strategy.
For each gene that passed the QC pipeline, I extracted its module membership. This information was used to fit the cluster-wise model:
gene_membership <- net$genes_and_modules %>%
mutate(across(c(Modules, Genes), as.factor))
head(gene_membership, n= 5)## Genes Modules
## 1 ENSG00000002822 6
## 2 ENSG00000003509 7
## 3 ENSG00000004534 1
## 4 ENSG00000005483 5
## 5 ENSG00000006744 54.3.2 Step 2: Partial Pooling Model
The next step aimed to inform the longitudinal model using WCGA modules, seeking a trade-off between the rigidity of complete pooling and the computational burden and interpretational complexity of a gene-wise no pooling approach.
4.3.2.1 Clusterwise Model - Hard Partial Pooling
I firstly decided to extend the complete pooling model including an interaction term between the \(\mathrm{modules}\) and the window term \(\sum_{k=1}^{4} \beta_k B_k(window)\). Also for the interaction terms, the priors specified in rstanarm were weakly informative: \(\mathcal{N}(\mathbf{0}, 2.5)\). These new coefficients were intended to capture trajectory variability between genes, parsimoniously accounted by the WGCA modules. However, this model did not include any gene-specific parameter, resulting in what I refer to as a hard partial pooling approach. The rstanarm code was:
rstanarm::stan_glmer(
expression ~ ns(Window, df=4) * Modules + Sex + Sequencing.Site + area + # fixed effects
(1 | subject), # random effects
data = raw_counts,
family = neg_binomial_2,
offset = offset,
chains = 4,
cores = cores,
adapt_delta = 0.8, # this is the rstan default value
control = list(max_treedepth =10), # this is the rstan default value
algorithm = "sampling"
)The mixing was appropriate, however the PPC analyses revealed a non-sufficient fit:
wrap_plots(list(median_clusterwise_interaction, sd_clusterwise_interaction, sum_0s_clusterwise_interaction)) +
plot_annotation(title = 'PPC in partial pooling model',
theme = theme(plot.title = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0)
))
Figure 8: Fit of the model that includes an interaction term between modules and the trajectory spline component.
This model captured gene variability through the modules; however, it did not allow for variability among genes within the same module. This lack of flexibility in modelling gene-specific variability may explain why the hard partial-pooling approach was not effective in reproducing the original data. In fact, the WGCA analyses are based on the co-expression matrix of the genes across the biological samples and clusters the genes based on their expression correlation. Thus, translating this to the BrainSpan context, genes assigned to the same cluster could exhibit similar developmental expression trajectories while differing in their mean expression levels (see Figure 9 for a conceptual example). For this reason, constraining all the genes of the module to follow the same trajectory, may have hindered the ability of this model to adapt to the data structure.
Figure 9: Conceptual illustration showing that co-expression modules group genes with similar developmental trajectories but different mean expression levels. This within-module heterogeneity may limit the ability of the hard partial-pooling model to reproduce the observed data.
4.3.2.2 Integrating Gene-Wise Variability - Soft Partial Pooling
I named the next model “soft partial pooling” as it was designed to account for the within module heterogeneity (Figure 9). In particular, the model specified a gene specific intercept thereby adding a new level of hierarchy. Each module was characterized by a specific shrinkage distribution governing the intercepts of its constituent genes. Thus, each gene has a specific intercept reagulated by the module-level shrinkage hyperprior: \[y_{ij} \sim \mathrm{NegBin}(\mu_{ij}, \phi) \\ \mu_{ij} = \exp(\eta_{ij}) \\ \eta_{ij} = \alpha + f(B_k)^\top\beta^{(m)} + \gamma^{(m)}_{g} + u_j + log(\mathrm{offset_{ij}})\] Where \(f(B_k)\beta^{(m)}\) is the spline-module interaction \(\sum_{k=1}^{4} B_k^\top \text{window}_{ijk} * \mathrm{Modules}\). The \(\gamma^{(m)}_{g}\) is the gene specific interaction, \(u_j\) is the random intercept for each \(j^{th}\) subject and \(log(\mathrm{offset_{ij}})\) is the normalization factor. \[\alpha \sim \mathcal{N} (0,2) \\ f(B_k)^\top\beta^{(m)} \sim \mathcal{N}(\mu_m, \sigma_m^2) \quad \\ \mu_m \sim \mathcal{N}(0,1)\\ \sigma_m^2 \sim \mathcal{N}^+(0,0.5) \\ \gamma^{(m)}_{g} \sim \mathcal{N}(0,\sigma^2_{g[m]}) \\ \sigma^2_{g[m]} \sim \mathcal{N}^+(0,0.5) \\ u_j \sim \mathcal{N}(0,1) \\ \phi \sim ~ \exp(1)\] The \(\gamma^{(m)}_g\) parameter is a second hierarchy nested within the modules. In particular, the shrinkage of the gene \(g\) in the module \(m\) depends by the shrinkage parameter \(\sigma^2_{g[m]}\). In this setting, modules with higher precision will shrink genes more than modules with lower precision. For simplicity, I did not include the control covariates (\(\mathrm{Sex}\) and \(\mathrm{Seq.Site}\)) that I had included in the former models.
This model specification added a further level of flexibility leading to an increase of fit, showing that the gene-level variability is crucial in determining a reasonable fit of the data.
Figure 10: PPC of the soft partial pooling model with a gene-wise intercept nested into modules. Adding a degree of flexibility increases the fit within a parsimonious setting.
More in depth, the PPC showed a non optimal, but sensibly increased, level of fit, reducing the distance between the \(\mathrm{y_{rep}}\) statistic distribution and the true \(\mathrm{y}\) statistic values. In fact, the difference between the \(y_{rep}\) median and the true median was much smaller than the same difference in the previous partial pooling model (see Figure 8). However, the model overestimated the sd and underestimated the Sum of 0s.
Another useful convergence check is to compare the density of each chain. In particular, if the chains converged on the stationary distribution their density should show substantial overlap. The overlay plot for the parameters of interest showed acceptable overlaps, with lower mixing for the \(\alpha\) parameter:
Figure 11: Density plot for each chain after the warmup.
Although some parameters exhibited ESS values close to the recommended threshold of 400 for four chains, as suggested in the Stan diagnostics guidelines (Stan Development Team 2025). In particular, selecting \(\alpha\) and \(f(B_k)^\top\beta{(m)}\) as parameters of interest, I could extract their ESS values using:
par <- c( # list of variables of interest
"alpha",
as.vector(outer(1:9, 1:4, \(i, j) paste0("beta_module[", i, ",", j, "]")))
)
fit$summary(variables = par) %>%
dplyr::filter(rhat > 1.01 | ess_bulk < 400)## # A tibble: 27 × 4
## variable rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl>
## 1 alpha 1.02 149. 362.
## 2 beta_module[1,2] 1.01 354. 906.
## 3 beta_module[2,2] 1.01 385. 1066.
## 4 beta_module[4,2] 1.01 384. 892.
## 5 beta_module[5,2] 1.01 368. 876.
## 6 beta_module[6,2] 1.01 348. 819.
## 7 beta_module[7,2] 1.01 340. 779.
## 8 beta_module[8,2] 1.01 364. 915.
## 9 beta_module[9,2] 1.01 393. 960.
## 10 beta_module[1,3] 1.02 350. 762.
## # ℹ 17 more rows
Figure 12: Plot of the Neff proportion for each variable. The lower bound for optimal model is 10%, marked as the first vertical dashed line.
The \(\mathrm{Neff}\) ratio plot was in line with the overlay plot (Figure 11), showing suboptimal overlap across chains by the \(\alpha\) parameter, which is reflected in a low effective sample size. In general, the relatively low ESS in the absence of divergences may be due to the difficulties in model identifiability which are particularly evident for \(\alpha\), as a gene specific intercept is included in the model (\(\gamma^{(m)}_{g}\)).
Briefly, the partial pooling strategy consisted of two steps: first, I identified the co-expression clusters; second, these clusters were incorporated into the Bayesian model. This strategy represented the best approach for accounting the complexity of BrainSpan dataset, with the soft-partial pooling approach substantially improving model fit. This finding is consistent with previous literature, suggesting that the gene-wise flexibility is a key point in accomodating the \(\mathrm{NegBin}\) distributional assumption. Moreover, the cluster-wise approach provided gene-specific flexibility while limiting the number of model parameters. This is achieved through the incorporation of prior biological knowledge, in this case represented by the gene’s module membership. Finally, the residual under-fitting of the model may be attributed by the fixed dispersion parameter \(\phi\) of the \(\mathrm{NegBin}\) distribution, which constrained all modules and genes to be equally dispersed.
Since no model showed a robust and reliable solution, I could not move to the second aim of the general objectives and evaluate the impact of the normalization strategies in the transcriptome trajectories definition.
5 Final Remarks and Future Implementations
This work had two objectives: first, to introduce a novel Bayesian framework for modelling transcriptomic trends and, then, evaluate the impact of the normalization strategies on the transcriptomic trajectories (General Objectives).
The modelling process revealed that the ROI structure was highly dependent and mostly accounted by the subject hierarchy, leading to non-convergence of the first model. This is in line with previous spatial transcriptomic works that highlighted the prevalence of temporal differentiation over regional specificities in the cortical development (Kang et al. 2011; Liu and Shimogori 2026). Moreover, the gene-wise flexibility was a key point in improving the model fit and preserving the \(\mathrm{NegBin}\) distributional assumption. Overall, identifying a good strategy to find a trade-off between model complexity and interpretability over the gene-wise specificity represented the main challenge for a useful and interpretable model. The current solution is not robust enough for an actual implementation but may represent a starting point for further investigations. Indeed, a good modelling strategy is an essential requirement for a thorough discussion on the impact of the normalization strategies in the trajectory context.
During this report, I focused on modelling the parameter \(\mu_{ij}\). However, future developments could also model the \(\phi\) dispersion parameter of the \(\mathrm{NegBin}\) distribution. Allowing module-specific dispersion would better account module variability that, in turn, might improve the fit of the soft partial pooling model.
Ultimately, different modelling strategies and biological information could be incorporated as prior information. One possible extension would use the gene hubness to inform the \(\sigma^2_{g[m]}\) dispersion parameter. Conceptually, genes that are more representative of their module (higher hubness) will be expected to have lower trajectory variability around the module-level trajectory. This could be modeled as: \[ \gamma^{(m)}_{g} \sim \mathcal{N}(0,\sigma^2_{g[m]}) \\ \sigma^2_{g[m]} \sim \mathcal{N}^+(0,\tau_{g[m]}^{-1}) \\ \tau_{g[m]} = 2 + h_{g[m]} \] Where \(\gamma^{(m)}_{g}\) represents the deviation of the gene \(g\) from the mean trajectory of module \(m\); \(\tau_{g[m]}\) controls the dispersion of the gene-specific variability \(\sigma^2_{g[m]}\), and \(h_{g[m]}\) represents the hubness of the gene \(g\) within the module \(m\). The \(2\) scale on the \(\tau_{g[m]}\) coefficient is intended to preserve the baseline shrinkage of the prior distribution of \(\sigma^2_{g[m]}\) parameter. Specifically, when \(h_{g[m]} = 0\) the prior reduces to the original \(\mathcal{N}^+(0,0.5)\), whereas increasing \(h_{g[m]}\) increases the precision of the distribution, resulting in stronger shrinkage and a more concentrated prior around zero - the module mean trajectory.
This represents only one of many possible modelling strategies for incorporating prior biological information into the model. I see biostatistics as an ideal playground where different types of data and analytical results could be used in the interplay between prior definition and model outcomes.