This tutorial runs the whole command-line workflow on the synthetic example that ships with rsemflow: six samples (three Control, three Drug, in two batches) and 200 genes. Every output shown below was produced by the commands on this page when the documentation was built.
The data are synthetic; the numbers illustrate the file formats, not biology.
0. The input
All outputs are plain TSV files. To look at them, this page uses a small helper that reads a table and shows its first rows:
show_table <- function(path, n = 6, cols = NULL) {
x <- utils::read.delim(path, check.names = FALSE)
if (!is.null(cols)) x <- x[, cols, drop = FALSE]
utils::head(x, n)
}.
├── data/
│ ├── C1.genes.results C2.genes.results C3.genes.results
│ └── D1.genes.results D2.genes.results D3.genes.results
└── metadata.tsv
show_table("metadata.tsv")
#> sample_id condition batch age
#> 1 C1 Control A 42
#> 2 C2 Control B 51
#> 3 C3 Control A 47
#> 4 D1 Drug B 49
#> 5 D2 Drug A 55
#> 6 D3 Drug B 53Each RSEM file is the standard gene-level output:
show_table("data/C1.genes.results", n = 3)
#> gene_id transcript_id(s) length effective_length expected_count
#> 1 ENSG00000000001.2 ENST00000000001.1 837 657 37
#> 2 ENSG00000000002.3 ENST00000000002.1 874 694 54
#> 3 ENSG00000000003.1 ENST00000000003.1 911 731 71
#> TPM FPKM
#> 1 3222.915 2578.332
#> 2 4504.587 3603.669
#> 3 5682.149 4545.7191. Import
#> Imported 6 samples and 200 genes into study/Import checks that every RSEM file has the required columns, that
file names and sample_id match one-to-one, and that all
files share the same gene table. It then reads the counts with tximport
and writes a study directory:
list.files("study")
#> [1] "effective_length.tsv" "expected_count.tsv" "metadata.tsv"
#> [4] "samples.tsv" "study-info.tsv" "study.rds"
#> [7] "tpm.tsv"
show_table("study/study-info.tsv", n = 10)
#> key value
#> 1 rsemflow_version 0.1.1
#> 2 format_version 1
#> 3 created_at 2026-10-06 00:37:32 UTC
#> 4 samples 6
#> 5 genes 200
#> 6 zero_effective_length_values 0expected_count.tsv, tpm.tsv, and
effective_length.tsv are gene × sample matrices.
study.rds is what later commands read.
age contains only numbers, so it was treated as numeric;
condition and batch became categorical. If a
numeric-looking column is really a set of groups (for example time
points 0, 6, 24), add --p-categorical time to
data import.
2. Expression summaries
#> Wrote expression summaries to results/summary/
show_table("results/summary/samples.tsv")
#> sample_id total_expected_count detected_genes undetected_genes
#> 1 C1 21880 200 0
#> 2 C2 21970 200 0
#> 3 C3 21960 200 0
#> 4 D1 27225 200 0
#> 5 D2 26585 200 0
#> 6 D3 26765 200 0
#> detected_fraction tpm_q25 tpm_median tpm_q75 tpm_sum
#> 1 1 2619.979 4326.391 6689.095 1e+06
#> 2 1 2636.108 4416.054 6665.577 1e+06
#> 3 1 2628.449 4448.008 6661.793 1e+06
#> 4 1 1477.125 3068.127 5310.146 1e+06
#> 5 1 1511.890 3073.025 5360.227 1e+06
#> 6 1 1479.227 2939.596 5321.024 1e+06
show_table("results/summary/genes.tsv", n = 3)
#> gene_id ensembl_gene_id mean_count median_count mean_tpm median_tpm
#> 1 ENSG00000000001.2 ENSG00000000001 186.125 169.625 12479.56 11436.16
#> 2 ENSG00000000002.3 ENSG00000000002 228.625 212.125 14745.73 13721.42
#> 3 ENSG00000000003.1 ENSG00000000003 271.125 254.625 16827.82 15821.05
#> detected_samples detected_fraction
#> 1 6 1
#> 2 6 1
#> 3 6 1These are expression-level summaries computed from the RSEM tables, not alignment QC. rsemflow does not flag or remove samples; that decision stays with you.
3. Sample structure: PCA and correlation
#> Warning: replacing previous import 'S4Arrays::makeNindexFromArrayViewport' by
#> 'DelayedArray::makeNindexFromArrayViewport' when loading 'SummarizedExperiment'
#> Wrote PCA tables to results/pca/By default PCA uses the variance-stabilized counts (DESeq2 VST) of the 500 most variable genes (here all 200), centred and not scaled.
show_table("results/pca/scores.tsv", cols = c("sample_id", "PC1", "PC2", "PC3"))
#> sample_id PC1 PC2 PC3
#> 1 C1 -15.586657 9.0209984 -4.491250
#> 2 C2 -13.842807 -0.6441269 2.820268
#> 3 C3 -8.224491 -9.2509593 4.632941
#> 4 D1 8.610575 -6.3443824 -7.339276
#> 5 D2 13.626016 0.5887981 -1.590494
#> 6 D3 15.417364 6.6296722 5.967811
show_table("results/pca/variance.tsv", n = 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.9037183loadings.tsv holds the gene loadings and
analysis-info.tsv records the transformation and number of
genes used.
#> Wrote sample correlation matrix to results/correlation.tsv
corr <- show_table("results/correlation.tsv")
corr[, -1] <- round(corr[, -1], 3)
corr
#> sample_id C1 C2 C3 D1 D2 D3
#> 1 C1 1.000 0.583 0.259 0.039 -0.093 -0.182
#> 2 C2 0.583 1.000 0.620 0.236 0.057 -0.075
#> 3 C3 0.259 0.620 1.000 0.471 0.255 0.080
#> 4 D1 0.039 0.236 0.471 1.000 0.812 0.649
#> 5 D2 -0.093 0.057 0.255 0.812 1.000 0.824
#> 6 D3 -0.182 -0.075 0.080 0.649 0.824 1.000Use --p-method spearman or
--p-transform log2-tpm to change the defaults.
4. Normalized expression
#> Wrote expression matrix to results/vst.tsv--p-method accepts vst,
normalized-counts (DESeq2 size factors), and
log2-tpm (log2(TPM + 1)).
5. Differential expression
The samples come from two batches, so batch goes into the model and the Drug-versus-Control effect is estimated within batch:
rsemflow differential deseq2 \
--i-study study/ \
--p-formula '~ batch + condition' \
--p-term condition \
--p-reference condition::Control \
--o-differential results/differential/#> Wrote 1 differential result table(s) to results/differential/
list.files("results/differential")
#> [1] "analysis-info.tsv" "condition_Drug_vs_Control.tsv"
#> [3] "contrasts.tsv" "design-matrix.tsv"
de <- utils::read.delim("results/differential/condition_Drug_vs_Control.tsv")
de <- de[order(de$padj), ]
utils::head(de[, c("gene_id", "contrast", "baseMean", "log2FoldChange", "lfcSE", "stat", "padj")])
#> gene_id contrast baseMean log2FoldChange lfcSE
#> 33 ENSG00000000033.1 Drug vs Control 38.40338 -1.671021 0.2315823
#> 23 ENSG00000000023.3 Drug vs Control 45.12578 -1.774831 0.2524512
#> 22 ENSG00000000022.2 Drug vs Control 33.54292 -1.607113 0.2404466
#> 34 ENSG00000000034.2 Drug vs Control 49.82866 -1.825825 0.3096257
#> 30 ENSG00000000030.1 Drug vs Control 105.94716 -4.876302 0.8843351
#> 32 ENSG00000000032.3 Drug vs Control 26.84058 -1.428867 0.2612781
#> stat padj
#> 33 -7.215667 1.073405e-10
#> 23 -7.030393 2.059522e-10
#> 22 -6.683868 1.551447e-09
#> 34 -5.896878 1.852214e-07
#> 30 -5.514088 1.402374e-06
#> 32 -5.468760 1.510649e-06The first six columns identify the comparison; the rest are the DESeq2 result columns, with unshrunken log2 fold changes. The three other files record exactly what was fitted:
show_table("results/differential/analysis-info.tsv", n = 10)
#> key value
#> 1 formula ~batch + condition
#> 2 term condition
#> 3 test wald
#> 4 reduced
#> 5 alpha 0.05
show_table("results/differential/contrasts.tsv")
#> name term contrast coefficient
#> 1 condition_Drug_vs_Control condition condition::Drug::Control NA
show_table("results/differential/design-matrix.tsv")
#> sample_id (Intercept) batchB conditionDrug
#> 1 C1 1 0 0
#> 2 C2 1 1 0
#> 3 C3 1 0 0
#> 4 D1 1 1 1
#> 5 D2 1 0 1
#> 6 D3 1 1 1For paired, multi-group, numeric, interaction, and time-course
designs, see vignette("designs", package = "rsemflow").
6. Preranked GSEA
GSEA ranks genes by the DESeq2 Wald statistic (stat).
With real human or mouse data, MSigDB collections are one option
away:
rsemflow enrichment gsea \
--i-differential results/differential/condition_Drug_vs_Control.tsv \
--p-species human --p-collection hallmark \
--o-enrichment results/gsea.tsvThe example genes are synthetic and match no MSigDB set, so here we use a custom GMT file: genes 1–20, genes 21–40, and genes 101–150.
ids <- sprintf("ENSG%011d", 1:200)
writeLines(c(
paste(c("GENES_001_020", "synthetic", ids[1:20]), collapse = "\t"),
paste(c("GENES_021_040", "synthetic", ids[21:40]), collapse = "\t"),
paste(c("GENES_101_150", "synthetic", ids[101:150]), collapse = "\t")
), "sets.gmt")rsemflow enrichment gsea \
--i-differential results/differential/condition_Drug_vs_Control.tsv \
--i-genesets sets.gmt \
--p-min-size 10 \
--o-enrichment results/gsea.tsv#> Wrote 3 GSEA results to results/gsea.tsv
gsea <- utils::read.delim("results/gsea.tsv")
gsea$leadingEdge <- substr(gsea$leadingEdge, 1, 40)
gsea[, c("pathway", "size", "NES", "pval", "padj", "leadingEdge")]
#> pathway size NES pval padj
#> 1 GENES_021_040 20 -1.9646165 1.272006e-09 3.816017e-09
#> 2 GENES_001_020 20 2.4456215 5.528771e-09 8.293156e-09
#> 3 GENES_101_150 50 -0.8387051 7.555310e-01 7.555310e-01
#> leadingEdge
#> 1 ENSG00000000033;ENSG00000000023;ENSG0000
#> 2 ENSG00000000001;ENSG00000000012;ENSG0000
#> 3 ENSG00000000104;ENSG00000000136;ENSG0000leadingEdge lists genes separated by ;.
p-values are estimated by sampling, so rsemflow fixes the seed
(--p-seed, default 42) and repeated runs give the same
table.
7. GSVA and ssGSEA scores
rsemflow enrichment gsva \
--i-study study/ \
--i-genesets sets.gmt \
--p-min-size 10 \
--o-scores results/gsva.tsv#> Warning: replacing previous import 'S4Arrays::makeNindexFromArrayViewport' by
#> 'DelayedArray::makeNindexFromArrayViewport' when loading 'HDF5Array'
#> Wrote 3 pathway score rows to results/gsva.tsv
gsva <- utils::read.delim("results/gsva.tsv", check.names = FALSE)
gsva <- gsva[, -(2:4)]
gsva[, -1] <- round(gsva[, -1], 3)
gsva
#> pathway C1 C2 C3 D1 D2 D3
#> 1 GENES_001_020 -0.017 -0.708 -0.810 0.870 0.690 0.200
#> 2 GENES_021_040 0.706 0.822 0.156 -0.234 -0.772 -0.817
#> 3 GENES_101_150 -0.066 -0.083 0.174 -0.173 0.010 0.097Scores are computed per sample from log2(TPM + 1); add
--p-method ssgsea for single-sample GSEA. rsemflow does not
test scores between groups.
8. Gene annotation
Annotation adds gene symbols and biotypes while keeping the original
RSEM gene_id. The most faithful source is the GTF used to
build the RSEM reference:
rsemflow annotation gtf --i-study study/ --i-gtf Homo_sapiens.GRCh38.gtf.gz --o-annotation results/annotation.tsvWithout the GTF, the Bioconductor OrgDb packages map stable Ensembl IDs:
rsemflow annotation ensembl --i-study study/ --p-species human --o-annotation results/annotation.tsvJoin the annotation to any result table by gene_id.
Everything that was written
list.files("results", recursive = TRUE)
#> [1] "correlation.tsv"
#> [2] "differential/analysis-info.tsv"
#> [3] "differential/condition_Drug_vs_Control.tsv"
#> [4] "differential/contrasts.tsv"
#> [5] "differential/design-matrix.tsv"
#> [6] "gsea.tsv"
#> [7] "gsva.tsv"
#> [8] "pca/analysis-info.tsv"
#> [9] "pca/loadings.tsv"
#> [10] "pca/scores.tsv"
#> [11] "pca/variance.tsv"
#> [12] "summary/genes.tsv"
#> [13] "summary/samples.tsv"
#> [14] "vst.tsv"