Differentiation timing score

Differentiation timing score

PCA-based scoring of differentiation timing from bulk RNA-seq time-course data

Author: Zoheb Khan

Created: 2025-08-29

Modified: 2026-10-12

Compiled: 2026-10-12

Environment: R 4.6.1; Bioconductor 3.23

Dependencies: AnnotationDbi 1.74.0; bookdown 0.47; circlize 0.4.18; clusterProfiler 4.20.0; ComplexHeatmap 2.28.0; DESeq2 1.52.0; edgeR 4.10.1; ggplot2 4.0.3; ggrepel 0.9.8; knitr 1.51; org.Hs.eg.db 3.23.1; patchwork 1.3.2; ragg 1.5.2; scales 1.4.0; systemfonts 1.3.2; viridis 0.6.5

Source: github.com/ZohebKhan1/pca-maturation-scoring

1 Introduction

Differentiation timing is estimated by projecting bulk RNA-seq samples onto a reference trajectory in principal component analysis (PCA) space. PCA is fitted to the mean expression of temporal genes at each reference timepoint, and the resulting centroids are connected in time order. Each sample’s position on the trajectory gives a predicted reference time and a differentiation score from 0 to 1. The projection distance measures fit to the reference.

The analysis uses a cardiomyocyte differentiation time course from GSE122380, with leave-one-cell-line-out validation to assess prediction across cell lines.

2 Function inputs

get_temporal_genes() requires integer-like raw counts and variance-stabilized (VST) expression with matching gene and sample identifiers. score_differentiation_timing() uses normalized expression and the selected temporal genes. In both matrices, rows are genes and columns are samples. Row and column names must be unique, with one metadata row per sample.

The default metadata columns are sample_id and day_numeric; other names can be specified through sample_id_col and time_col. The reference must span at least two distinct, finite numeric timepoints. Non-reference samples can have missing timepoints.

The scorer uses all supplied samples as the reference unless reference_samples specifies a subset. For comparisons across experimental conditions, gene selection and trajectory fitting are restricted to the reference cohort.

Gene selection requires DESeq2 and edgeR:

install.packages('BiocManager')
BiocManager::install(c('DESeq2', 'edgeR'))

3 Reference dataset

GSE122380 contains 192 bulk RNA-seq samples from 13 human induced pluripotent stem cell lines collected over days 1 through 15 of cardiomyocyte differentiation. The metadata, raw counts, and VST expression are included in the repository.

The code runs from the repository root. load_GSE122380_data() aligns the count and VST matrices to the metadata sample order:

source('data/load_GSE122380_data.R')
source('functions/get_temporal_genes.R')
source('functions/score_differentiation_timing.R')

reference_data <- load_GSE122380_data()
metadata <- reference_data$metadata
counts <- reference_data$counts
vst <- reference_data$vst

Table 3.1: Samples and represented cell lines by differentiation day.

Metric	D1	D2	D3	D4	D5	D6	D7	D8	D9	D10	D11	D12	D13	D14	D15
Samples	13	12	12	12	13	13	13	13	13	13	13	13	13	13	13
Cell lines	13	12	12	12	13	13	13	13	13	13	13	13	13	13	13

PCA of the most variable genes and correlations between day-mean expression profiles show the temporal structure of the reference (Figure 1).

Figure 1. GSE122380 reference time course and day-level expression correlation. a. PCA using the top 10% most variable VST genes, colored by differentiation day. b. Pearson correlation between day-mean profiles from the same variable-gene set.

The scoring model uses temporal genes and their mean expression at each timepoint.

4 Scoring workflow

4.1 Temporal-gene selection

get_temporal_genes() applies three sequential filters within the reference cohort. Raw counts are used for the expression filter and DESeq2 likelihood-ratio test (LRT), and VST expression is used for the range filter. The trajectory model includes genes that pass all three:

1. Expression filter. Counts are normalized by the trimmed mean of M-values (TMM) method and expressed as counts per million (CPM). Genes are retained if their highest mean CPM across reference timepoints reaches the expression threshold.
2. DESeq2 likelihood-ratio test. The full model includes categorical time and any requested adjustment covariates. The reduced model contains the adjustment covariates without time. Genes with an adjusted LRT p-value below the selected threshold continue to the final filter.
3. VST-range filter. Genes are retained if the range of their mean VST expression across reference timepoints reaches the selected threshold.

