An Introduction to MMAPPR2

In MMAPPR2: Mutation Mapping Analysis Pipeline for Pooled RNA-Seq

```{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
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):

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 and make sure any differences are accounted for in the VEPFlags object.

Note: If you have any trouble installing VEP, using their Docker image 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:

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.

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.

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:

mmapprData <- mmappr(param)

The MMAPPR2 pipeline can also be run a step at a time:

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:

## 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:

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. You can configure the VEPFlags object like this:

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:

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:

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

sessionInfo()

Sys.setenv(PATH=origPath)

MMAPPR2 documentation built on Nov. 8, 2020, 6:53 p.m.