Differential-expression designs • rsemflow

rsemflow differential deseq2 fits one DESeq2 model and writes one table per comparison. Three options define the analysis:

Option	Meaning
`--p-formula`	the full model, e.g. `'~ batch + condition'`
`--p-term`	which term’s results to write
`--p-reference COLUMN::LEVEL`	the reference (denominator) level of a factor; repeatable

Quote the formula so the shell does not split it. Each section below shows the metadata, the command, the equivalent R call, and the tables produced.

Two groups

sample_id   condition
C1  Control
C2  Control
C3  Control
D1  Drug
D2  Drug
D3  Drug

rsemflow differential deseq2 --i-study study/ \
  --p-formula '~ condition' --p-term condition \
  --p-reference condition::Control \
  --o-differential results/differential/

differential_deseq2(study, ~ condition, "condition", reference = "condition::Control")

Output: condition_Drug_vs_Control.tsv. Positive log2FoldChange means higher in Drug.

Adjusting for batch

Put the nuisance variable before the variable of interest:

--p-formula '~ batch + condition' --p-term condition --p-reference condition::Control

Batch and condition must not be confounded. If every Control sample was in batch A and every Drug sample in batch B, the two effects cannot be separated and rsemflow stops with “The design matrix is not full rank”.

Paired samples

sample_id   patient condition
P01_pre P01 Pre
P01_post    P01 Post
P02_pre P02 Pre
P02_post    P02 Post

--p-formula '~ patient + condition' --p-term condition --p-reference condition::Pre

Including patient compares Post with Pre within each patient. Output: condition_Post_vs_Pre.tsv.

More than two groups

--p-formula '~ treatment' --p-term treatment --p-reference treatment::Vehicle

Every level is compared with the reference: treatment_DrugA_vs_Vehicle.tsv, treatment_DrugB_vs_Vehicle.tsv. Add other pairs with --p-contrast COLUMN::NUMERATOR::DENOMINATOR (repeatable):

--p-contrast treatment::DrugB::DrugA

differential_deseq2(
  study, ~ treatment, "treatment",
  reference = "treatment::Vehicle",
  contrast = "treatment::DrugB::DrugA"
)

A numeric covariate

When the term is numeric, the table reports the change in log2 expression per unit of the variable:

--p-formula '~ sex + age + response_score' --p-term response_score

Output: response_score.tsv, with contrast = per_unit. Columns are numeric when every value parses as a number; force a type at import with --p-categorical or --p-numeric. For example, sampling times 0, 6, and 24 h analysed as three groups:

rsemflow data import --i-rsem-dir data/ --m-metadata metadata.tsv \
  --p-categorical time --o-study study/

Interaction

Does the drug response differ between genotypes?

--p-formula '~ genotype * treatment' \
--p-term genotype:treatment \
--p-reference genotype::WT --p-reference treatment::Vehicle

One table per interaction coefficient, named after the DESeq2 coefficient (genotypeKO.treatmentDrug.tsv). Its log2FoldChange is the drug effect in KO minus the drug effect in WT.

In an interaction model, --p-term treatment on its own gives the drug effect at the reference genotype (WT). rsemflow warns about this.

The effect within one group (simple effect)

To get the drug effect within KO from the same interaction model, add one --p-contrast and one --p-at:

--p-formula '~ genotype * treatment' \
--p-term treatment \
--p-reference genotype::WT --p-reference treatment::Vehicle \
--p-contrast treatment::Drug::Vehicle \
--p-at genotype::KO

differential_deseq2(
  study, ~ genotype * treatment, "treatment",
  reference = c("genotype::WT", "treatment::Vehicle"),
  contrast = "treatment::Drug::Vehicle",
  at = "genotype::KO"
)

Output: treatment_Drug_vs_Vehicle_at_genotype_KO.tsv, computed as the main effect plus the KO interaction coefficient (recorded in contrasts.tsv). The package tests check that this equals the drug effect obtained by refitting with KO as the reference genotype.

In version 0.1.0, --p-at accepts exactly one --p-at and one --p-contrast, and the contrast denominator must be the reference level.

Any effect of time (likelihood-ratio test)

To ask whether expression changes over time at all, compare a full and a reduced model:

--p-formula '~ patient + time' \
--p-test lrt --p-reduced '~ patient' \
--p-term time

differential_deseq2(
  study, ~ patient + time, "time",
  test = "lrt", reduced = ~ patient
)

Output: time_omnibus.tsv. The reduced formula must be nested in the full one and must drop the term. Because the test has no single direction, log2FoldChange and lfcSE are NA, and the table cannot be used for preranked GSEA.

What every run records

File	Content
`design-matrix.tsv`	the model matrix actually fitted, one row per sample
`contrasts.tsv`	one row per result table: name, term, contrast, DESeq2 coefficient
`analysis-info.tsv`	formula, term, test, reduced formula, alpha

What rsemflow checks, and what it does not

Checked before fitting: model variables exist and have no missing values, the term is part of the formula, reference and contrast levels exist, and the design matrix is full rank.

Not checked: whether you have enough biological replicates, whether the model is appropriate for your experiment, or whether samples should be excluded. Log fold changes are not shrunk.