Time Course RNA-Seq: Running the Analysis on Your Own Samples

By the end of this you will have a gene-by-timepoint expression matrix from your own blood draws. You will also have a list of genes whose expression changes over time at a controlled false discovery rate, fitted trajectories for each of those genes, and a set of clusters that group genes with similar shapes. You need the following.

• FASTQs from at least six timepoints (eight or more is better)
• a transcriptome index
• ~32 GB RAM and 16 cores for quantification
• R ≥ 4.3 with these packages:
  • tximport
  • DESeq2
  • limma
  • splines
  • DEGreport
• a metadata table for every sample, recording the draw time to the minute, the library prep batch, the RIN, and a same-day complete blood count if you have one

That last item matters more than most of the pipeline.

1. Fix the sampling schedule before you spend money on libraries

Time course design is where most of the statistical power is won or lost, and no amount of modeling recovers a badly spaced series. Three decisions:

Resolution. You cannot resolve a process faster than roughly twice your sampling interval. If you are watching an acute perturbation such as a vaccine, an infection, a diet change, or a training block, the interesting transcriptional response happens in hours to days. Daily or twice-daily sampling for the first week then beats weekly sampling for two months. If you are building a baseline for slow drift, monthly for a year is reasonable. Mixed designs work: dense early, sparse late. The dense-early pattern matches the biology of most perturbations, where an immediate response and a later adapted state are different phenomena and get separated only if you sample both regimes ¹.

Number of timepoints vs. replicates. With a single individual, timepoints are your sample size. Going from 6 to 12 timepoints buys you more than duplicating libraries at 6. Duplicate one or two timepoints anyway, as a technical-variance anchor you can use to sanity-check dispersion estimates.

Clock time. Draw at the same time of day, ideally within a 30-minute window, fasted or consistently non-fasted. Circadian rhythm is one of the largest time-dependent signals in the blood transcriptome and is precisely the kind of periodic structure time series methods are built to detect ². If you sample at 7am on Monday and 4pm on Thursday, you have built a circadian experiment and labeled it something else.

Record everything in a flat TSV now:

sample  fastq_r1            fastq_r2            time_h  draw_clock  batch  rin   neut_pct  lymph_pct
T00     T00_R1.fastq.gz     T00_R2.fastq.gz     0       07:12       A      8.9   58.2      31.1
T06     T06_R1.fastq.gz     T06_R2.fastq.gz     6       13:05       A      8.7   62.0      28.4
...

2. Quantify with salmon, with decoys and bias correction on

Use selective alignment against a decoy-aware index. Build once:

# GENCODE v46, human
grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) \
  | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v46.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz

salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_idx_v46 \
  -k 31 -p 16 --gencode

Then quantify each sample identically. Do not vary flags across timepoints, because a flag change is a batch effect.

salmon quant -i salmon_idx_v46 -l A \
  -1 T00_R1.fastq.gz -2 T00_R2.fastq.gz \
  --gcBias --seqBias --posBias \
  --numGibbsSamples 20 \
  --validateMappings -p 16 \
  -o quants/T00

