## ----setup, echo=FALSE--------------------------------------------------- suppressPackageStartupMessages({ library(DESeq2) library(org.Hs.eg.db) }) colDataFile <- system.file(package="Sydney2016", "extdata", "airway-colData.tab") assayFile <- system.file(package="Sydney2016", "extdata", "airway-assay.tab") ## ----dependencies, eval=FALSE-------------------------------------------- # source("https://bioconductor.org/biocLite.R") # biocLite(c("DESeq2", "org.Hs.eg.db")) ## ----Sydney2016, eval=FALSE---------------------------------------------- # biocLite("Bioconductor/Syndey2016") ## ----Sydney2016-filepaths, eval=FALSE------------------------------------ # system.file(package="Sydney2016", "extdata") ## ---- eval=FALSE--------------------------------------------------------- # colDataFile <- file.chooose() # find 'airway-colData.tab' ## ------------------------------------------------------------------------ colData <- read.table(colDataFile) colData ## ---- eval=FALSE--------------------------------------------------------- # assayFile <- file.chooose() # find 'airway-assay.tab' ## ------------------------------------------------------------------------ assay <- read.table(assayFile) head(assay) ## ------------------------------------------------------------------------ colSums(assay) ## ------------------------------------------------------------------------ plot(density(rowMeans(asinh(assay)))) ## ------------------------------------------------------------------------ d <- dist(t(asinh(assay))) plot(cmdscale(d), pch=19, cex=2) plot(cmdscale(d), pch=19, cex=2, col=colData$cell) plot(cmdscale(d), pch=19, cex=2, col=colData$dex) ## ------------------------------------------------------------------------ library(DESeq2) dds <- DESeqDataSetFromMatrix(assay, colData, ~ cell + dex) ## ------------------------------------------------------------------------ dds <- DESeq(dds) dds ## ------------------------------------------------------------------------ result <- results(dds) result ridx <- head(order(result$padj), 10) top = result[ridx,] top ## ------------------------------------------------------------------------ library(org.Hs.eg.db) top$Symbol <- mapIds(org.Hs.eg.db, rownames(top), "SYMBOL", "ENSEMBL") top ## ------------------------------------------------------------------------ sessionInfo()