## ----setup_knitr, include = FALSE, cache = FALSE------------------------------ library("Rmmquant") library("BiocStyle") library("S4Vectors") library("SummarizedExperiment") library("knitr") library("rmarkdown") library("TBX20BamSubset") library("TxDb.Mmusculus.UCSC.mm9.knownGene") library("org.Mm.eg.db") library("DESeq2") opts_chunk$set(message = FALSE, error = FALSE, warning = FALSE, cache = FALSE, fig.width = 5, fig.height = 5) ## ---- first------------------------------------------------------------------- dir <- system.file("extdata", package="Rmmquant", mustWork = TRUE) gtfFile <- file.path(dir, "test.gtf") bamFile <- file.path(dir, "test.bam") output <- RmmquantRun(gtfFile, bamFile) ## ---- matrix------------------------------------------------------------------ assays(output)$counts ## ---- counts------------------------------------------------------------------ assays(output)$counts ## ---- stats------------------------------------------------------------------- colData(output) ## ---- bamfiles---------------------------------------------------------------- bamFiles <- getBamFileList() sampleNames <- names(bamFiles) ## ---- annotation-------------------------------------------------------------- gr <- genes(TxDb.Mmusculus.UCSC.mm9.knownGene, filter=list(tx_chrom="chr19")) ## ---- ensembl----------------------------------------------------------------- ensemblIds <- sapply(as.list(org.Mm.egENSEMBL[mappedkeys(org.Mm.egENSEMBL)]) [mcols(gr)$gene_id], `[[`, 1) gr <- gr[! is.na(names(ensemblIds)), ] names(gr) <- unlist(ensemblIds) ## ---- deseq2------------------------------------------------------------------ output <- RmmquantRun(genomicRanges=gr, readsFiles=bamFiles, sampleNames=sampleNames, sorts=FALSE) coldata <- data.frame(condition=factor(c(rep("control", 3), rep("KO", 3)), levels=c("control", "KO")), row.names=sampleNames) dds <- DESeqDataSetFromMatrix(countData=assays(output)$counts, colData =coldata, design =~ condition) dds <- DESeq(dds) res <- lfcShrink(dds, coef=2) res$padj <- ifelse(is.na(res$padj), 1, res$padj) res[res$padj < 0.05, ] ## ---- session_info------------------------------------------------------------ devtools::session_info()