inst/doc/derfinder-quickstart.R

In derfinder: Annotation-agnostic differential expression analysis of RNA-seq data at base-pair resolution via the DER Finder approach

## ----'installDer', eval = FALSE-----------------------------------------------
#  if (!requireNamespace("BiocManager", quietly = TRUE)) {
#        install.packages("BiocManager")
#    }
#  
#  BiocManager::install("derfinder")
#  
#  ## Check that you have a valid Bioconductor installation
#  BiocManager::valid()

## ----'citation'---------------------------------------------------------------
## Citation info
citation("derfinder")

## ----vignetteSetup, echo=FALSE, message=FALSE, warning = FALSE----------------
## Track time spent on making the vignette
startTime <- Sys.time()

## Bib setup
library("RefManageR")

## Write bibliography information
bib <- c(
    derfinder = citation("derfinder")[1],
    BiocStyle = citation("BiocStyle"),
    knitr = citation("knitr")[3],
    rmarkdown = citation("rmarkdown")[1],
    brainspan = RefManageR::BibEntry(
        bibtype = "Unpublished",
        key = "brainspan",
        title = "Atlas of the Developing Human Brain [Internet]. Funded by ARRA Awards 1RC2MH089921-01, 1RC2MH090047-01, and 1RC2MH089929-01.",
        author = "BrainSpan", year = 2011, url = "http://www.brainspan.org/"
    ),
    originalder = citation("derfinder")[2],
    R = citation(),
    IRanges = citation("IRanges"),
    sessioninfo = citation("sessioninfo"),
    testthat = citation("testthat"),
    GenomeInfoDb = RefManageR::BibEntry(
        bibtype = "manual",
        key = "GenomeInfoDb",
        author = "Sonali Arora and Martin Morgan and Marc Carlson and H. Pagès",
        title = "GenomeInfoDb: Utilities for manipulating chromosome and other 'seqname' identifiers",
        year = 2017, doi = "10.18129/B9.bioc.GenomeInfoDb"
    ),
    GenomicRanges = citation("GenomicRanges"),
    ggplot2 = citation("ggplot2"),
    bumphunter = citation("bumphunter")[1],
    TxDb.Hsapiens.UCSC.hg19.knownGene = citation("TxDb.Hsapiens.UCSC.hg19.knownGene"),
    AnnotationDbi = RefManageR::BibEntry(
        bibtype = "manual",
        key = "AnnotationDbi",
        author = "Hervé Pagès and Marc Carlson and Seth Falcon and Nianhua Li",
        title = "AnnotationDbi: Annotation Database Interface",
        year = 2017, doi = "10.18129/B9.bioc.AnnotationDbi"
    ),
    BiocParallel = citation("BiocParallel"),
    derfinderHelper = citation("derfinderHelper"),
    GenomicAlignments = citation("GenomicAlignments"),
    GenomicFeatures = citation("GenomicFeatures"),
    GenomicFiles = citation("GenomicFiles"),
    Hmisc = citation("Hmisc"),
    qvalue = citation("qvalue"),
    Rsamtools = citation("Rsamtools"),
    rtracklayer = citation("rtracklayer"),
    S4Vectors = RefManageR::BibEntry(
        bibtype = "manual", key = "S4Vectors",
        author = "Hervé Pagès and Michael Lawrence and Patrick Aboyoun",
        title = "S4Vectors: S4 implementation of vector-like and list-like objects",
        year = 2017, doi = "10.18129/B9.bioc.S4Vectors"
    ),
    bumphunterPaper = RefManageR::BibEntry(
        bibtype = "article",
        key = "bumphunterPaper",
        title = "Bump hunting to identify differentially methylated regions in epigenetic epidemiology studies",
        author = "Jaffe, Andrew E and Murakami, Peter and Lee, Hwajin and Leek, Jeffrey T and Fallin, M Daniele and Feinberg, Andrew P and Irizarry, Rafael A",
        year = 2012, journal = "International Journal of Epidemiology"
    ),
    derfinderData = citation("derfinderData"),
    RefManageR = citation("RefManageR")[1]
)

## ----'ultraQuick', eval = FALSE-----------------------------------------------
#  ## Load libraries
#  library("derfinder")
#  library("derfinderData")
#  library("GenomicRanges")
#  
#  ## Determine the files to use and fix the names
#  files <- rawFiles(system.file("extdata", "AMY", package = "derfinderData"),
#      samplepatt = "bw", fileterm = NULL
#  )
#  names(files) <- gsub(".bw", "", names(files))
#  
#  ## Load the data from disk -- On Windows you have to load data from Bam files
#  fullCov <- fullCoverage(files = files, chrs = "chr21", verbose = FALSE)
#  
#  ## Get the region matrix of Expressed Regions (ERs)
#  regionMat <- regionMatrix(fullCov, cutoff = 30, L = 76, verbose = FALSE)
#  
#  ## Get pheno table
#  pheno <- subset(brainspanPheno, structure_acronym == "AMY")
#  
#  ## Identify which ERs are differentially expressed, that is, find the DERs
#  library("DESeq2")
#  
#  ## Round matrix
#  counts <- round(regionMat$chr21$coverageMatrix)
#  
#  ## Round matrix and specify design
#  dse <- DESeqDataSetFromMatrix(counts, pheno, ~ group + gender)
#  
#  ## Perform DE analysis
#  dse <- DESeq(dse, test = "LRT", reduced = ~gender, fitType = "local")
#  
#  ## Extract results
#  mcols(regionMat$chr21$regions) <- c(mcols(regionMat$chr21$regions), results(dse))
#  
#  ## Save info in an object with a shorter name
#  ers <- regionMat$chr21$regions
#  ers

