--- title: "An Introduction to MMAPPR2" author: - name: Kyle Johnsen affiliation: Brigham Young University, Provo, UT package: MMAPPR2 output: BiocStyle::html_document: toc_float: true BiocStyle::pdf_document: default abstract: | Instructions on mapping mutations from forward genetic screens using MMAPPR2. vignette: | %\VignetteIndexEntry{An Introduction to MMAPPR2} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{bash, echo = FALSE} # Cooking show code to account for cases when samtools is not installed echo '#!/bin/bash' > /tmp/samtools chmod 755 /tmp/samtools ``` ```{r, echo = FALSE} dataDir <- system.file('extdata', package = 'MMAPPR2data') WTpileupFile <- file.path(dataDir, 'exwt.plp') MTpileupFile <- file.path(dataDir, 'exmut.plp') samtoolsScript <- file('/tmp/samtools', "a") writeLines(c( 'if [[ ${@:$#} == *"wt.bam"* ]];', 'then', paste('cat', WTpileupFile), 'else', paste('cat', MTpileupFile), 'fi' ), samtoolsScript) close(samtoolsScript) origPath <- Sys.getenv('PATH') Sys.setenv(PATH = paste(origPath, '/tmp', sep = ':')) ``` # Getting Started ## Data You'll need BAM file(s) for your wild-type pool, BAM file(s) for your mutant pool, and the reference genome for your species in fasta format. We recommend that each pool contain at least 20 individuals to ensure a good number of recombinations to measure. ## Installing Dependencies MMAPPR2 depends on two system tools to function: Samtools and VEP. Both must be installed and in the PATH to be found by the appropriate functions. ### Installing Samtools Instructions to install samtools can be found at https://github.com/samtools/samtools and installation instructions are in the INSTALL file included with samtools. ### Installing VEP You'll need Ensembl VEP, which you can install like this, replacing `my_species` with your species (e.g., `danio_rerio`): ```{r installVEP, eval=FALSE} git clone https://github.com/Ensembl/ensembl-vep.git cd ensembl-vep perl INSTALL.pl -a ac -s {my_species} ``` This installs the most recent VEP and allows you to create a cache for your desired species, which is what MMAPPR2 expects by default. If you depart from the installation shown here, or if things don't go smoothly, see [Ensembl's instructions](http://www.ensembl.org/info/docs/tools/vep/script/vep_download.html#installer) and make sure any differences are accounted for in the [`VEPFlags`](#configure-vepflags-object) object. *Note:* If you have any trouble installing VEP, using [their Docker image](http://www.ensembl.org/info/docs/tools/vep/script/vep_download.html#docker) may save you a lot of hassle. *Note:* We have found that R sometimes has issues finding VEP, especially when perlbrew is used. If you encounter errors at the path to your perl installation to the .Rprofile file. For example: ```{r, eval = FALSE} Sys.setenv(PATH=paste("/Path/to/Perlbrew", Sys.getenv("PATH"), sep=":")) ``` # Basic Use ## Setting Parameters For our example, we will use just the golden gene from the GRCz11 zebrafish reference genome. Here we also configure the VEPFlags object to use our example fasta and GTF files. [See below for more info.](#configure-vepflags-object) **Make sure your reference genome is the same you'll use with VEP! This will be the most recent assembly available on Ensembl unless you customize. You should use the same genome in aligning your sequencing data as well.** ```{r param} BiocParallel::register(BiocParallel::MulticoreParam()) ## see below for explanation of BiocParallel library(MMAPPR2, quietly = TRUE) library(MMAPPR2data, quietly = TRUE) library(Rsamtools, quietly = TRUE) # This is normally configured automatically: vepFlags <- ensemblVEP::VEPFlags(flags = list( format = 'vcf', # <-- this is necessary vcf = FALSE, # <-- as well as this species = 'danio_rerio', database = FALSE, # <-- these three arguments allow us to run VEP offline, fasta = goldenFasta(), # <-╯| which you probably won't need gff = goldenGFF(), # <------╯ filter_common = TRUE, coding_only = TRUE # assuming RNA-seq data )) param <- MmapprParam(refFasta = goldenFasta(), wtFiles = exampleWTbam(), mutFiles = exampleMutBam(), species = 'danio_rerio', vepFlags = vepFlags, ## optional outputFolder = tempOutputFolder()) ## optional ``` ## Running MMAPPR2 With parameters set, running the pipeline should be as simple as the following: ```{r mmappr} mmapprData <- mmappr(param) ``` The MMAPPR2 pipeline can also be run a step at a time: ```{r mmappr-steps} md <- new('MmapprData', param = param) ## calculateDistance() takes a MmapprData object postCalcDistMD <- calculateDistance(md) postLoessMD <- loessFit(postCalcDistMD) postPrePeakMD <- prePeak(postLoessMD) postPeakRefMD <- peakRefinement(postPrePeakMD) postCandidatesMD <- generateCandidates(postPeakRefMD) outputMmapprData(postCandidatesMD) ``` If the pipeline fails midway, the `MmapprData` object is saved, which you can then load and use for inspection and debugging: ```{r recover-md} ## Contents of output folder: cat(paste(system2('ls', outputFolder(param(mmapprData)), stdout = TRUE)), sep = '\n') mdFile <- file.path(outputFolder(param(mmapprData)), 'mmappr_data.RDS') md <- readRDS(mdFile) md ``` ## Results If everything goes well you should be able to track down your mutation using the `candidates` slot of your `MmapprData` object or by looking at files in the output folder: ```{r results} head(candidates(mmapprData)$`18`, n=2) outputTsv <- file.path(outputFolder(param(mmapprData)), '18.tsv') cat(paste(system2('head', outputTsv, stdout = TRUE)), sep = '\n') ``` # Advanced Configuration ## Configure VEPFlags Object MMAPPR2 uses the `r BiocStyle::Biocpkg("ensemblVEP")` Bioconductor package to predict the effect of variants in the peak region. To customize this process, you'll need to configure a `VEPFlags` object. Look at [Ensembl's website for script options](http://uswest.ensembl.org/info/docs/tools/vep/script/vep_options.html). You can configure the `VEPFlags` object like this: ```{r vepFlags} library(ensemblVEP, quietly = TRUE) vepFlags <- VEPFlags(flags = list( ### DEFAULT SETTINGS format = 'vcf', # <-- this is necessary vcf = FALSE, # <-- as well as this species = 'danio_rerio', database = FALSE, cache = TRUE, filter_common = TRUE, coding_only = TRUE # assuming RNA-seq data ### YOU MAY FIND THESE INTERESTING: # everything = TRUE # enables many optional analyses, such as Polyphen and SIFT # per_gene = TRUE # will output only the most severe variant per gene # pick = TRUE # will output only one consequence per variant )) ``` ## BiocParallel Configuration MMAPPR2 simply uses the default `bpparam` registered. You can change this (for example, if `r BiocStyle::Biocpkg("BiocParallel")` isn't working correctly) with the `BiocParallel::register` command. For example: ```{r bpparam} library(BiocParallel, quietly = TRUE) register(SerialParam()) register(MulticoreParam(progressbar=TRUE)) registered() ``` The last param registered becomes the default. ## Reference Genome The variant calling step requires a `BiocStyle::Biocpkg("gmapR")` `GmapGenome`, which is normally automatically generated from the `refFasta` parameter. If for some reason you want to generate your own, the process is like this: ```{r refGenome, eval=FALSE} refGenome <- gmapR::GmapGenome(goldenFasta(), name='slc24a5', create=TRUE) ``` ## Whole-genome Sequencing (WGS) MMAPPR2, like its predecessor, was designed for and tested using RNA-Seq data. However, the principles at work should still apply for WGS data. ----------------- # Session Info ```{r sessionInfo} sessionInfo() ``` ```{r, echo = FALSE} Sys.setenv(PATH=origPath) ```