Practical 9: Genomics, systems biology and multi-omics
Learning objectives
Align counts and metadata; distinguish known exposure from estimated normalization; fit a restricted count model and interpret log2 effects, dispersion and multiplicity.
Biological question and data
Does expression differ between the two diets in a hypothetical 48-animal experiment?
The fixture simulates negative-binomial counts for 24 genes, with independent animals, crossed diet/batch and known multiplicative exposure. The first four genes have a supplied diet effect. Gene01 also has a genotype effect. A quasi-Poisson fit below is an intentionally restricted mean-model illustration; its single variance multiplier per gene is not the generating negative-binomial variance function.
No biological download is required. Numerical inputs are hypothetical, not observed animal data. The runnable examples use base R. Count and abundance units are specified where used.
1. Load and align the count panel
Count rows are genes and columns are animals. Metadata order is verified before modelling.
The +1 transformation supports a plot; it is not the response used in the count model. Real RNA-seq normalization must consider composition and cannot assume this known exposure.
3. Fit the mean model and report uncertainty
Fit each gene with diet, batch and a log-exposure offset. The diet coefficient uses the natural-log link; divide it and its standard error by log(2) to obtain log2 units.
The intervals and p-values are conditional on this restricted working model. Genotype variation for gene01 remains unmodelled here. A fitted association and an approximate p-value are not a causal-mechanism claim.
4. Ask what the design can estimate
Check the model matrix before interpreting a coefficient. In this hypothetical design, both diets occur in each batch. Construct a deliberately confounded alternative and compare its rank.
Show R code
X <-model.matrix(~diet+batch,meta)confounded <- metaconfounded$batch <-factor(ifelse(confounded$diet=="control","A","B"))X_bad <-model.matrix(~diet+batch,confounded)data.frame(design=c("crossed","confounded"),columns=c(ncol(X),ncol(X_bad)),rank=c(qr(X)$rank,qr(X_bad)$rank))
The second design cannot separate diet from batch. Do not solve that biological ambiguity by dropping batch and then calling the remaining coefficient a treatment effect. A simpler model can produce a number while the ambiguity persists.
5. Explore samples on a descriptive scale
PCA uses animals as observations and genes as variables. Here we use log2 exposure-adjusted counts with a pseudocount of one and centre each gene without scaling it to unit variance. This is a teaching transformation rather than a fitted variance-stabilizing transformation.
Describe what you see before inspecting p-values. Do the diets overlap? Does batch appear related to either axis? Neither clustering nor separation provides a test of the fitted diet contrast. Known exposure adjustment does not remove the simulated batch effect.
6. Inspect a fitted gene
Fit gene01 again and examine its fitted means and Pearson residuals. A residual is a discrepancy between an observation and its fitted expectation. We plot residuals rather than silently excluding an unusual animal.
ord <-order(abs(pearson),decreasing=TRUE)head(data.frame(animal=meta$animal[ord],observed=counts["gene01",ord],fitted=fit_mean[ord],residual=pearson[ord]),4)
Look for systematic changes in residual spread or isolated discrepancies. The quasi-Poisson variance is proportional to the mean, while the generator also contains a squared-mean term. A visually acceptable plot does not establish correct p-value calibration. Gene01 additionally has an omitted genotype effect; the model should be assessed against the biological question rather than only the plot.
For comparison, fit a Poisson model to the same mean structure. Its standard error assumes less biological variation. Both fits use the same observations and mean model, so this comparison isolates the variance assumption.
The ratio of these standard errors is the square root of the estimated quasi-Poisson variance multiplier. It need not exceed one for every conceivable dataset; compare the displayed values for this realization.
7. Report effects alongside uncertainty
Plot all genes, maintaining the feature order. These are approximate pointwise intervals from the restricted working model. They are not simultaneous intervals across all 24 genes.
Show R code
plot(result$log2_effect,seq_len(nrow(result)),pch=16,xlim=range(c(result$lower,result$upper)),yaxt="n",xlab="Supplement versus control: fitted log2 abundance ratio",ylab="Gene")segments(result$lower,seq_len(nrow(result)), result$upper,seq_len(nrow(result)))axis(2,at=seq_len(nrow(result)),labels=rownames(result),las=2,cex.axis=.65)abline(v=0,lty=2,col="grey50")
A negative effect denotes lower fitted abundance under supplement; a positive effect denotes higher abundance. An interval crossing zero on the log2 scale crosses one on the ratio scale. Failure to pass a threshold does not prove identical expression.
For this fixture we know the generating diet effects, which is unusual in real data. Compare estimates with that supplied truth without judging the method from one simulation alone.
Show R code
truth <-c(rep(.65/log(2),4),rep(0,20))comparison <-data.frame(gene=rownames(result),truth=truth,estimate=result$log2_effect,selected=result$BH<.05)comparison
Sampling variation, the omitted genotype term and the working variance model all matter. A method-validation study would require repeated simulations spanning relevant designs, with bias, coverage and error rates evaluated explicitly. Counting discoveries in this single panel is an interpretation exercise.
8. Write the analysis record
For your report, state the experimental unit, sample numbers, tested features, supplied exposure, design formula, contrast direction, variance model and multiplicity procedure. Report one effect with its interval and distinguish the fitted tissue-level comparison from a mechanism claim. Record sessionInfo() when reproducing the practical.
A full real-data extension should retain this reporting discipline while replacing the simplified analysis with an appropriate documented RNA-seq workflow. Keep raw-count input, descriptive transformations and normalized displays conceptually separate.
Interpretation and limits
This practical teaches design, exposure, contrasts and uncertainty; it is not a production differential-expression pipeline. The DESeq2 workflow adds composition-aware normalization and a gene-level negative-binomial dispersion framework. Gene filtering, relatedness, batch design and larger-scale validation need explicit decisions for real data.
Exercises
What supplied log2 effect corresponds to .65 on the natural-log scale?
Why is Poisson variance too restrictive for these generated counts?
Change the seed and explain why estimates and adjusted p-values change.
What would a complete treatment–batch confounding problem prevent?
Why might normalized counts still show batch separation in PCA?
Convert the gene01 estimate and interval to abundance ratios and write a result sentence.
Does selecting a gene at BH < .05 establish a 95% probability that its diet effect is real?
Why should an unusual residual trigger inspection of records before exclusion?
Suggested answers
.65/log(2), about .94; a finite-sample estimate need not equal it.
The generating variance is mu + mu squared/12, greater than mu.
Sampling realization changes, while the generating parameters remain the same.
Separate estimation of diet and batch effects.
Exposure adjustment changes the count scale; it does not remove systematic batch differences.
Use 2^estimate, 2^lower and 2^upper; describe the supplement-versus-control tissue-level comparison conditional on batch and the working model.
No. The BH procedure addresses false discovery rate across a family of tests under its assumptions, not an individual posterior probability.
An unusual observation may be biological, technical or a model discrepancy. Exclusion needs an independently justified reason and a sensitivity assessment.
