## ----setup, echo=F-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- knitr::opts_chunk$set(fig.width=6, fig.height=4, cache=TRUE, autodep=TRUE) # To examine objects: # devtools::load_all(".", export_all=F) ; qwraps2::lazyload_cache_dir("vignettes/3_shift_cache/html") options(width=1000) ## ----load, message=F, warning=F------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- library(tidyverse) library(reshape2) library(limma) library(topconfects) library(weitrix) # Produce consistent results set.seed(12345) # 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() ) peaks <- system.file("GSE83162", "peaks.csv.gz", package="weitrix") %>% read_csv() counts <- system.file("GSE83162", "peak_count.csv.gz", package="weitrix") %>% read_csv() %>% column_to_rownames("name") %>% as.matrix() genes <- system.file("GSE83162", "genes.csv.gz", package="weitrix") %>% read_csv() %>% column_to_rownames("name") samples <- data.frame(sample=I(colnames(counts))) %>% extract(sample, c("strain","time"), c("(.+)-(.+)"), remove=FALSE) %>% mutate( strain=factor(strain,unique(strain)), time=factor(time,unique(time))) rownames(samples) <- samples$sample groups <- dplyr::select(peaks, group=gene_name, name=name) # Note the order of this data frame is important ## ----examine_raw---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- samples head(groups, 10) counts[1:10,1:5] ## ----shift---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- wei <- counts_shift(counts, groups) colData(wei) <- cbind(colData(wei), samples) rowData(wei) <- cbind(rowData(wei), genes[match(rownames(wei), rownames(genes)),]) ## ----comp, message=F------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ comp_seq <- weitrix_components_seq(wei, p=10, design=~0) ## ----scree---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- components_seq_screeplot(comp_seq) ## ----exam, fig.height=6--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- comp <- comp_seq[[4]] matrix_long(comp$col, row_info=samples, varnames=c("sample","component")) %>% ggplot(aes(x=time, y=value, color=strain, group=strain)) + geom_hline(yintercept=0) + geom_line() + geom_point(alpha=0.75, size=3) + facet_grid(component ~ .) + labs(title="Sample scores for each component", y="Sample score", x="Time", color="Strain") ## ----calibrate------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ cal_comp <- weitrix_calibrate_trend(wei, comp) metadata(cal_comp)$weitrix$trend_fit ## ----calibrate_show------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- rowData(cal_comp) %>% as.data.frame() %>% ggplot(aes(x=total_reads, y=dispersion_before)) + facet_wrap(~ cut(per_read_var, 9)) + geom_point(size=0.1) + geom_point(aes(y=dispersion_trend), color="red", size=0.1) + scale_x_log10() + scale_y_log10() + labs(y="Dispersion (log scale)", x="Total reads (log scale)") ## ----limma---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- fit_comp <- cal_comp %>% weitrix_elist() %>% lmFit(comp$col) ## ----C1------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- limma_confects(fit_comp, "C1", full=TRUE, fdr=0.05) ## ----C2------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- limma_confects(fit_comp, "C2", full=TRUE, fdr=0.05) ## ----C3------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- limma_confects(fit_comp, "C3", full=TRUE, fdr=0.05) ## ----C4------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- limma_confects(fit_comp, "C4", full=TRUE, fdr=0.05) ## ----examiner, message=F, warning=F, fig.height=3------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- examine <- function(gene_wanted, title) { peak_names <- filter(peaks, gene_name==gene_wanted)$name counts[peak_names,] %>% melt(value.name="reads", varnames=c("peak","sample")) %>% left_join(samples, by="sample") %>% ggplot(aes(x=factor(as.integer(peak)), y=reads)) + facet_grid(strain ~ time) + geom_col() + labs(x="Peak",y="Reads",title=title) } examine("YLR058C", "SHM2, from C1") examine("YLR333C", "RPS25B, from C2") examine("YDR077W", "SED1, from C3") examine("YIL015W", "BAR1, from C4") examine("tK(CUU)M", "tK(CUU)M, from C4")