Genomic association and regulatory-variation analysis

Practical 29: Genomics, systems biology and multi-omics

Learning objectives

Define a counted-allele dosage; fit a restricted eQTL-style mean model with covariates; compare correlated candidate evidence; distinguish association, prediction and mediation.

Biological question and data

Does marker1 associate with the fictional gene01 expression profile after accounting for diet and batch?

The generator supplies a marker1 effect on gene01 log expected expression, as well as diet and batch effects. Animals are unrelated/independent in the model. The supplied exposure is known. The quasi-Poisson working model is not a calibrated GWAS/eQTL pipeline for real related animals.

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. Align the counted allele and expression

The 0/1/2 dosage counts a fictional designated allele. Keep sample IDs and the exposure offset explicit.

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()
meta <- d$metadata
meta$dosage <- d$genotypes[meta$animal,"marker1"]
meta$count <- as.vector(d$counts["gene01",meta$animal])
meta$exposure <- unname(d$exposure[meta$animal])
stopifnot(all(meta$dosage%in%0:2),qr(model.matrix(~diet+batch+dosage,meta))$rank==4)
table(meta$dosage)

 0  1  2 
24 17  7

A dosage cannot be called favourable or causal from its numeric coding.

2. Fit the supplied expression-association model

Use diet, batch, dosage and exposure. Translate the dosage coefficient into expected expression ratio per counted allele.

Show R code 
fit <- glm(count~diet+batch+dosage+offset(log(exposure)),data=meta,family=quasipoisson())
tab <- coef(summary(fit))
tab
                Estimate Std. Error   t value     Pr(>|t|)
(Intercept)    3.5615999 0.10122556 35.184788 7.373035e-34
dietsupplement 0.6993352 0.09064197  7.715357 1.037504e-09
batchB         0.2401195 0.08424766  2.850163 6.625961e-03
dosage         0.2904990 0.05306185  5.474724 1.983275e-06
Show R code 
c(log_effect_per_allele=coef(fit)["dosage"],
 mean_expression_ratio_per_allele=exp(coef(fit)["dosage"]),
 supplied_generator_ratio=exp(.3))
           log_effect_per_allele.dosage mean_expression_ratio_per_allele.dosage 
                               0.290499                                1.337094 
               supplied_generator_ratio 
                               1.349859

The fitted estimate need not equal exp(.3) in one count realization. The uncertainty belongs to a restricted working variance model, not a qualified biological eQTL analysis.

3. Create a correlated candidate and examine ambiguity

A nearly identical marker can tag the same signal. Genotype correlation is not an independent confirmation.

Show R code 
candidate <- meta$dosage
candidate[1] <- (candidate[1]+1)%%3
c(dosage_correlation=cor(meta$dosage,candidate))
dosage_correlation 
         0.9802964
Show R code 
second <- glm(count~diet+batch+candidate+offset(log(exposure)),data=meta,family=quasipoisson())
c(marker1_effect=coef(fit)["dosage"],tagged_candidate_effect=coef(second)["candidate"])
            marker1_effect.dosage tagged_candidate_effect.candidate 
                        0.2904990                         0.2893723
Show R code 
stopifnot(length(candidate)==48)

A correlated candidate can inherit association. Two conditional models here do not calculate a fine-mapping posterior or identify a causal allele.

4. Verify the allele-coding equivalence

A simple linear comparison shows that opposite counted alleles reverse the effect sign while preserving fitted values. This isolated illustration does not replace the expression count model above.

Show R code 
dosage_A <- c(0,0,1,1,2,2)
phenotype <- c(2,3,4,5,6,7)
dosage_B <- 2-dosage_A
fit_A <- lm(phenotype~dosage_A);fit_B <- lm(phenotype~dosage_B)
c(A_slope=coef(fit_A)[2],B_slope=coef(fit_B)[2])
A_slope.dosage_A B_slope.dosage_B 
               2               -2
Show R code 
stopifnot(max(abs(fitted(fit_A)-fitted(fit_B)))<1e-12)

5. State what candidate evidence leaves unresolved

Correlated markers can tag the same signal. A nearest-gene assignment, eQTL association and trait association do not jointly establish expression mediation. Check assembly, alleles, tissue, LD population and variant coverage before integrating results.

Fine-mapping and colocalization are additional model-based analyses, not interpretations automatically available from the toy fit. Related livestock observations require an appropriate dependence model. A predictive model also needs independent evaluation rather than a ranking by association p-value.

Interpretation challenge

A marker has a positive effect in one report and a negative effect in another. What should be checked before declaring a biological contradiction?

Counted allele, strand, assembly, phenotype direction, tissue and model. Opposite allele coding can reverse a sign for the same fitted comparison.

Interpretation and limits

Real studies need allele/reference alignment, genotype and expression QC, relatedness/structure treatment, adequate replication, cis/trans definitions, multiple testing and model diagnostics. A TWAS association is not automatically evidence that expression mediates the variant–trait effect. Fine-mapping requires an explicit candidate/model/LD framework.

Exercises

  1. What is the expected ratio per allele in the generator?
  2. Does a correlated candidate represent another independent biological effect?
  3. What would an expression-mediation claim additionally require?
  1. exp(.3), about 1.35, conditional on the supplied model.
  2. No; association can be tagged through correlated dosage.
  3. Temporal/causal assumptions and evidence addressing shared genetic and environmental drivers, not merely paired association tests.

Continue learning

B13 Association, B06 Transcriptomics and B12 Statistics.

Continue with integrated interpretation. Return to the practical index or the unified map.