The benefit of sharing information across genes
John Blischak
2018-08-20
Last updated: 2019-02-26
Checks: 6 0
Knit directory: dc-bioc-limma/analysis/
This reproducible R Markdown analysis was created with workflowr (version 1.2.0.9000). The Report tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(12345) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Untracked files:
Untracked: analysis/table-s1.txt
Untracked: analysis/table-s2.txt
Untracked: code/tb-scratch.R
Untracked: data/counts_per_sample.txt
Untracked: docs/style.css
Untracked: docs/table-s1.txt
Untracked: docs/table-s2.txt
Untracked: factorial-dox.rds
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | 2372aa1 | John Blischak | 2019-01-09 | Build site. |
| html | f440a87 | John Blischak | 2018-08-20 | Build site. |
| html | 4976490 | John Blischak | 2018-08-20 | Build site. |
| Rmd | 1920833 | John Blischak | 2018-08-20 | Refactor first edition of chapter 1 into distinct lessons. |
Introduction
The simulation and visualizations below demonsrate the differences in the results due to limma sharing information across genes to shrink the estimates of the variance.
Setup
library("cowplot")
library("dplyr")
library("ggplot2")
theme_set(theme_classic(base_size = 16))
library("knitr")
opts_chunk$set(fig.width = 10, fig.height = 5, message = FALSE)
library("stringr")
library("tidyr")Simulation
Create some synthetic data for illustrating concepts. The simulated gene expression matrix has 100 genes and 6 samples (3 treatment and 3 control).
set.seed(12345)
create_exp_mat <- function(n1, n2, ng,
alpha_mean, beta_mean, epsilon_sd) {
status <- c(rep(0, n1), rep(1, n2))
ns <- length(status)
status <- matrix(status, nrow = 1)
alpha <- rnorm(ng, mean = alpha_mean, sd = 1)
beta <- matrix(rnorm(ng, mean = beta_mean, sd = 1), ncol = 1)
epsilon <- matrix(rnorm(ng * ns, mean = 0, sd = epsilon_sd),
nrow = ng, ncol = ns)
Yg <- alpha + beta %*% status + epsilon
return(Yg)
}
gexp <- rbind(
# 30 non-DE genes with high variance
create_exp_mat(n1 = 3, n2 = 3, ng = 30, alpha_mean = 10, beta_mean = -1:1, epsilon_sd = 3),
# 30 non-DE genes with low variance
create_exp_mat(n1 = 3, n2 = 3, ng = 30, alpha_mean = 10, beta_mean = -1:1, epsilon_sd = 1),
# 10 upregulated DE genes with low variance
create_exp_mat(n1 = 3, n2 = 3, ng = 10, alpha_mean = 10, beta_mean = 5, epsilon_sd = 1),
# 10 upregulated DE genes with high variance
create_exp_mat(n1 = 3, n2 = 3, ng = 10, alpha_mean = 10, beta_mean = 5, epsilon_sd = 3),
# 10 downregulated DE genes with low variance
create_exp_mat(n1 = 3, n2 = 3, ng = 10, alpha_mean = 10, beta_mean = -5, epsilon_sd = 1),
# 10 downregulated DE genes with high variance
create_exp_mat(n1 = 3, n2 = 3, ng = 10, alpha_mean = 10, beta_mean = -5, epsilon_sd = 3)
)
# Add names for samples
group <- rep(c("con", "treat"), each = ncol(gexp) / 2)
samples <- paste0(group, 1:3)
colnames(gexp) <- samples
# Add names for genes
genes <- sprintf("gene%02d", 1:nrow(gexp))
rownames(gexp) <- genes
heatmap(gexp)
Standard linear model
Find differentially expressed genes using a standard linear model.
lm_beta <- numeric(length = nrow(gexp))
lm_se <- numeric(length = nrow(gexp))
lm_p <- numeric(length = nrow(gexp))
for (i in 1:length(lm_p)) {
mod <- lm(gexp[i, ] ~ group)
result <- summary(mod)
lm_beta[i] <- result$coefficients[2, 1]
lm_se[i] <- result$coefficients[2, 2]
lm_p[i] <- result$coefficients[2, 4]
}
hist(lm_p, xlab = "p-values", main = "Standard linear model")
limma linear model
Find differentially expressed genes using limma.
library("limma")
design <- model.matrix(~group)
colnames(design) <- c("Intercept", "treat")
fit <- lmFit(gexp, design)
head(fit$coefficients) Intercept treat
gene01 11.316083 -2.4577980
gene02 9.833304 2.7130980
gene03 12.653098 -0.2048963
gene04 12.275601 0.2934781
gene05 8.617135 2.3383110
gene06 5.878178 3.9361382fit <- eBayes(fit)
results <- decideTests(fit[, 2])
summary(results) treat
Down 15
NotSig 71
Up 14stats <- topTable(fit, coef = "treat", number = nrow(fit), sort.by = "none")
hist(stats[, "P.Value"], xlab = "p-values", main = "limma linear model")
Comparison
Compare the p-values from lm and limma (both adjusted for multiple testing with the BH FDR).
stats <- cbind(stats,
sd = apply(gexp, 1, sd),
var = apply(gexp, 1, var),
lm_beta, lm_se,
lm_p = p.adjust(lm_p, method = "BH"))
stats$labels_pre <- c(rep("non-DE; high-var", 30),
rep("non-DE; low-var", 30),
rep("DE-up; low-var", 10),
rep("DE-up; high-var", 10),
rep("DE-down; low-var", 10),
rep("DE-down; high-var", 10))
stats$labels <- rep("non-DE", nrow(stats))
stats$labels[stats$adj.P.Val < 0.05 & stats$lm_p < 0.05] <- "DE"
stats$labels[stats$adj.P.Val < 0.05 & stats$lm_p >= 0.05] <- "limma-only"
stats$labels[stats$adj.P.Val >= 0.05 & stats$lm_p < 0.05] <- "lm-only"
table(stats$labels)
DE limma-only lm-only non-DE
22 7 3 68 table(stats$labels, stats$labels_pre)
DE-down; high-var DE-down; low-var DE-up; high-var
DE 0 7 3
limma-only 1 3 1
lm-only 0 0 0
non-DE 9 0 6
DE-up; low-var non-DE; high-var non-DE; low-var
DE 9 0 3
limma-only 1 1 0
lm-only 0 0 3
non-DE 0 29 24stopifnot(stats$logFC == stats$lm_beta)
de <- data.frame(effect_size = stats$lm_beta,
std_dev = stats$sd,
lm = stats$lm_p < 0.05,
limma = stats$adj.P.Val < 0.05)
head(de) effect_size std_dev lm limma
1 -2.4577980 3.750555 FALSE FALSE
2 2.7130980 3.353835 FALSE FALSE
3 -0.2048963 1.715542 FALSE FALSE
4 0.2934781 3.511502 FALSE FALSE
5 2.3383110 1.955378 FALSE FALSE
6 3.9361382 3.326517 FALSE FALSE# View the number of discrepancies
table(de$lm, de$limma)
FALSE TRUE
FALSE 68 7
TRUE 3 22# Plot effect size (y-axis) vs. standard deviation (x-axis)
ggplot(de, aes(x = std_dev, y = effect_size, color = limma)) +
geom_point()
ggplot(stats, aes(x = sd, y = logFC, color = labels)) +
geom_point()
ggplot(stats, aes(x = logFC, y = -log10(P.Value), color = labels)) +
geom_point()
Example genes
Visualize example genes with boxplots. Note that the limma-only gene has higher variance compared to the lm-only gene.
# Find a good example of a DE gene
index <- which(stats$labels_pre == "DE-up; low-var" & stats$labels == "DE")[1]
single_gene <- gexp %>% as.data.frame %>%
slice(index) %>%
gather(key = "group", value = "gene") %>%
mutate(group = str_extract(group, "[a-z]*")) %>%
as.data.frame()
# Find a gene that is DE for both, DE for lm-only, and DE for limma-only
de_not <- de_lm <- which(stats$labels == "non-DE" &
stats$labels_pre == "non-DE; high-var" &
stats$logFC > 0)[1]
de_both <- which(stats$labels == "DE" &
stats$labels_pre == "DE-up; low-var")[1]
de_lm <- which(stats$labels == "lm-only" &
stats$labels_pre == "non-DE; low-var" &
stats$logFC > 0)[1]
de_limma <- which(stats$labels == "limma-only" &
stats$labels_pre == "DE-up; high-var")[1]
compare <- gexp %>%
as.data.frame() %>%
slice(c(de_not, de_both, de_lm, de_limma)) %>%
mutate(type = c("neither", "both", "lm-only", "limma-only")) %>%
gather(key = "group", value = "gene", con1:treat3) %>%
mutate(group = str_extract(group, "[a-z]*")) %>%
as.data.frame()
head(compare) type group gene
1 neither con 6.681872
2 both con 8.555641
3 lm-only con 9.959914
4 limma-only con 11.391149
5 neither con 8.144218
6 both con 7.977472# Plot gene expression (gene; y-axis) vs. group (x-axis)
ggplot(compare, aes(x = group, y = gene)) +
geom_boxplot() +
facet_wrap(~type, nrow = 1)
sessionInfo()R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] limma_3.38.3 tidyr_0.8.2 stringr_1.4.0 knitr_1.21 dplyr_0.8.0.1
[6] cowplot_0.9.4 ggplot2_3.1.0
loaded via a namespace (and not attached):
[1] Rcpp_1.0.0 compiler_3.5.2 pillar_1.3.1
[4] git2r_0.24.0 plyr_1.8.4 workflowr_1.2.0.9000
[7] tools_3.5.2 digest_0.6.18 evaluate_0.13
[10] tibble_2.0.1 gtable_0.2.0 pkgconfig_2.0.2
[13] rlang_0.3.1 yaml_2.2.0 xfun_0.5
[16] withr_2.1.2 fs_1.2.6 rprojroot_1.2
[19] grid_3.5.2 tidyselect_0.2.5 glue_1.3.0
[22] R6_2.4.0 rmarkdown_1.11 purrr_0.3.0
[25] magrittr_1.5 whisker_0.3-2 backports_1.1.3
[28] scales_1.0.0 htmltools_0.3.6 assertthat_0.2.0
[31] colorspace_1.4-0 labeling_0.3 stringi_1.3.1
[34] lazyeval_0.2.1 munsell_0.5.0 crayon_1.3.4