## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----eval = F----------------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # BiocManager::install("doseR") ## ----------------------------------------------------------------------------- library(doseR) library(SummarizedExperiment) ## ----------------------------------------------------------------------------- data(hmel.data.doser) ## ----------------------------------------------------------------------------- reps <- c("Male", "Male", "Male", "Female", "Female", "Female") ## ----------------------------------------------------------------------------- annotxn <- data.frame("Chromosome" = factor(hmel.dat$chromosome, levels = 1:21)) annotxn$ZA <- factor(ifelse(hmel.dat$chromosome == 21, "Z", "A"), levels = c("A", "Z")) ## ----------------------------------------------------------------------------- counts <- hmel.dat$readcounts colData <- S4Vectors::DataFrame(Treatment=as.factor(reps), row.names=colnames(hmel.dat$readcounts)) rowData <- S4Vectors::DataFrame(annotation = annotxn, seglens = hmel.dat$trxLength, row.names=rownames(hmel.dat$readcounts) ) se2 <- SummarizedExperiment(assays=list(counts=counts), colData=colData, rowData=rowData) ## ----------------------------------------------------------------------------- SummarizedExperiment::colData(se2)$Libsizes <- getLibsizes3(se2, estimationType = "total") ## ----------------------------------------------------------------------------- SummarizedExperiment::assays(se2)$rpkm <- make_RPKM(se2) # se2 now equals se... data(hmel.se) # loads hmel.dat list # MD5 checksum equivalence: #library(digest) digest::digest(se) == digest::digest(se2) ## ----fig.align="center", fig.width=6------------------------------------------ plotMA.se(se, samplesA ="Male", samplesB = "Female", cex = .2 , pch = 19, col = rgb(0,0,0, .2), xlab = "Log2(Average RPKM)", ylab = "Log2(Male:Female)") ## ----------------------------------------------------------------------------- f_se <- simpleFilter(se, mean_cutoff = 0.01, counts = FALSE) ## ----fig.align="center", fig.width=5------------------------------------------ plotExpr(f_se, groupings = "annotation.ZA", clusterby_grouping = FALSE, col=c("grey80","red","grey80","red"), notch=TRUE, ylab = "Log2(RPKM)") ## ----------------------------------------------------------------------------- se.male <- f_se[, colData(f_se)$Treatment == "Male"] male_ZvA <- generateStats(se.male , groupings = "annotation.ZA", LOG2 = FALSE) male_ZvA$summary # distributional summary statistics male_ZvA$kruskal # htest class output from kruskal.test() lapply(male_ZvA$data, head) # a record of values used for statistics. ## ----fig.align="center"------------------------------------------------------- plotExpr(f_se, groupings = "annotation.Chromosome", col=c(rep("grey80", 20), "red"), notch=TRUE, ylab = "Log2(RPKM)", las = 2, treatment = "Male", clusterby_grouping = TRUE ) ## ----fig.align="center", fig.width=4------------------------------------------ plotRatioBoxes(f_se, groupings = "annotation.ZA", treatment1 = "Male", treatment2 = "Female", outline = FALSE, col = c("grey80", "red"), ylab = "Log(Male:Female)" ) ## ----fig.align="center", fig.width = 5---------------------------------------- plotRatioDensity(f_se, groupings = "annotation.ZA", treatment1 = "Male", treatment2 = "Female", type = "l", xlab = "Log(Male:Female)", ylab = "Density") ## ----fig.align="center", fig.width=10----------------------------------------- par(mfrow = c(1,2)) plotRatioBoxes(f_se, groupings = "annotation.Chromosome", treatment1 = "Male", treatment2 = "Female", outline = FALSE, col=c(rep("grey80", 20), "red"), ylab = "Log(Male:Female)", xlab = "Chromosome" ) plotRatioDensity(f_se, groupings = "annotation.Chromosome", treatment1 = "Male", treatment2 = "Female", type = "l", xlab = "Log(Male:Female)", ylab = "Density", col=c(rep("grey80", 20), "red"), lty = 1) ## ----------------------------------------------------------------------------- za.ratios.test <- test_diffs(f_se, groupings = "annotation.ZA", treatment1 = "Male", treatment2 = "Female", LOG2 = FALSE ) za.ratios.test$summary # summary statistics for each grouping za.ratios.test$kruskal # htest class output from kruskal.test() lapply(za.ratios.test$data, head) # values used for summaries and tests