Skip to contents

Get started with rsemflow

Source: vignettes/rsemflow.Rmd
rsemflow.Rmd

rsemflow takes the *.genes.results files written by RSEM and a sample metadata table, and produces the standard downstream results of a bulk RNA-seq experiment as tab-separated tables:

  • expression summaries per sample and per gene,
  • normalized and variance-stabilized expression,
  • PCA and sample correlation,
  • DESeq2 differential expression for a design you specify,
  • preranked GSEA and GSVA/ssGSEA pathway scores,
  • gene annotation.

It does not draw figures. Every table can be opened in R, Python, Excel, or Prism.

What you need 

project/
├── data/
│   ├── S01.genes.results
│   ├── S02.genes.results
│   └── ...
└── metadata.tsv

metadata.tsv is an ordinary table with a sample_id column that matches the file names (S01.genes.results → S01), plus any columns that describe the samples:

sample_id   condition   batch
S01 Control A
S02 Control B
S03 Drug    A
S04 Drug    B

Three commands

Once installed (see vignette("installation", package = "rsemflow")), a first analysis is three commands:

rsemflow data import --i-rsem-dir data/ --m-metadata metadata.tsv --o-study study/
rsemflow expression pca --i-study study/ --o-pca results/pca/
rsemflow differential deseq2 \
  --i-study study/ \
  --p-formula '~ batch + condition' \
  --p-term condition \
  --p-reference condition::Control \
  --o-differential results/differential/

data import checks the input and stores it as a study directory. Every later command reads that study, so the RSEM files are parsed and validated once.

The same analysis in R

Each command is an ordinary R function. Here they run on the synthetic example shipped with the package (6 samples, 200 genes):

library(rsemflow)

ex <- system.file("extdata", "example", package = "rsemflow")
study <- read_rsem_study(
  rsem_dir = file.path(ex, "data"),
  metadata = file.path(ex, "metadata.tsv")
)
study
#> <rsemflow_study>
#>   samples:  6
#>   genes:    200
#>   metadata: condition, batch, age
#>   created:  2026-10-06 00:37:20 UTC (rsemflow 0.1.1)
pca <- compute_pca(study)
#> Warning: replacing previous import 'S4Arrays::makeNindexFromArrayViewport' by
#> 'DelayedArray::makeNindexFromArrayViewport' when loading 'SummarizedExperiment'
pca$variance[1:3, ]
#>   component eigenvalue variance_fraction cumulative_variance
#> 1       PC1  199.94297        0.64846122           0.6484612
#> 2       PC2   50.38480        0.16340953           0.8118708
#> 3       PC3   28.31976        0.09184752           0.9037183
de <- differential_deseq2(
  study,
  formula = ~ batch + condition,
  term = "condition",
  reference = "condition::Control"
)
names(de$results)
#> [1] "condition_Drug_vs_Control"
head(de$results$condition_Drug_vs_Control[, c("gene_id", "baseMean", "log2FoldChange", "padj")])
#>             gene_id baseMean log2FoldChange         padj
#> 1 ENSG00000000001.2 169.1780       2.397434 0.0007999406
#> 2 ENSG00000000002.3 208.8794       2.232664 0.0038441366
#> 3 ENSG00000000003.1 248.5807       2.127022 0.0126338835
#> 4 ENSG00000000004.2 288.2820       2.053341 0.0522060612
#> 5 ENSG00000000005.3 327.9833       1.998960 0.1207709981
#> 6 ENSG00000000006.1 367.6847       1.957152 0.1696442033

The example data are synthetic, so the numbers carry no biological meaning.

How rsemflow treats your design

rsemflow never chooses a model for you. The formula (~ batch + condition), the term to report (condition), and the reference level (Control) are always explicit. In exchange, it stops instead of returning results that cannot be trusted: when sample IDs do not match the metadata, when gene tables differ between files, when a model variable has missing values, or when the design is confounded (for example, every Control sample in batch A and every Drug sample in batch B).

Next steps