The LRT tests whether adding categorical time improves the count model. If the full and reduced models have maximized log likelihoods ell_full and ell_reduced, the likelihood-ratio statistic is

\[ \mathrm{LR} = 2(\ell_\mathrm{full} - \ell_\mathrm{reduced}). \]

For GSE122380, the full model is ~ cell_line + day_numeric and the reduced model is ~ cell_line. Time is converted to a categorical factor, allowing nonlinear expression changes across days. Including cell_line in both models accounts for differences among cell lines. With no adjustment covariates, the reduced model is ~ 1.

source('functions/get_temporal_genes.R')

temporal_selection <- get_temporal_genes(
  raw_counts = counts,
  vst_expression = vst,
  metadata = metadata,
  sample_id_col = 'sample_id',
  time_col = 'day_numeric',
  adjustment_covariates = 'cell_line'
)

temporal_genes <- temporal_selection$temporal_genes

reference_group_col and reference_group_value restrict all three filters to a specified group of reference samples.

The default thresholds are a maximum day-mean TMM CPM of at least 10, an adjusted LRT p-value below 1e-7, and a day-mean VST range of at least 0.6. The number of genes retained at each step is shown in Table 4.1. Thresholds for another study depend on its library sizes, sample size, timepoint coverage, and experimental design.

Table 4.1: Sequential temporal-gene filtering.

Metric	Value
Genes in the aligned count/VST input	13,615
Genes with maximum day-mean TMM CPM >= 10	11,411
Genes passing LRT adjusted p-value < 1e-7	10,863
Genes passing all filters and used for scoring	8,215

Figure 2 shows the first 1,500 temporal genes in LRT result order. PCA, scoring, and validation use the complete temporal-gene set.

Figure 2. Heatmap of the 1,500 displayed temporal genes across the reference samples, ordered by differentiation day.

The heatmap genes form four temporal expression clusters (Figure 3). Gene Ontology (GO) Biological Process enrichment is calculated for each cluster using the 1,500 displayed genes as the enrichment universe. These analyses describe the temporal expression patterns; scoring uses all selected temporal genes together.

Figure 3. Temporal clusters and GO enrichment for the 1,500 heatmap genes. Each row pairs a cluster trajectory with its top enriched GO Biological Process terms.

4.2 Reference PCA

The reference PCA is fitted to mean temporal-gene expression at each timepoint. Each timepoint contributes one expression profile to the fit, regardless of its sample count.

Let \(\mathbf{x}_i\) be the normalized expression vector for sample \(i\), restricted to the selected temporal genes. For a reference timepoint \(t\), let \(\bar{\mathbf{x}}_t\) be the mean vector across its reference samples:

\[ \bar{\mathbf{x}}_t = \frac{1}{n_t}\sum_{i:\,t_i=t} \mathbf{x}_i. \tag{1} \]

PCA is centered without variance scaling. The smallest number of principal components (PCs) explaining at least 99% of the variance between timepoint means is retained. With \(\boldsymbol{\mu}\) denoting the mean of those profiles and \(W\) the retained loading vectors, the PCA coordinates for sample \(i\) are

\[ \mathbf{z}_i = (\mathbf{x}_i - \boldsymbol{\mu})W. \tag{2} \]

Reference and experimental samples are projected using the same centering and loadings. Seven PCs meet the 99% variance target in GSE122380; the retained number depends on the reference cohort and temporal-gene set.

Figure 4. Individual samples projected into centered, unscaled PCA spaces fitted to reference timepoint means. The full temporal-gene set defines the scoring model. The C1+C2 and C3+C4 subsets show how early and late expression patterns contribute to the time-course structure.

Figure 5 shows genes with the largest positive and negative PC1 and PC2 loadings, with GO enrichment for the strongest PC1 loadings.