--gcBias matters when libraries were prepped on different days, which in a longitudinal series they always were. --numGibbsSamples 20 gives you per-transcript inferential replicates, which you need if you ever want transcript-level or differential-transcript-usage analysis (fishpond/swish). Check quants/*/logs/salmon_quant.log for mapping rate. Whole blood with globin present will run 60–75%. PBMC or globin-depleted whole blood should be above 80%. A single sample 15 points below the rest is a red flag and not noise.

3. Import to gene level and filter

library(tximport); library(DESeq2); library(readr)
meta <- read_tsv("metadata.tsv")
files <- file.path("quants", meta$sample, "quant.sf"); names(files) <- meta$sample
tx2gene <- read_tsv("tx2gene_v46.tsv")  # TXNAME, GENEID

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

dds <- DESeqDataSetFromTximport(txi, colData = as.data.frame(meta),
                                design = ~ 1)
keep <- rowSums(counts(dds) >= 10) >= ceiling(0.5 * ncol(dds))
dds <- dds[keep, ]

Filter with the time structure in mind. A gene expressed only at three consecutive timepoints and zero elsewhere is a real transient, not noise, and a blanket “expressed in ≥50% of samples” filter deletes exactly the genes you sampled densely to find. We use two filters and take the union: ≥10 counts in ≥50% of samples, or ≥10 counts in ≥3 consecutive timepoints.

Before modeling, look at the data. vst(dds, blind=TRUE) then PCA, colored by time and shaped by batch. In a clean series, PC1 or PC2 traces a visible arc through time. If PC1 separates batches and time is buried at PC4, fix that before continuing (see step 5).

4. Choose a model that matches the shape you expect

There are four families in common use, and the choice matters more than the software ³.

Categorical time (time as a factor). Fit ~ batch + factor(time_h) and test with a likelihood ratio test against ~ batch. This makes no assumption about trajectory shape and is the right default when you have few timepoints and several replicates per timepoint. With one replicate per timepoint it is useless: the model saturates and there are zero residual degrees of freedom.

Natural splines. This is our default for a single-subject series with one library per timepoint. You spend df parameters on time instead of n_timepoints - 1, which leaves residual degrees of freedom for dispersion estimation.

library(splines)
tt <- meta$time_h
design(dds) <- ~ batch + ns(tt, df = 3)
dds <- DESeq(dds, test = "LRT", reduced = ~ batch)
res <- results(dds, alpha = 0.05)
sum(res$padj < 0.05, na.rm = TRUE)

Pick df from the sampling density, not from the result: df = 3 for 8–12 timepoints, df = 4–5 for 15–25. If time is unevenly spaced, place knots at quantiles of the observed times and not at even intervals: ns(tt, knots = quantile(tt, c(.33, .66))). Over-flexible splines fit RNA extraction noise and produce trajectories that wiggle beautifully and mean nothing.

Parametric impulse models. If the perturbation is a single stimulus, the expected shape is a sigmoid rise and a sigmoid decay back toward a new plateau. ImpulseDE2 fits exactly that six-parameter form and gains power over splines when the shape assumption holds, and it will also classify genes as transient vs. monotonic (a genuinely useful annotation for a stimulus experiment).

Nonparametric Gaussian processes. A GP over time with a squared-exponential kernel handles irregular spacing, missing timepoints, and noisy observations without committing to a shape, and it gives a posterior over the trajectory rather than a point estimate. GP-based tests were developed and validated on exactly this problem, including human Th17 differentiation time series where they identified dynamically regulated genes that per-timepoint pairwise testing missed ⁴. The cost is compute and the need to think about kernel hyperparameters. We reach for GPs when the sampling is irregular enough that splines are awkward.

A broader survey of these families, including hidden Markov and state-space approaches for longer series, is worth reading before you commit ²⁵.

5. Handle the confounders that will otherwise dominate

Cell composition. In whole blood or PBMC, most of what looks like a transcriptional response is a shift in cell proportions. Neutrophil fraction alone swings the apparent expression of thousands of genes. If you have a same-day CBC, put the proportions in the model:

design(dds) <- ~ batch + neut_pct + lymph_pct + ns(tt, df = 3)

If you do not, deconvolve: CIBERSORTx with LM22, or granulator, on the VST matrix, then include the top one or two estimated fractions as covariates. Adjusting removes the composition signal from your gene-level result, so keep the fractions as their own time series and analyze them separately. Longitudinal immune profiling platforms treat composition as a first-class readout rather than a nuisance, and that is the right instinct ⁶.

Batch and prep date. If you sequenced three timepoints in January and five in June, batch is partly confounded with time and no covariate fully separates them. Either randomize prep so each batch spans the time range, or freeze RNA and prep all libraries together. The latter is what we do for personal series: collect into PAXgene or Tempus tubes, freeze at −80°C, prep everything in one batch at the end. RNA is stable in those tubes for years at −80.

Globin. Whole-blood RNA is dominated by HBB, HBA1, HBA2. Without depletion, a large share of reads are globin and effective library depth collapses. Use a globin depletion kit at prep, or accept the loss and sequence deeper. If globin fraction varies across timepoints, it drives normalization artifacts, so compute per-sample globin percentage and check whether it correlates with time.

Degradation. RIN below 7 changes 3′ coverage bias and shows up as a spurious expression-vs-time trend if RIN drifted over the collection period. Regress on RIN if it varies, and record it always.

6. Cluster trajectories, do not eyeball them

Once you have a few thousand significant genes, the useful output is a small number of shapes. Standardize each gene’s fitted trajectory (z-score across time) and cluster.

DEGreport::degPatterns is the fastest reasonable path:

library(DEGreport)
vsd <- vst(dds, blind = FALSE)
sig <- rownames(res)[which(res$padj < 0.05)]
clusters <- degPatterns(assay(vsd)[sig, ], metadata = as.data.frame(meta),
                        time = "time_h", col = NULL, minc = 15)

For anything larger, TCseq (fuzzy c-means on standardized trajectories) is better, because soft assignment tells you which genes sit between two shapes instead of forcing them. The maSigPro package fits polynomial regression per gene and then groups, and remains a solid choice for multi-condition designs.

Matrix factorization is the alternative worth knowing. Instead of assigning each gene to one cluster, decompose the time-by-gene matrix into patterns that genes load onto continuously. CoGAPS applied to time course data separated an immediate therapeutic response from a later acquired-resistance program in the same series, which is a distinction hard clustering blurs ¹. If your experiment has a plausible early phase and late phase, run both and compare.

Then annotate: clusterProfiler::enrichGO or fgsea on each cluster’s gene set against the background of expressed genes. Do not use the whole genome as the background. Ranked GSEA on the spline LRT statistic is generally more informative than over-representation on a thresholded list.

7. Connect the transcript series to your other time series

The reason to run RNA-seq longitudinally rather than once is that trajectories can be aligned with other trajectories. If you have proteomics on the same draws, correlate each gene’s spline fit against its protein’s fit with a lag, since protein follows transcript by hours to days and the lag itself is informative. If you have CGM, summarize each day into features (mean glucose, MAGE, time in a band) and correlate those daily series against transcript clusters. If you have serial CBCs and metabolic panels, they anchor the composition adjustment from step 5.

Methods for this are still moving, and most integration papers use one of three approaches. They concatenate features and apply a joint factorization, model each modality separately and correlate the latent trajectories, or fit an explicit dynamic model with cross-modality edges ⁵⁷. For a first pass on one person’s data, correlate latent trajectories. It is interpretable and it fails loudly. Multi-omic time series in other systems have been analyzed this way with the different molecular layers resolved onto a common time scale before integration ⁸.

None of this output is a clinical result. Differential expression in your own blood over a few timepoints is a hypothesis about your biology. Anything that looks like it might matter for your health belongs in front of a physician with the raw values. A self-directed decision is not the place for it.

Common problems

Every gene is significant. Usually batch-time confounding, or a spline with too many degrees of freedom and too few residual degrees of freedom. Check ns() df against ncol(dds). If df + n_covariates + 1 is close to your sample count, dispersion estimates collapse and p-values are meaningless. Also re-plot the PCA colored by neutrophil fraction.

Nothing is significant. Check the mapping rates and the library sizes first. Then check whether you sampled during the window when something changed. A weekly series through an event that resolves in 48 hours will show a flat line.

Trajectories that wiggle at high frequency. Over-flexible model, or one outlier timepoint pulling a knot. Refit with df = 3, and look at Cook’s distance: assays(dds)[["cooks"]].

Irregular or missing timepoints. Splines and GPs handle this natively, and factor-time models do not. Do not interpolate counts to fill gaps. Put the real times in ns(tt, ...) and let the model do the work.

Uncertain about ordering across a long series. If samples span months and the biological process is not clock-synchronized, consider pseudotime ordering as a check on the nominal ordering before you trust the calendar axis ³.

Low-count genes dominating the significant list. Add lfcShrink(dds, coef = ..., type = "apeglm") for effect sizes, and require a minimum fold change across the fitted trajectory (say, 1.5× between the trajectory minimum and maximum) on top of the FDR threshold.

Footnotes

1. Genevieve Stein-O’Brien, Luciane T. Kagohara, Sijia Li, et al. Integrated time course omics analysis distinguishes immediate therapeutic response from acquired resistance. Genome Medicine, 2018. https://doi.org/10.1186/s13073-018-0545-2 ↩ ↩²

2. Vera-Khlara S. Oh, Robert W. Li. Temporal Dynamic Methods for Bulk RNA-Seq Time Series Data. Genes, 2021. https://doi.org/10.3390/genes12030352 ↩ ↩²

3. Daniel Spies, Constance Ciaudo. Dynamics in Transcriptomics: Advancements in RNA-seq Time Course and Downstream Analysis. Computational and Structural Biotechnology Journal, 2015. https://doi.org/10.1016/j.csbj.2015.08.004 ↩ ↩²

4. Tarmo Äijö, Vincent Butty, Zhi Chen, et al. Methods for time series analysis of RNA-seq data with application to human Th17 cell differentiation. Bioinformatics, 2014. https://doi.org/10.1093/bioinformatics/btu274 ↩

5. Yulan Liang, Arpad Kelemen. Dynamic modeling and network approaches for omics time course data: overview of computational approaches and applications. Briefings in Bioinformatics, 2017. https://doi.org/10.1093/bib/bbx036 ↩ ↩²

6. Suhas V. Vasaikar, Adam K. Savage, Qiuyu Gong, et al. A comprehensive platform for analyzing longitudinal multi-omics data. Nature Communications, 2023. https://doi.org/10.1038/s41467-023-37432-w ↩

7. Ali Reza Taheriyoun, Allen Ross, Abolfazl Safikhani, et al. Longitudinal omics data analysis: approaches and applications. Computational and Structural Biotechnology Journal, 2026. https://doi.org/10.1016/j.csbj.2026.01.001 ↩

8. Joan Márquez‐Molins, Pascual Villalba‐Bermell, Julia Corell‐Sierra, et al. Integrative time‐scale and multi‐omics analysis of host responses to viroid infection. Plant, Cell & Environment, 2023. https://doi.org/10.1111/pce.14647 ↩