Dream_and_contrasts_demo.R

# Gabriel Hoffman
# October 19, 2020
#
# Examples illustrating using contrasts with dream

library(variancePartition)
register(SnowParam(4))
data(varPartData)

# Since Tissue is a factor, extract the levels
# The first level will be used as the baseline
# If you want to use a different level as your baseline,
# you can use factor() or relevel()
levels(info$Tissue)
#> [1] "A" "B" "C"

# Example 1: Implied baseline
#----------------------------

# specify formula
form = ~ Tissue + (1|Individual)

# run dream
fit <- dream( geneExpr[1:4,], form, info)

# Look at the coefficient estimates:
# TissueA is not explicitly estimated.
# There is an (Intercept) term, TissueB (relative to baseline (i.e. TissueA))
# and TissueC (relative to baseline)
head(coef(fit))
#>       (Intercept)   TissueB    TissueC
#> gene1  -9.8813655 -1.216217 -0.4601605
#> gene2   0.3256370 -2.144469 -2.1305913
#> gene3   0.7924753 -1.425857 -0.3404834
#> gene4  -3.1557621 -1.839594 -2.2436210

# Note for gene1, the effect of TissueB compared to TissueA (implicitly) is -1.216217

# If you want to compare expression of Tissues A and B, 
# you only need to test the coefficient for TissueB, 
# since the comparison to baseline (TissueA) is already built in
topTable(fit, coef="TissueB")
#>           logFC     AveExpr         t      P.Value    adj.P.Val     z.std
#> gene2 -2.144469  -1.1281610 -7.119066 5.594229e-10 2.237691e-09 -6.201458
#> gene4 -1.839594  -4.5359748 -6.661737 3.947194e-09 7.894389e-09 -5.886391
#> gene3 -1.425857   0.1702122 -5.106590 2.449325e-06 3.265767e-06 -4.712303
#> gene1 -1.216217 -10.4664549 -4.092563 1.069817e-04 1.069817e-04 -3.874186

# Note that (Intercept) does **not** represent baseline,
# it is the global mean after removing all other effects,
# so there is almost never a reason to include (Intercept) in a contrast


# Example 2: Explicit baseline
#----------------------------

# specify formula, omitting intercept
form = ~ 0 + Tissue + (1|Individual)

# run dream
fit <- dream( geneExpr[1:4,], form, info)

# There is no (Intercept) here, so all 3 Tissue levels are estimated explicitly
# Note the difference compared to Example 1:
# Here TissueB is the estimated expression of each gene in Tissue B (after accounting for other variables in the formula)
# The comparison to baseline is not included. 
# So in order to test the difference between Tissues A and B, we need to explicitly specify the contrasts
head(coef(fit))
#>          TissueA     TissueB     TissueC
#> gene1 -9.8813655 -11.0975823 -10.3415260
#> gene2  0.3256370  -1.8188318  -1.8049544
#> gene3  0.7924753  -0.6333813   0.4519919
#> gene4 -3.1557621  -4.9953563  -5.3993831

# Explicitly specify contrasts using desired baseline
# This evaluates TissueB - TissueA
L1 = getContrast( geneExpr[1:4,], form, info, c('TissueB', 'TissueA') )

# evaluate model with contrasts
fit <- dream( geneExpr[1:4,], form, info, L=L1)

# the model now evalutes each coefficient and the custom contrasts
head(coef(fit))
#>              L1    TissueA     TissueB     TissueC
#> gene1 -1.216217 -9.8813655 -11.0975823 -10.3415260
#> gene2 -2.144469  0.3256370  -1.8188318  -1.8049544
#> gene3 -1.425857  0.7924753  -0.6333813   0.4519919
#> gene4 -1.839594 -3.1557621  -4.9953563  -5.3993831

# Note that for gene1, the effect of L1 is -1.216217
# This is exactly the same as the model in Example 1
# As expected, we get the same result

# Now perform hypothesis testing on L1
topTable(fit, coef="L1")
#>           logFC     AveExpr         t      P.Value    adj.P.Val     z.std
#> gene2 -2.144469  -1.1281610 -7.119066 5.594229e-10 2.237691e-09 -6.201458
#> gene4 -1.839594  -4.5359748 -6.661737 3.947194e-09 7.894389e-09 -5.886391
#> gene3 -1.425857   0.1702122 -5.106590 2.449325e-06 3.265767e-06 -4.712303
#> gene1 -1.216217 -10.4664549 -4.092563 1.069817e-04 1.069817e-04 -3.874186