Figure 5. PCA loading genes and GO enrichment. a-b. The ten largest positive and negative PC1 and PC2 loading genes. Ambiguous Ensembl-to-symbol mappings retain the Ensembl ID. c-d. GO enrichment among the 500 strongest negative and positive PC1 loadings.

PC1 is oriented toward increasing differentiation day in the figures. This choice of axis sign preserves distances and timing scores.

4.3 Trajectory fitting and sample projection

The reference timepoint means define centroids in the retained PCA space. Connecting adjacent centroids in time order produces a polyline trajectory that can bend across multiple PCs.

For one segment, let \(\mathbf{a}\) and \(\mathbf{b}\) be its start and end centroids, and let \(\mathbf{z}\) be a sample’s PCA coordinate. The segment direction is \(\mathbf{v}=\mathbf{b}-\mathbf{a}\). The closest point on the finite segment is

\[ \alpha = \operatorname{clip}_{[0,1]} \left( \frac{\langle \mathbf{z}-\mathbf{a},\mathbf{v}\rangle} {\langle \mathbf{v},\mathbf{v}\rangle} \right), \qquad \mathbf{q}=\mathbf{a}+\alpha\mathbf{v}. \tag{3} \]

The projection is clipped to each finite segment by constraining \(\alpha\) between 0 and 1. The segment minimizing \(\lVert\mathbf{z}-\mathbf{q}\rVert^2\) determines the sample’s assignment, and its endpoint times are used to interpolate predicted time:

\[ \widehat{t}=t_a+\alpha(t_b-t_a), \qquad s=\frac{\widehat{t}-t_{\min}}{t_{\max}-t_{\min}}. \tag{4} \]

Here, \(t_a\) and \(t_b\) are the segment endpoint times, \(\widehat{t}\) is predicted time, and \(s\) is the differentiation score. Scores are bounded between 0 and 1, corresponding to the earliest and latest reference times.

Zero-length segments are omitted. If the trajectory folds back near itself, nearby segments can yield ambiguous time assignments even when projection distances are small.

source('functions/score_differentiation_timing.R')

timing_fit <- score_differentiation_timing(
  expression_matrix = vst,
  metadata = metadata,
  temporal_genes = temporal_genes,
  sample_id_col = 'sample_id',
  time_col = 'day_numeric'
)

timing_fit$scores[, c(
  'sample_id',
  'observed_time',
  'predicted_time',
  'differentiation_score',
  'squared_distance'
)]

Figure 6 displays the first two PCs of the fitted trajectory. Projection distances are calculated across all retained PCs.

Figure 6. Differentiation trajectory in reference PCA space. Black points are day centroids; arrows connect them in day order. The day 1 and day 15 centroids anchor scores 0 and 1.

Time is interpolated separately within each segment, allowing distances per unit time to vary across intervals. Figure 7 shows the reference samples projected back onto the fitted trajectory.

Figure 7. Reference-sample differentiation timing scores by observed day. Dashed lines mark the bounded reference endpoints at 0 and 1.

These scores describe fit within the reference cohort. Prediction in held-out samples is assessed by leave-one-cell-line-out validation.

5 Leave-one-cell-line-out validation

Each cell line is held out in turn. Temporal-gene selection, including the DESeq2 LRT, and PCA trajectory fitting are repeated using the remaining cell lines. The held-out samples are then projected onto that fold’s reference trajectory.

All folds share the supplied VST matrix because the upstream transformation procedure is unavailable in the repository. Validation therefore assesses transfer across cell lines within this processed cohort, conditional on the shared transformation. Transfer to an independently processed study requires separate validation.

Within each fold, genes with positive and negative Spearman correlations with time are labeled maturation and progenitor genes, respectively. Reported accuracy uses the complete temporal-gene set selected in that fold.

source('functions/get_temporal_genes.R')
source('functions/score_differentiation_timing.R')
source('functions/run_loo_maturation.R')

loo_validation <- run_loo_maturation(
  counts = counts,
  vst = vst,
  metadata = metadata
)

loo_validation$summary

The summary includes held-out predictions for all 13 cell lines. Figure 8 displays 12 cell lines in a 3 by 4 layout; predictions for cell line 19190 are included in the summary and the remaining validation figures.