## ----'start', message=FALSE---------------------------------------------------
## Load libraries
library("derfinder")
library("derfinderData")
library("GenomicRanges")

## ----'locateAMYfiles'---------------------------------------------------------
## Determine the files to use and fix the names
files <- rawFiles(system.file("extdata", "AMY", package = "derfinderData"),
    samplepatt = "bw", fileterm = NULL
)
names(files) <- gsub(".bw", "", names(files))

## ----'getData', eval = .Platform$OS.type != "windows"-------------------------
## Load the data from disk
fullCov <- fullCoverage(
    files = files, chrs = "chr21", verbose = FALSE,
    totalMapped = rep(1, length(files)), targetSize = 1
)

## ----'getDataWindows', eval = .Platform$OS.type == "windows", echo = FALSE----
#  ## Load data in Windows case
#  foo <- function() {
#      load(system.file("extdata", "fullCov", "fullCovAMY.RData",
#          package = "derfinderData"
#      ))
#      return(fullCovAMY)
#  }
#  fullCov <- foo()

## ----'regionMatrix'-----------------------------------------------------------
## Get the region matrix of Expressed Regions (ERs)
regionMat <- regionMatrix(fullCov, cutoff = 30, L = 76, verbose = FALSE)

## ----'exploreRegMatRegs'------------------------------------------------------
## regions output
regionMat$chr21$regions

## ----'exploreRegMatBP'--------------------------------------------------------
## Base-level coverage matrices for each of the regions
## Useful for plotting
lapply(regionMat$chr21$bpCoverage[1:2], head, n = 2)

## Check dimensions. First region is 565 long, second one is 138 bp long.
## The columns match the number of samples (12 in this case).
lapply(regionMat$chr21$bpCoverage[1:2], dim)

## ----'exploreRegMatrix'-------------------------------------------------------
## Dimensions of the coverage matrix
dim(regionMat$chr21$coverageMatrix)

## Coverage for each region. This matrix can then be used with limma or other pkgs
head(regionMat$chr21$coverageMatrix)

## ----'phenoData'--------------------------------------------------------------
## Get pheno table
pheno <- subset(brainspanPheno, structure_acronym == "AMY")

## ----'identifyDERsDESeq2'-----------------------------------------------------
## Identify which ERs are differentially expressed, that is, find the DERs
library("DESeq2")

## Round matrix
counts <- round(regionMat$chr21$coverageMatrix)

## Round matrix and specify design
dse <- DESeqDataSetFromMatrix(counts, pheno, ~ group + gender)

## Perform DE analysis
dse <- DESeq(dse, test = "LRT", reduced = ~gender, fitType = "local")

## Extract results
mcols(regionMat$chr21$regions) <- c(
    mcols(regionMat$chr21$regions),
    results(dse)
)

## Save info in an object with a shorter name
ers <- regionMat$chr21$regions
ers

## ----'vennRegions', fig.cap = "Venn diagram showing ERs by annotation class."----
## Find overlaps between regions and summarized genomic annotation
annoRegs <- annotateRegions(ers, genomicState$fullGenome, verbose = FALSE)

library("derfinderPlot")
venn <- vennRegions(annoRegs,
    counts.col = "blue",
    main = "Venn diagram using TxDb.Hsapiens.UCSC.hg19.knownGene annotation"
)

## ----'nearestGene'------------------------------------------------------------
## Load database of interest
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- keepSeqlevels(TxDb.Hsapiens.UCSC.hg19.knownGene, "chr21")

## Find nearest feature
library("bumphunter")
genes <- annotateTranscripts(txdb)
annotation <- matchGenes(ers, subject = genes)

## Restore seqlevels
txdb <- restoreSeqlevels(txdb)

## View annotation results
head(annotation)
## You can use derfinderPlot::plotOverview() to visualize this information

## ----'firstfive', fig.cap = paste0("Base-pair resolution plot of differentially expressed region ", 1:5, ".")----
## Extract the region coverage
regionCov <- regionMat$chr21$bpCoverage
plotRegionCoverage(
    regions = ers, regionCoverage = regionCov,
    groupInfo = pheno$group, nearestAnnotation = annotation,
    annotatedRegions = annoRegs, whichRegions = seq_len(5), txdb = txdb,
    scalefac = 1, ask = FALSE, verbose = FALSE
)

## ----createVignette, eval=FALSE-----------------------------------------------
#  ## Create the vignette
#  library("rmarkdown")
#  system.time(render("derfinder-quickstart.Rmd", "BiocStyle::html_document"))
#  
#  ## Extract the R code
#  library("knitr")
#  knit("derfinder-quickstart.Rmd", tangle = TRUE)

## ----reproducibility1, echo=FALSE---------------------------------------------
## Date the vignette was generated
Sys.time()

## ----reproducibility2, echo=FALSE---------------------------------------------
## Processing time in seconds
totalTime <- diff(c(startTime, Sys.time()))
round(totalTime, digits = 3)

## ----reproducibility3, echo=FALSE-------------------------------------------------------------------------------------
## Session info
library("sessioninfo")
options(width = 120)
session_info()

## ----vignetteBiblio, results = 'asis', echo = FALSE, warning = FALSE, message = FALSE---------------------------------
## Print bibliography
PrintBibliography(bib, .opts = list(hyperlink = "to.doc", style = "html"))