Multi-omics associations with an animal-production trait

Practical 22: Genomics, systems biology and multi-omics

Learning objectives

Estimate conditional feature associations; compare marginal and adjusted patterns; report effect scale and multiple testing; distinguish prediction from mechanism.

Biological question and data

Which hypothetical molecular features are associated with the supplied animal response after accounting for diet and batch?

Animals are independent in this fixture. Response units are arbitrary; features are already log2 abundance. The shared latent factor generates association without specifying measured features as causes. Linear-model inference remains conditional on the selected covariates and residual assumptions.

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. Prepare matched measurements and the design

Keep animal rows aligned and examine the crossed design.

Show R code 
helper_dir <- if (dir.exists("R")) "R" else "../R"
stopifnot(file.exists(file.path(helper_dir,"teaching_data.R")))
source(file.path(helper_dir,"teaching_data.R"))
d <- teaching_data()
X <- cbind(t(d$protein),t(d$metabolite));meta <- d$metadata
stopifnot(identical(rownames(X),meta$animal))
xtabs(~diet+batch,meta)
            batch
diet          A  B
  control    12 12
  supplement 12 12
Show R code 
summary(meta$response)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  8.906  10.449  11.580  11.753  12.859  16.172

The target is a feature association conditional on measured covariates, not a treatment effect or breeding-value estimate.

2. Fit a prespecified conditional model per feature

For each feature use response ~ feature + diet + batch. Report the association per one log2 unit of the feature.

Show R code 
analyse <- function(feature) {
 z <- meta;z$feature <- feature
 fit <- lm(response~feature+diet+batch,z)
 tab <- coef(summary(fit))
 c(effect=tab["feature","Estimate"],SE=tab["feature","Std. Error"],
 p=tab["feature","Pr(>|t|)"])
}
result <- as.data.frame(t(apply(X,2,analyse)))
result$BH <- p.adjust(result$p,"BH")
result[order(result$p),]
              effect        SE            p           BH
protein07 0.96467712 0.1785370 2.520293e-06 0.0000352841
protein08 0.77904941 0.1697253 3.685172e-05 0.0001774899
feature05 0.93931832 0.2050779 3.803354e-05 0.0001774899
feature06 0.73810685 0.1810339 1.880842e-04 0.0006094365
protein03 1.53838147 0.3817141 2.176559e-04 0.0006094365
protein06 0.93020081 0.2360360 2.867551e-04 0.0006116257
protein05 1.03244542 0.2633799 3.058129e-04 0.0006116257
protein04 1.02969041 0.2711056 4.435272e-04 0.0007761726
feature04 0.78667260 0.2141877 6.467806e-04 0.0010061032
feature02 0.84730914 0.2449126 1.213522e-03 0.0016989309
feature03 0.75370169 0.2603853 5.888535e-03 0.0074944996
protein02 0.92341706 0.4475763 4.502963e-02 0.0525345672
protein01 0.44569039 0.4257286 3.008680e-01 0.3240117042
feature01 0.05167415 0.2879219 8.583904e-01 0.8583904429
Show R code 
stopifnot(all(is.finite(as.matrix(result))),all(result$BH>=result$p-1e-12))

The adjusted p-values concern the supplied set of tests under the model and applicable dependence conditions. They do not identify causal biomarkers.

3. Compare an unadjusted and conditional relationship

Use the same feature and animals for both fits.

Show R code 
z <- meta;z$feature <- X[,1]
a <- lm(response~feature,z)
b <- lm(response~feature+diet+batch,z)
c(marginal=coef(a)["feature"],conditional=coef(b)["feature"])
   marginal.feature conditional.feature 
          0.5791273           0.4456904
Show R code 
plot(z$feature,z$response,pch=19,col=as.integer(z$diet),
 xlab="Feature log2 abundance",ylab="Response (arbitrary units)")

Adjustment changes the estimand; it is not an automatic guarantee that every nuisance variable has been removed.

4. Check how units change a coefficient

Rescaling a predictor changes its slope and standard error in inverse proportion, while fitted values and its test statistic remain equivalent in this simple model.

Show R code 
set.seed(2202)
x_unit <- rnorm(40);response_unit <- 1+.6*x_unit+rnorm(40)
a_unit <- lm(response_unit~x_unit)
x_ten <- 10*x_unit;b_unit <- lm(response_unit~x_ten)
c(original_slope=coef(a_unit)[2],rescaled_slope=coef(b_unit)[2])
original_slope.x_unit  rescaled_slope.x_ten 
           0.78490807            0.07849081
Show R code 
stopifnot(max(abs(fitted(a_unit)-fitted(b_unit)))<1e-12)

5. Specify adjustment before interpreting direction

A conditional coefficient compares animals at the same values of included covariates under the model. Adding a covariate can change the target as well as the estimate. A variable affected by exposure should not automatically be called a confounder.

Report feature units, tissue/time, sample subset, covariates, uncertainty and multiplicity. Check influential observations and residual patterns before interpreting a ranking. A feature associated after adjustment is not automatically a validated predictor, and predictive value is not proof of a causal role.

Interpretation challenge

Two feature slopes have different magnitudes because one predictor is in log2 units and another in standardized units. Can their raw magnitudes be ranked as biological importance?

No. Slopes depend on scale and the conditional model. Explain a meaningful contrast or standardized comparison and retain uncertainty and measurement context.

Interpretation and limits

Marker relatedness, repeated animals, nonlinear associations, measurement error and selection require suitable analysis in real data. A feature selected on all outcomes cannot then be evaluated on those same animals as an independently validated predictor.

Exercises

  1. What units does the fitted feature coefficient have?
  2. Why can a feature predict without being a cause?
  3. Where should feature selection occur in a prediction workflow?
  1. Arbitrary response units per one log2 unit of the supplied abundance.
  2. It may track a shared biological or technical driver.
  3. Within training-only validation; then evaluate the complete selected pipeline on held-out animals.

Continue learning

B12 Statistics, B13 Association and B18 Synthesis.

Continue with prediction or biomarker validation. Return to the practical index or the unified map.