## ----echo=F--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # To test # devtools::load_all(".",export_all=F) ; rmarkdown::render("vignettes/airway.Rmd") knitr::opts_chunk$set(fig.width=6, fig.height=4) ## ----setup, warning=F, message=F------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ library(weitrix) library(EnsDb.Hsapiens.v86) library(edgeR) library(limma) library(reshape2) library(tidyverse) library(airway) set.seed(1234) # BiocParallel supports multiple backends. # If the default hangs or errors, try others. BiocParallel::register( BiocParallel::SnowParam() ) # The most reliable backed is to use serial processing #BiocParallel::register( BiocParallel::SerialParam() ) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- data("airway") airway ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- counts <- assay(airway,"counts") design <- model.matrix(~ dex + cell, data=colData(airway)) good <- filterByExpr(counts, design=design) table(good) airway_elist <- DGEList(counts[good,]) %>% calcNormFactors() %>% voom(design, plot=TRUE) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- airway_weitrix <- as_weitrix(airway_elist) # Include row and column information colData(airway_weitrix) <- colData(airway) rowData(airway_weitrix) <- mcols(genes(EnsDb.Hsapiens.v86))[rownames(airway_weitrix),c("gene_name","gene_biotype")] airway_weitrix ## ----message=F------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ comp_seq <- weitrix_components_seq(airway_weitrix, p=6) comp_seq ## ----message=F------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ rand_weitrix <- weitrix_randomize(airway_weitrix) rand_comp <- weitrix_components(rand_weitrix, p=1) components_seq_screeplot(comp_seq, rand_comp) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- comp <- comp_seq[[4]] comp$col[,-1] %>% melt(varnames=c("Run","component")) %>% left_join(as.data.frame(colData(airway)), by="Run") %>% ggplot(aes(y=cell, x=value, color=dex)) + geom_vline(xintercept=0) + geom_point(alpha=0.5, size=3) + facet_grid(~ component) + labs(title="Sample scores for each component", x="Sample score", y="Cell line", color="Treatment") comp$row[,-1] %>% melt(varnames=c("name","component")) %>% ggplot(aes(x=comp$row[name,"(Intercept)"], y=value)) + geom_point(cex=0.5, alpha=0.5) + facet_wrap(~ component) + labs(title="Gene loadings for each component vs average log2 RPM", x="average log2 RPM", y="gene loading") ## ----message=F------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ comp_nonvarimax <- weitrix_components(airway_weitrix, p=4, use_varimax=FALSE) comp_nonvarimax$col[,-1] %>% melt(varnames=c("Run","component")) %>% left_join(as.data.frame(colData(airway)), by="Run") %>% ggplot(aes(y=cell, x=value, color=dex)) + geom_vline(xintercept=0) + geom_point(alpha=0.5, size=3) + facet_grid(~ component) + labs(title="Sample scores for each component, no varimax rotation", x="Sample score", y="Cell line", color="Treatment") ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- airway_elist <- weitrix_elist(airway_weitrix) fit <- lmFit(airway_elist, comp$col) %>% eBayes() fit$df.prior fit$s2.prior topTable(fit, "C1") all_top <- topTable(fit, "C1", n=Inf, sort.by="none") plotMD(fit, "C1", status=all_top$adj.P.Val <= 0.01) ## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- library(topconfects) limma_confects(fit, "C1")