Table 5.1: Leave-one-cell-line-out validation summary.

Metric	Value
Held-out cell lines	13
Held-out samples	192
Correlation between actual and predicted day	0.979
Squared Pearson correlation	0.959
Mean absolute error	0.55 days
Median absolute error	0.30 days
Predictions within 1 day	82.3%
Predictions within 2 days	97.9%

Held-out predictions have a Pearson correlation of 0.979 with observed days, a median absolute error of 0.30 days, and 82.3% of predictions within one day. These estimates apply to GSE122380 under the processing and validation scheme described above.

Figure 8. Held-out predictions for 12 cell lines. Each panel reports Pearson \(r\) and squared Pearson correlation, \(r^2\), to summarize linear association between observed and predicted days.

Pearson \(r\) and \(r^2\) summarize association, while absolute prediction error measures accuracy on the day scale. Figures 9 and 10 show the pooled predictions and errors by day.

Figure 9. Predicted versus observed differentiation day for all held-out samples. The dashed line is identity; the solid line is an ordinary least-squares summary.

For held-out sample \(i\), absolute prediction error is

\[ e_i = |\widehat{t}_i - t_i|. \tag{5} \]

Figure 10 shows mean absolute error by differentiation day, with error bars representing one standard deviation across held-out samples.

Figure 10. Absolute held-out prediction error by day.

6 Interpreting results

6.1 Scores and projection distances

timing_fit$scores contains the sample-level estimates:

• predicted_time gives the sample’s position on the original reference-time scale.
• differentiation_score rescales that position so the earliest and latest reference times are 0 and 1.
• squared_distance measures how close the sample is to its nearest trajectory segment in the retained PCA space.

Scores are specific to the reference lineage, differentiation protocol, and sampled time range. Comparisons across conditions require a shared reference model.

A large projection distance relative to reference samples indicates poor fit to the trajectory and limits the interpretation of predicted time. Differences in cell composition, treatment response, or batch can produce these deviations.

6.2 Reference selection for another study

For comparisons across conditions, temporal-gene selection and trajectory fitting use an independent reference group. The following example identifies that group through a metadata column named group, with reference samples labeled reference. Here, counts, vst, and metadata contain aligned reference and experimental samples from the new study, with the same sample_id and day_numeric columns used above:

temporal_selection <- get_temporal_genes(
  raw_counts = counts,
  vst_expression = vst,
  metadata = metadata,
  reference_group_col = 'group',
  reference_group_value = 'reference'
)

timing_fit <- score_differentiation_timing(
  expression_matrix = vst,
  metadata = metadata,
  temporal_genes = temporal_selection$temporal_genes,
  reference_samples = temporal_selection$reference_samples
)

Experimental samples can have missing timepoints. Adjustment covariates can be included in gene selection when supported by the reference design.

The reference should span the relevant differentiation stages and use gene identifiers and normalization compatible with the experimental samples. Differences in protocol, cell composition, or batch can alter scores and projection distances. Application to a substantially different dataset requires validation in that setting.

7 References

1. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15:550. doi:10.1186/s13059-014-0550-8
2. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–140. doi:10.1093/bioinformatics/btp616
3. Strober BJ, Elorbany R, Rhodes K, et al. Dynamic genetic regulation of gene expression during cellular differentiation. Science. 2019;364(6447):1287–1290. doi:10.1126/science.aaw0040
4. NCBI GEO accession GSE122380

