## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----------------------------------------------------------------------------- #| label: load-packages #| include: false library(diffwrap) ## ----------------------------------------------------------------------------- #| label: example-data set.seed(20240721) n_genes <- 400L; n_true_de <- 60L samples <- sprintf("S%02d", 1:8) group <- rep(c("control", "treated"), each = 4) subject <- rep(c("P1", "P2", "P3", "P4"), times = 2) # each subject: one of each group ## baseline mean expression per gene, with a low-expression tail base_mu <- exp(rnorm(n_genes, mean = 7.5, sd = 1.5)) base_mu[(n_genes - 39):n_genes] <- runif(40, 0, 14) ## true log2 fold changes: first 60 genes differentially expressed, half up half down lfc <- numeric(n_genes) lfc[seq_len(n_true_de)] <- rep(c(1, -1), length.out = n_true_de) * runif(n_true_de, 1.2, 3) ## per-subject offset, so the paired/blocked analyses have a real effect to remove subj_off <- stats::setNames(runif(4, 0.85, 1.18), c("P1", "P2", "P3", "P4")) counts_mat <- sapply(seq_along(samples), function(j) { mu <- base_mu * subj_off[subject[j]] if (group[j] == "treated") mu <- mu * 2^lfc rnbinom(n_genes, mu = mu, size = 1 / 0.15) # negative binomial, dispersion 0.15 }) dimnames(counts_mat) <- list(sprintf("ENSG%011d", seq_len(n_genes)), samples) ## the htseq-count summary rows that filtering is expected to strip special <- c("__no_feature", "__ambiguous", "__too_low_aQual", "__not_aligned", "__alignment_not_unique") counts_mat <- rbind(counts_mat, matrix(sample(20000:90000, length(special) * 8, replace = TRUE), nrow = length(special), dimnames = list(special, samples))) samp_info_raw <- data.frame(SampleName = samples, Group = group, Subject = subject, stringsAsFactors = FALSE) dim(counts_mat) head(counts_mat[, 1:4]) samp_info_raw ## ----------------------------------------------------------------------------- #| label: read-from-file counts_file <- file.path(tempdir(), "example_counts.tsv") write.table(counts_mat, counts_file, sep = "\t", quote = FALSE, col.names = NA) samp.info <- diff_expr_get_samp_info(samp_info_raw, samples = "SampleName", groups = "Group") counts <- diff_expr_read_counts(counts_file, samp.info) dim(counts) # the htseq-count summary rows are still present here head(counts[, 1:4]) ## ----------------------------------------------------------------------------- #| label: full-run #| eval: false # out.dir <- file.path(tempdir(), "diffwrap_example") # dir.create(out.dir, showWarnings = FALSE) # # res <- diffExpr(expr.dat = counts_file, # samp.info = samp_info_raw, # samples = "SampleName", # groups = "Group", # control = "control", # analysis.name = "demo", # out.dir = out.dir, # enr.do = FALSE) ## ----------------------------------------------------------------------------- #| label: paired-run #| eval: false # # 'Subject' pairs each control sample with a treated sample from the same subject # res_paired <- diffExpr(expr.dat = counts_file, # samp.info = samp_info_raw, # samples = "SampleName", # groups = "Group", # pairs = "Subject", # fixed effect, intercept design # control = "control", # analysis.name = "demo_paired", # out.dir = out.dir, # enr.do = FALSE) # # # the same column used as a correlation block instead (forces voom) # res_blocked <- diffExpr(expr.dat = counts_file, # samp.info = samp_info_raw, # samples = "SampleName", # groups = "Group", # pairs = "Subject", # 'pairs' is the blocking variable # block = TRUE, # duplicateCorrelation # control = "control", # analysis.name = "demo_blocked", # out.dir = out.dir, # enr.do = FALSE) ## ----------------------------------------------------------------------------- #| label: manual-filter counts_f <- diff_expr_filter_counts(counts, samp.info, strict = TRUE) dim(counts_f) # low-expression genes and the __-rows have gone ## ----------------------------------------------------------------------------- #| label: manual-model groups <- stats::relevel(samp.info$Groups, ref = "control") d <- edgeR::DGEList(counts = counts_f, group = groups) d <- edgeR::calcNormFactors(d) design <- diff_expr_make_design(samp.info = samp.info, groups = groups) contrasts <- diff_expr_make_contrasts(design = design, groups = groups) design contrasts ## ----------------------------------------------------------------------------- #| label: manual-fit fit.l <- diff_expr_fit(counts = counts_f, d = d, design = design, do.voom = FALSE, quasi.likelihood = TRUE) de <- edgeR::glmQLFTest(fit.l$fit, contrast = contrasts[, 1]) head(edgeR::topTags(de, n = 10)) ## ----------------------------------------------------------------------------- #| label: check-truth tt <- as.data.frame(edgeR::topTags(de, n = Inf)) top50 <- rownames(tt)[seq_len(50)] n_true <- sum(as.integer(sub("^ENSG0*", "", top50)) <= 60) cat("Of the 50 top-ranked genes,", n_true, "are truly differentially expressed\n") ## ----------------------------------------------------------------------------- #| label: reshape-table ## identical in structure to res$contrasts[["treated-control"]] from a full run de_tab <- merge(edgeR::cpm(fit.l$d2, log = TRUE), tt, by = "row.names") names(de_tab)[1] <- "ID" de_tab$gene_symbol <- de_tab$ID # offline demo; a real run already carries symbols ## ----------------------------------------------------------------------------- #| label: reshape-heatmap #| fig-width: 7 #| fig-height: 8 grDevices::pdf(tempfile(fileext = ".pdf")) # swallow the companion correlograms hm <- pheatmap_plots(de_tab, id = "ID", samp.info = samp.info, samples = "SampleNames", groups = groups, fdr.thr = 0.01, topn = 30, color.blind.pal = "RdBu") # default palette is "PuOr" invisible(grDevices::dev.off()) hm$fdr$regular # the FDR-filtered heatmap, new settings ## ----------------------------------------------------------------------------- #| label: reshape-volcano #| eval: false # v <- diff_expr_volcano_plot(de_tab, id = "ID", base.size = 11) # 16 suits the 15-inch PDF # v$FDR ## ----------------------------------------------------------------------------- #| label: logging #| eval: false # res <- diffExpr(..., out.dir = out.dir, verbose = FALSE) # quiet console # res <- diffExpr(..., out.dir = out.dir, verbose = "all") # everything echoed # res <- diffExpr(..., out.dir = out.dir, # log.file = file.path(out.dir, "my_run.log")) # custom log location ## ----------------------------------------------------------------------------- #| label: enrichment #| eval: false # res <- diffExpr(expr.dat = counts_file, # samp.info = samp_info_raw, # samples = "SampleName", # groups = "Group", # control = "control", # analysis.name = "demo_annotated", # out.dir = out.dir, # biom.use = TRUE, # biom.data.set = "hsapiens_gene_ensembl", # enr.do = TRUE, # enr.methods = c("clusterProfilerGO", "gProfileR")) ## ----------------------------------------------------------------------------- #| label: session-info sessionInfo()