## ----style, echo = FALSE, results = 'asis'-------------------------------------------------------- options(width=100) knitr::opts_chunk$set( eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")), cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE"))) ## ----setup, echo=FALSE---------------------------------------------------------------------------- suppressPackageStartupMessages({ library(ENAR2016) library(ggplot2) library(airway) library(DESeq2) library(org.Hs.eg.db) library(TxDb.Hsapiens.UCSC.hg19.knownGene) ## library(BSgenome.Hsapiens.UCSC.hg19) library(AnnotationHub) }) ## ------------------------------------------------------------------------------------------------- rnorm(10) # sample 10 random normal deviates ## ----vectors-------------------------------------------------------------------------------------- x <- rnorm(1000) y <- x + rnorm(1000) plot(y ~ x) ## ----classes-------------------------------------------------------------------------------------- df <- data.frame(Y = y, X = x) fit <- lm(Y ~ X, df) plot(Y ~ X, df) abline(fit, col="red", lwd=2) anova(fit) ## ----ggplot2, warning=FALSE----------------------------------------------------------------------- library(ggplot2) ggplot(df, aes(x=X, y=Y)) + geom_point() + geom_smooth(color="blue") + geom_smooth(method="lm", color="red") ## ------------------------------------------------------------------------------------------------- library(airway) data("airway") se <- airway ## ------------------------------------------------------------------------------------------------- head(assay(se)) ## ------------------------------------------------------------------------------------------------- colData(se) design = ~ cell + dex se$dex <- relevel(se$dex, "untrt") # 'untrt' as reference level ## ------------------------------------------------------------------------------------------------- rowRanges(se) ## ------------------------------------------------------------------------------------------------- library(DESeq2) dds <- DESeqDataSet(se, design = ~ cell + dex) ## ------------------------------------------------------------------------------------------------- dds <- DESeq(dds) ## ------------------------------------------------------------------------------------------------- res <- results(dds) head(res) mcols(res) # what does each column mean? head(res[order(res$padj),]) # 'top table' by p-adj ## ------------------------------------------------------------------------------------------------- library(org.Hs.eg.db) ttbl <- head(res[order(res$padj),]) # 'top table' by p-adj (ensid <- rownames(ttbl)) mapIds(org.Hs.eg.db, ensid, "SYMBOL", "ENSEMBL") select(org.Hs.eg.db, ensid, c("SYMBOL", "GENENAME"), "ENSEMBL") columns(org.Hs.eg.db) ## ------------------------------------------------------------------------------------------------- ## gene models, organized by entrez identifier library(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene ## from Ensembl to Entrez identifiers egid <- mapIds(org.Hs.eg.db, ensid, "ENTREZID", "ENSEMBL") transcriptsBy(txdb, "gene")[egid] ## ---- eval=FALSE---------------------------------------------------------------------------------- # library(AnnotationHub) # hub = AnnotationHub() # query(hub, c("ensembl", "gtf", "release-83")) ## ------------------------------------------------------------------------------------------------- gr <- GRanges("chr1", IRanges(c(1000, 2000, 3000), width=100), strand=c("+", "-", "*"), score=c(10, 20, 30)) gr ## ------------------------------------------------------------------------------------------------- library(airway) data(airway) airway ## ------------------------------------------------------------------------------------------------- result <- integer() for (i in 1:10) result[i] = sqrt(i) ## ------------------------------------------------------------------------------------------------- result <- integer(10) for (i in 1:10) result[[i]] = sqrt(i) ### same as the much more compact, expressive, and robust... result <- sapply(1:10, sqrt) ## ------------------------------------------------------------------------------------------------- result <- sqrt(1:10) ## ----bioc-parallel-------------------------------------------------------------------------------- library(BiocParallel) system.time({ res <- lapply(1:8, function(i) { Sys.sleep(1); i }) }) system.time({ res <- bplapply(1:8, function(i) { Sys.sleep(1); i }) }) ## ---- eval=FALSE---------------------------------------------------------------------------------- # X <- list(1, "2", 3) # res <- bptry(bplapply(X, sqrt)) # X <- list(1, 2, 3) # res <- bplapply(X, sqrt, BPREDO=res) ## ----r-code--------------------------------------------------------------------------------------- x <- rnorm(1000) plot(x) ## ----sessionInfo---------------------------------------------------------------------------------- sessionInfo()