8 Session information

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## time zone: America/New_York
## tzcode source: internal
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets
## [8] methods   base
## other attached packages:
##  [1] viridis_0.6.5               viridisLite_0.4.3
##  [3] systemfonts_1.3.2           scales_1.4.0
##  [5] ragg_1.5.2                  patchwork_1.3.2
##  [7] org.Hs.eg.db_3.23.1         ggrepel_0.9.8
##  [9] ggplot2_4.0.3               edgeR_4.10.1
## [11] limma_3.68.4                DESeq2_1.52.0
## [13] SummarizedExperiment_1.42.0 MatrixGenerics_1.24.0
## [15] matrixStats_1.5.0           GenomicRanges_1.64.0
## [17] Seqinfo_1.2.0               ComplexHeatmap_2.28.0
## [19] clusterProfiler_4.20.0      circlize_0.4.18
## [21] AnnotationDbi_1.74.0        IRanges_2.46.0
## [23] S4Vectors_0.50.1            Biobase_2.72.0
## [25] BiocGenerics_0.58.1         generics_0.1.4
## [27] bookdown_0.47
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3      rstudioapi_0.19.0       jsonlite_2.0.0
##   [4] shape_1.4.6.1           tidydr_0.0.6            magrittr_2.0.5
##   [7] magick_2.9.1            ggtangle_0.1.2          rmarkdown_2.31
##  [10] farver_2.1.2            GlobalOptions_0.1.4     fs_2.1.0
##  [13] vctrs_0.7.3             Cairo_1.7-0             memoise_2.0.1
##  [16] ggtree_4.2.0            htmltools_0.5.9         S4Arrays_1.12.0
##  [19] SparseArray_1.12.2      gridGraphics_0.5-1      sass_0.4.10
##  [22] bslib_0.11.0            htmlwidgets_1.6.4       plyr_1.8.9
##  [25] httr2_1.3.0             cachem_1.1.0            igraph_2.3.3
##  [28] lifecycle_1.0.5         iterators_1.0.14        pkgconfig_2.0.3
##  [31] Matrix_1.7-5            R6_2.6.1                fastmap_1.2.0
##  [34] gson_0.2.0              clue_0.3-68             digest_0.6.39
##  [37] aplot_0.3.1             enrichplot_1.32.0       colorspace_2.1-3
##  [40] ggnewscale_0.5.2        aisdk_1.4.12            textshaping_1.0.5
##  [43] RSQLite_3.53.3          labeling_0.4.3          mgcv_1.9-4
##  [46] httr_1.4.8              polyclip_1.10-7         abind_1.4-8
##  [49] compiler_4.6.1          bit64_4.8.2             fontquiver_0.2.1
##  [52] withr_3.0.3             doParallel_1.0.17       S7_0.2.2
##  [55] BiocParallel_1.46.0     DBI_1.3.0               ggforce_0.5.0
##  [58] MASS_7.3-65             rappdirs_0.3.4          DelayedArray_0.38.2
##  [61] rjson_0.2.23            tools_4.6.1             otel_0.2.0
##  [64] ape_5.8-1               scatterpie_0.2.6        glue_1.8.1
##  [67] callr_3.8.0             nlme_3.1-169            GOSemSim_2.38.3
##  [70] cluster_2.1.8.2         reshape2_1.4.5          gtable_0.3.6
##  [73] tidyr_1.3.2             XVector_0.52.0          foreach_1.5.2
##  [76] pillar_1.11.1           stringr_1.6.0           yulab.utils_0.2.4
##  [79] splines_4.6.1           dplyr_1.2.1             tweenr_2.0.3
##  [82] treeio_1.36.1           lattice_0.22-9          bit_4.6.0
##  [85] tidyselect_1.2.1        locfit_1.5-9.12         fontLiberation_0.1.0
##  [88] GO.db_3.23.1            Biostrings_2.80.1       knitr_1.51
##  [91] gridExtra_2.3.1         fontBitstreamVera_0.1.1 xfun_0.60
##  [94] statmod_1.5.2           stringi_1.8.7           lazyeval_0.2.3
##  [97] ggfun_0.2.1             yaml_2.3.12             evaluate_1.0.5
## [100] codetools_0.2-20        gdtools_0.5.1           tibble_3.3.1
## [103] qvalue_2.44.0           BiocManager_1.30.27     ggplotify_0.1.3
## [106] cli_3.6.6               jquerylib_0.1.4         processx_3.9.0
## [109] dichromat_2.0-0.1       Rcpp_1.1.2              png_0.1-9
## [112] parallel_4.6.1          blob_1.3.0              DOSE_4.6.0
## [115] tidytree_0.4.8          ggiraph_0.9.6           enrichit_0.2.0
## [118] purrr_1.2.2             crayon_1.5.3            GetoptLong_1.1.1
## [121] rlang_1.3.0             KEGGREST_1.52.2