Nothing
## ----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]
)
## ----'start', message=FALSE-------------------------------------------------------------------------------------------
## Load libraries
library("derfinder")
library("derfinderData")
library("GenomicRanges")
## ----'phenoData', results = 'asis'------------------------------------------------------------------------------------
library("knitr")
## Get pheno table
pheno <- subset(brainspanPheno, structure_acronym == "AMY")
## Display the main information
p <- pheno[, -which(colnames(pheno) %in% c(
"structure_acronym",
"structure_name", "file"
))]
rownames(p) <- NULL
kable(p, format = "html", row.names = TRUE)
## ----'getData', eval = .Platform$OS.type != "windows"-----------------------------------------------------------------
## 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
system.time(fullCov <- fullCoverage(
files = files, chrs = "chr21",
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()
## ----'webData', eval = FALSE------------------------------------------------------------------------------------------
# ## Determine the files to use and fix the names
# files <- pheno$file
# names(files) <- gsub(".AMY", "", pheno$lab)
#
# ## Load the data from the web
# system.time(fullCov <- fullCoverage(
# files = files, chrs = "chr21",
# totalMapped = rep(1, length(files)), targetSize = 1
# ))
## ----'exploreFullCov'-------------------------------------------------------------------------------------------------
## Lets explore it
fullCov
## ----'filterCov'------------------------------------------------------------------------------------------------------
## Filter coverage
filteredCov <- lapply(fullCov, filterData, cutoff = 2)
## ----'exploreFilteredCov'---------------------------------------------------------------------------------------------
## Similar to fullCov but with $position
filteredCov
## ----'compareCov'-----------------------------------------------------------------------------------------------------
## Compare the size in Mb
round(c(
fullCov = object.size(fullCov),
filteredCov = object.size(filteredCov)
) / 1024^2, 1)
## ----'regionMatrix'---------------------------------------------------------------------------------------------------
## Use regionMatrix()
system.time(regionMat <- regionMatrix(fullCov, cutoff = 30, L = 76))
## Explore results
class(regionMat)
names(regionMat$chr21)
## ----'exploreRegMatRegs'----------------------------------------------------------------------------------------------
## regions output
regionMat$chr21$regions
## Number of regions
length(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)
## ----'identifyDERsDESeq2'---------------------------------------------------------------------------------------------
## Required
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
deseq <- regionMat$chr21$regions
mcols(deseq) <- c(mcols(deseq), results(dse))
## Explore the results
deseq
## ----'buildLimmaModels'-----------------------------------------------------------------------------------------------
## Build models
mod <- model.matrix(~ pheno$group + pheno$gender)
mod0 <- model.matrix(~ pheno$gender)
## ----'transformForLimma'----------------------------------------------------------------------------------------------
## Transform coverage
transformedCov <- log2(regionMat$chr21$coverageMatrix + 32)
## ----'identifyDERsLimma'----------------------------------------------------------------------------------------------
## Example using limma
library("limma")
## Run limma
fit <- lmFit(transformedCov, mod)
fit0 <- lmFit(transformedCov, mod0)
## Determine DE status for the regions
## Also in https://github.com/LieberInstitute/jaffelab with help and examples
getF <- function(fit, fit0, theData) {
rss1 <- rowSums((fitted(fit) - theData)^2)
df1 <- ncol(fit$coef)
rss0 <- rowSums((fitted(fit0) - theData)^2)
df0 <- ncol(fit0$coef)
fstat <- ((rss0 - rss1) / (df1 - df0)) / (rss1 / (ncol(theData) - df1))
f_pval <- pf(fstat, df1 - df0, ncol(theData) - df1, lower.tail = FALSE)
fout <- cbind(fstat, df1 - 1, ncol(theData) - df1, f_pval)
colnames(fout)[2:3] <- c("df1", "df0")
fout <- data.frame(fout)
return(fout)
}
ff <- getF(fit, fit0, transformedCov)
## Get the p-value and assign it to the regions
limma <- regionMat$chr21$regions
limma$fstat <- ff$fstat
limma$pvalue <- ff$f_pval
limma$padj <- p.adjust(ff$f_pval, "BH")
## Explore the results
limma
## ----limmaVSdeseq2----------------------------------------------------------------------------------------------------
table(limma$padj < 0.05, deseq$padj < 0.05)
## ----'createMeanBW', eval = .Platform$OS.type != "windows"------------------------------------------------------------
## Calculate the mean: this step takes a long time with many samples
meanCov <- Reduce("+", fullCov$chr21) / ncol(fullCov$chr21)
## Save it on a bigwig file called meanChr21.bw
createBw(list("chr21" = DataFrame("meanChr21" = meanCov)),
keepGR =
FALSE
)
## ----'railMatrix', eval = .Platform$OS.type != "windows"--------------------------------------------------------------
## Identify files to use
summaryFile <- "meanChr21.bw"
## We had already found the sample BigWig files and saved it in the object 'files'
## Lets just rename it to sampleFiles for clarity.
sampleFiles <- files
## Get the regions
system.time(
regionMat.rail <- railMatrix(
chrs = "chr21", summaryFiles = summaryFile,
sampleFiles = sampleFiles, L = 76, cutoff = 30, maxClusterGap = 3000L
)
)
## ----'checkDifferences', eval = .Platform$OS.type != "windows"--------------------------------------------------------
## Overall not identical due to small rounding errors
identical(regionMat, regionMat.rail)
## Actual regions are the same
identical(ranges(regionMat$chr21$regions), ranges(regionMat.rail$chr21$regions))
## When you round, the small differences go away
identical(
round(regionMat$chr21$regions$value, 4),
round(regionMat.rail$chr21$regions$value, 4)
)
identical(
round(regionMat$chr21$regions$area, 4),
round(regionMat.rail$chr21$regions$area, 4)
)
## ----'libSize'--------------------------------------------------------------------------------------------------------
## Get some idea of the library sizes
sampleDepths <- sampleDepth(collapseFullCoverage(fullCov), 1)
sampleDepths
## ----'makeModels'-----------------------------------------------------------------------------------------------------
## Define models
models <- makeModels(sampleDepths,
testvars = pheno$group,
adjustvars = pheno[, c("gender")]
)
## Explore the models
lapply(models, head)
## ----'analyze'--------------------------------------------------------------------------------------------------------
## Create a analysis directory
dir.create("analysisResults")
originalWd <- getwd()
setwd(file.path(originalWd, "analysisResults"))
## Perform differential expression analysis
system.time(results <- analyzeChr(
chr = "chr21", filteredCov$chr21, models,
groupInfo = pheno$group, writeOutput = TRUE, cutoffFstat = 5e-02,
nPermute = 20, seeds = 20140923 + seq_len(20), returnOutput = TRUE
))
## ----'exploreResults'-------------------------------------------------------------------------------------------------
## Explore
names(results)
## ----'exploreOptionsStats'--------------------------------------------------------------------------------------------
## Explore optionsStats
names(results$optionsStats)
## Call used
results$optionsStats$analyzeCall
## ----'exploreCovPrep'-------------------------------------------------------------------------------------------------
## Explore coveragePrep
names(results$coveragePrep)
## Group means
results$coveragePrep$groupMeans
## ----'exploreFstats'--------------------------------------------------------------------------------------------------
## Explore optionsStats
results$fstats
## Note that the length matches the number of bases used
identical(length(results$fstats), sum(results$coveragePrep$position))
## ----'exploreRegs'----------------------------------------------------------------------------------------------------
## Explore regions
names(results$regions)
## ----'exploreRegs2'---------------------------------------------------------------------------------------------------
## Permutation summary information
results$regions[2:4]
## ----'exploreRegs3'---------------------------------------------------------------------------------------------------
## Candidate DERs
results$regions$regions
## ----'sensitivity'----------------------------------------------------------------------------------------------------
## Width of potential DERs
summary(width(results$regions$regions))
sum(width(results$regions$regions) > 50)
## Width of candidate DERs
sig <- as.logical(results$regions$regions$significant)
summary(width(results$regions$regions[sig]))
sum(width(results$regions$regions[sig]) > 50)
## ----'exploreAnnotation'----------------------------------------------------------------------------------------------
## Nearest annotation
head(results$annotation)
## ----'exploreTime', fig.cap = "Seconds used to run each step in analyzeChr()."----------------------------------------
## Time spent
results$timeinfo
## Use this information to make a plot
timed <- diff(results$timeinfo)
timed.df <- data.frame(Seconds = as.numeric(timed), Step = factor(names(timed),
levels = rev(names(timed))
))
library("ggplot2")
ggplot(timed.df, aes(y = Step, x = Seconds)) +
geom_point()
## ----'mergeResults'---------------------------------------------------------------------------------------------------
## Go back to the original directory
setwd(originalWd)
## Merge results from several chromosomes. In this case we only have one.
mergeResults(
chrs = "chr21", prefix = "analysisResults",
genomicState = genomicState$fullGenome,
optionsStats = results$optionsStats
)
## Files created by mergeResults()
dir("analysisResults", pattern = ".Rdata")
## ----'optionsMerge'---------------------------------------------------------------------------------------------------
## Options used to merge
load(file.path("analysisResults", "optionsMerge.Rdata"))
## Contents
names(optionsMerge)
## Merge call
optionsMerge$mergeCall
## ----'exploreFullRegs'------------------------------------------------------------------------------------------------
## Load all the regions
load(file.path("analysisResults", "fullRegions.Rdata"))
## Metadata columns
names(mcols(fullRegions))
## ----'exploreFullAnnoRegs'--------------------------------------------------------------------------------------------
## Load annotateRegions() output
load(file.path("analysisResults", "fullAnnotatedRegions.Rdata"))
## Information stored
names(fullAnnotatedRegions)
## Take a peak
lapply(fullAnnotatedRegions, head)
## ----'extra'----------------------------------------------------------------------------------------------------------
## Find overlaps between regions and summarized genomic annotation
annoRegs <- annotateRegions(fullRegions, genomicState$fullGenome)
## Indeed, the result is the same because we only used chr21
identical(annoRegs, fullAnnotatedRegions)
## Get the region coverage
regionCov <- getRegionCoverage(fullCov, fullRegions)
## Explore the result
head(regionCov[[1]])
## ----'derfinderPlot', eval = FALSE------------------------------------------------------------------------------------
# library("derfinderPlot")
#
# ## Overview of the candidate DERs in the genome
# plotOverview(
# regions = fullRegions, annotation = results$annotation,
# type = "fwer"
# )
#
# suppressPackageStartupMessages(library("TxDb.Hsapiens.UCSC.hg19.knownGene"))
# txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
#
# ## Base-levle coverage plots for the first 10 regions
# plotRegionCoverage(
# regions = fullRegions, regionCoverage = regionCov,
# groupInfo = pheno$group, nearestAnnotation = results$annotation,
# annotatedRegions = annoRegs, whichRegions = 1:10, txdb = txdb, scalefac = 1,
# ask = FALSE
# )
#
# ## Cluster plot for the first region
# plotCluster(
# idx = 1, regions = fullRegions, annotation = results$annotation,
# coverageInfo = fullCov$chr21, txdb = txdb, groupInfo = pheno$group,
# titleUse = "fwer"
# )
## ----'featureLevel'---------------------------------------------------------------------------------------------------
## Get the exon-level matrix
system.time(exonCov <- coverageToExon(fullCov, genomicState$fullGenome, L = 76))
## Dimensions of the matrix
dim(exonCov)
## Explore a little bit
tail(exonCov)
## ----'regionMatAnnotate'----------------------------------------------------------------------------------------------
## Annotate regions as exonic, intronic or intergenic
system.time(annoGenome <- annotateRegions(
regionMat$chr21$regions,
genomicState$fullGenome
))
## Note that the genomicState object included in derfinder only has information
## for chr21 (hg19).
## Identify closest genes to regions
suppressPackageStartupMessages(library("bumphunter"))
suppressPackageStartupMessages(library("TxDb.Hsapiens.UCSC.hg19.knownGene"))
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
genes <- annotateTranscripts(txdb)
system.time(annoNear <- matchGenes(regionMat$chr21$regions, genes))
## ----'static-vis', eval = FALSE---------------------------------------------------------------------------------------
# ## Identify the top regions by highest total coverage
# top <- order(regionMat$chr21$regions$area, decreasing = TRUE)[1:100]
#
# ## Base-level plots for the top 100 regions with transcript information
# library("derfinderPlot")
# plotRegionCoverage(regionMat$chr21$regions,
# regionCoverage = regionMat$chr21$bpCoverage,
# groupInfo = pheno$group, nearestAnnotation = annoNear,
# annotatedRegions = annoGenome, whichRegions = top, scalefac = 1,
# txdb = txdb, ask = FALSE
# )
## ----'epivizr', eval = FALSE------------------------------------------------------------------------------------------
# ## Load epivizr, it's available from Bioconductor
# library("epivizr")
#
# ## Load data to your browser
# mgr <- startEpiviz()
# ders_dev <- mgr$addDevice(
# fullRegions[as.logical(fullRegions$significantFWER)], "Candidate DERs"
# )
# ders_potential_dev <- mgr$addDevice(
# fullRegions[!as.logical(fullRegions$significantFWER)], "Potential DERs"
# )
# regs_dev <- mgr$addDevice(regionMat$chr21$regions, "Region Matrix")
#
# ## Go to a place you like in the genome
# mgr$navigate(
# "chr21", start(regionMat$chr21$regions[top[1]]) - 100,
# end(regionMat$chr21$regions[top[1]]) + 100
# )
#
# ## Stop the navigation
# mgr$stopServer()
## ----'exportBigWig', eval = .Platform$OS.type != "windows"------------------------------------------------------------
## Subset only the first sample
fullCovSmall <- lapply(fullCov, "[", 1)
## Export to BigWig
bw <- createBw(fullCovSmall)
## See the file. Note that the sample name is used to name the file.
dir(pattern = ".bw")
## Internally createBw() coerces each sample to a GRanges object before
## exporting to a BigWig file. If more than one sample was exported, the
## GRangesList would have more elements.
bw
## ----'exampleNameStyle', eval = FALSE---------------------------------------------------------------------------------
# ## Set global species and chrsStyle options
# options(species = "arabidopsis_thaliana")
# options(chrsStyle = "NCBI")
#
# ## Then proceed to load and analyze the data
## ----'analyzeNonHuman', eval = FALSE----------------------------------------------------------------------------------
# ## Load transcript database information
# library("TxDb.Athaliana.BioMart.plantsmart28")
#
# ## Set organism options
# options(species = "arabidopsis_thaliana")
# options(chrsStyle = "NCBI")
#
# ## Run command with more arguments
# analyzeChr(txdb = TxDb.Athaliana.BioMart.plantsmart28)
## ----'runExample'-----------------------------------------------------------------------------------------------------
## Find some regions to work with
example("loadCoverage", "derfinder")
example("getRegionCoverage", "derfinder")
## ----'loadWhich'------------------------------------------------------------------------------------------------------
## Illustrate reading data from a set of regions
test <- loadCoverage(
files = files, chr = "21", cutoff = NULL, which = regions,
protectWhich = 0, fileStyle = "NCBI"
)
## Some reads were ignored and thus the coverage is lower as can be seen below:
sapply(test$coverage, max) - sapply(genomeDataRaw$coverage, max)
## ----'loadWhich2'-----------------------------------------------------------------------------------------------------
## Illustrate reading data from a set of regions
test2 <- loadCoverage(
files = files, chr = "21", cutoff = NULL,
which = regions, protectWhich = 3e4, fileStyle = "NCBI"
)
## Adding some padding to the regions helps get the same coverage
identical(sapply(test2$coverage, max), sapply(genomeDataRaw$coverage, max))
## A more detailed test reveals that the coverage matches at every base
all(mapply(function(x, y) {
identical(x, y)
}, test2$coverage, genomeDataRaw$coverage))
## ----Figure1, out.width="100%", fig.align="center", fig.cap = "derfinder F-statistics flow chart.", echo = FALSE------
knitr::include_graphics("http://lcolladotor.github.io/derfinder/fig/DERpathway.png")
## ----Figure2, out.width="100%", fig.align="center", fig.cap = "analyzeChr() flow chart.", echo = FALSE----------------
knitr::include_graphics("http://lcolladotor.github.io/derfinder/fig/analyzeChr.png")
## ----Figure3, out.width="100%", fig.align="center", fig.cap = "regionMatrix() flow chart.", echo = FALSE--------------
knitr::include_graphics("http://lcolladotor.github.io/derfinder/fig/regionMatrix.png")
## ----'derfinder-analysis', eval = FALSE-------------------------------------------------------------------------------
# ## Run derfinder's analysis steps with timing info
#
# ## Load libraries
# library("getopt")
#
# ## Available at http:/bioconductor.org/packages/derfinder
# library("derfinder")
#
# ## Specify parameters
# spec <- matrix(c(
# "DFfile", "d", 1, "character", "path to the .Rdata file with the results from loadCoverage()",
# "chr", "c", 1, "character", "Chromosome under analysis. Use X instead of chrX.",
# "mcores", "m", 1, "integer", "Number of cores",
# "verbose", "v", 2, "logical", "Print status updates",
# "help", "h", 0, "logical", "Display help"
# ), byrow = TRUE, ncol = 5)
# opt <- getopt(spec)
#
# ## Testing the script
# test <- FALSE
# if (test) {
# ## Speficy it using an interactive R session and testing
# test <- TRUE
# }
#
# ## Test values
# if (test) {
# opt <- NULL
# opt$DFfile <- "/ProjectDir/derCoverageInfo/chr21Cov.Rdata"
# opt$chr <- "21"
# opt$mcores <- 1
# opt$verbose <- NULL
# }
#
# ## if help was asked for print a friendly message
# ## and exit with a non-zero error code
# if (!is.null(opt$help)) {
# cat(getopt(spec, usage = TRUE))
# q(status = 1)
# }
#
# ## Default value for verbose = TRUE
# if (is.null(opt$verbose)) opt$verbose <- TRUE
#
# if (opt$verbose) message("Loading Rdata file with the output from loadCoverage()")
# load(opt$DFfile)
#
# ## Make it easy to use the name later. Here I'm assuming the names were generated using output='auto' in loadCoverage()
# eval(parse(text = paste0("data <- ", "chr", opt$chr, "CovInfo")))
# eval(parse(text = paste0("rm(chr", opt$chr, "CovInfo)")))
#
# ## Just for testing purposes
# if (test) {
# tmp <- data
# tmp$coverage <- tmp$coverage[1:1e6, ]
# library("IRanges")
# tmp$position[which(tmp$pos)[1e6 + 1]:length(tmp$pos)] <- FALSE
# data <- tmp
# }
#
# ## Load the models
# load("models.Rdata")
#
# ## Load group information
# load("groupInfo.Rdata")
#
#
# ## Run the analysis with lowMemDir
# analyzeChr(
# chr = opt$chr, coverageInfo = data, models = models,
# cutoffFstat = 1e-06, cutoffType = "theoretical", nPermute = 1000,
# seeds = seq_len(1000), maxClusterGap = 3000, groupInfo = groupInfo,
# subject = "hg19", mc.cores = opt$mcores,
# lowMemDir = file.path(tempdir(), paste0("chr", opt$chr), "chunksDir"),
# verbose = opt$verbose, chunksize = 1e5
# )
#
# ## Done
# if (opt$verbose) {
# print(proc.time())
# print(sessionInfo(), locale = FALSE)
# }
## ----createVignette, eval=FALSE---------------------------------------------------------------------------------------
# ## Create the vignette
# library("rmarkdown")
# system.time(render("derfinder-users-guide.Rmd", "BiocStyle::html_document"))
#
# ## Extract the R code
# library("knitr")
# knit("derfinder-users-guide.Rmd", tangle = TRUE)
## ----createVignette2--------------------------------------------------------------------------------------------------
## Clean up
file.remove("derfinderUsersGuideRef.bib")
unlink("analysisResults", recursive = TRUE)
file.remove("HSB113.bw")
file.remove("meanChr21.bw")
## ----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"))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.