The polymerase kinetics data generated by PacBio's SMRT sequencing technology allows simultaneous detection of common microbial base modifications in native genomic DNA in resequencing experiments. SMRTanalysis provides two tools for interpreting the results: kineticsTools to identify modified bases, and MotifMaker to detect consensus motifs associated with modified bases. These are available in SMRTlink as an integrated pipeline.

• Overview
• Manual
  • Running with SMRTLink
  • Running on the Command Line
  • Running on the Command Line with pbsmrtpipe
• Advanced Analysis Options
  • SMRTLink/pbsmrtpipe BaseMod Options
  • PBAlign Options
  • kineticsTools Options
  • MotifMaker Options
• Files
  • PBAlign Files
  • kineticsTools Files
  • MotifMaker Files
• Algorithm Modules
• Glossary

Analyses are performed in four stages: PBAlign, kineticsTools, MotifMaker Find, and MotifMaker Reprocess.

• PBAlign
  • PBAlign maps PacBio sequences to references using an algorithm selected from a selection of supported command-line alignment algorithms. Input can be a fasta, pls.h5, bas.h5 or ccs.h5 file or a fofn (file of file names). Output can be in CMP.H5, SAM or BAM format. If output is BAM format, aligner can only be BLASR and QVs will be loaded automatically.
• kineticsTools
  • kineticsTools takes the AlignmentSet (mapped BAM) output of pbalign and iterates over all nucleotides in the genome. From the reads aligned at each site, a profile of polymerase kinetics is obtained in the form of inter-pulse durations (IPDs), which are compared to expectations for the given chemistry. The ratio of observed-to-expected IPD is interpreted to decide whether a base is modified and (optionally) determine the specific modification type.
• MotifMaker
  • MotifMaker is a tool for identify motifs associated with DNA modifications in prokaryotic genomes. Modified DNA in prokaryotes commonly arises from restriction-modification systems that methylate a specific base in a specific sequence motif. The canonical example the m6A methylation of adenine in GATC contexts in E.coli. Prokaryotes may have a very large number of active restriction-modification systems present, leading to a complicated mixture of sequence motifs. PacBio SMRT sequencing is sensitive to the presence of methylated DNA at single base resolution, via shifts in the polymerase kinetics observed in the real-time sequencing traces. See our publication for more background on modification detection.

To run Isoseq using SMRTLink, follow the usual steps for analysing data on SMRTLink. TODO: Link to document explaining SMRTLink.

On the command line, the analysis is performed in 4 steps:

1. Run PBAlign on your unaligned subreads, generating an aligned BAM.
2. Run kineticsTools on your aligned BAM, generating a GFF and CSV containing base modification information.
3. Run MotifMaker on your GFF generated from kineticsTools, generating a CSV containing motif information.
4. Run MotifMaker on your base modification GFF and motif CSV, generating a motif GFF.

Step 1. PBAlign

First, align your sequences to your chosen reference.

 pbalign --concordant --hitPolicy=randombest --minAccuracy 70 --minLength 50 --algorithmOptions="-minMatch 12 -bestn 10 -minPctIdentity 70.0" subreads.bam reference.fasta aligned_subreads.bam

Where reference.fasta contains your reference sequences.

Where subreads.bam contains your unaligned reads.

Where aligned_subreads.bam is where your aligned reads will be stored.

Step 2. kineticsTools

Next, analyze your aligned sequences for base modifications.

 ipdSummary aligned_subreads.bam --reference reference.fasta --gff basemods.gff --csv basemods.csv --pvalue 0.001 --numWorkers 16 --identify m4C,m6A

Note: You will need to have an index file for your reference.fasta and for your aligned_subreads.bam.

To index your reference.fasta, use samtools faidx reference.fasta.

To index your aligned_subreads.bam, use pbindex aligned_subreads.bam.

Step 3. MotifMaker Find

Next, identify consensus motifs.

 motifMaker find -f reference.fasta -g basemods.gff -o motifs.csv

Step 4. MotifMaker Reprocess

Finally, generate a GFF file of all modifications that are part of motifs.

 motifMaker reprocess -f reference.fasta -g basemods.gff -m motifs.csv -o motifs.gff

####Install pbsmrtpipe pbsmrtpipe is a part of smrtanalysis-3.0 package and will be installed if smrtanalysis-3.0 has been installed on your system. Or you can download pbsmrtpipe and install.

You can verify that pbsmrtpipe is running OK by:

pbsmrtpipe --help

Now create an XML file from your subreads.

dataset create --type SubreadSet my.subreadset.xml subreads1.bam subreads2.bam ...

This will create a file called my.subreadset.xml.

Create a global options XML file which contains SGE related, job chunking and job distribution options that you may modify by:

 pbsmrtpipe show-workflow-options -o global_options.xml

Create a basemod options XML file which contains basemod-related options that you may modify by:

 pbsmrtpipe show-template-details pbsmrtpipe.pipelines.ds_modification_motif_analysis -o basemod_options.xml

The entries in the options XML files have the format:

 <option id="pbtranscript.task_options.min_seq_len">
            <value>300</value>
        </option>

And you can modify options using your favorite text editor, such as vim.

Once you have set your options, you are ready to run basemod via pbsmrtpipe:

pbsmrtpipe pipeline-id pbsmrtpipe.pipelines.ds_modification_motif_analysis -e eid_ref_dataset:reference.fasta -e eid_subread:my.subreadset.xml --preset-xml=basemod_options.xml --preset-xml=global_options.xml

MotifMaker options, except for --help, require the specification of a command first- either find or reprocess.

Unaligned Reads File (subreads.*)

(Required) This file contains your input reads. It can be a BAM, FASTA, XML, CCS.H5 or BAS.H5 file.

Reference File(reference.fasta)

(Required) This file contains the reference sequences that your input reads will be aligned against.

Configuration File

(Optional) If you prefer to specify all arguments in a config file rather than on the command line, use the --configFile argument.

Example config text file:

# This is a config file for pbalign where users can specify
# values for an arbitary subset of optional arguments for pbalign.
# Lines which start with '#' are comments. Otherwise each line
# should specify the value for exactly one optinal argument.
# Note that:
#   [1] Positional arguments, including inputFileName, outputFileName,
#       and referencePath, can not be specified in a config file.
#   [2] Sepcial arguments, including --verbose, --version (-v..),
#       --debug and --profile, in a config file will be ignored.
#   [3] Arguments specified in a config file will be overwritten
#       by arguments on command-line.

# Aligner's options.
--maxHits       = 20
--minAnchorSize = 1
--minLength     = 100
--algorithmOptions = "-noSplitSubreads -maxMatch 30 -nCandidates 30"

# SamFilter's filtering criteria and hit policy.
--hitPolicy     = randombest
--maxDivergence = 40

# Miscellaneous
--seed = 10

Pulse File

(Optional) Path to *.bas.h5 or *.pls.h5 when input type is *.fasta, and output type is *.cmp.h5. This is necessary to load pulse metrics into the alignment file (*.cmp.h5) for subsequent consumption by variantCaller.

Region Table

(Optional) path to *.rgn.h5. When input file is of type *.bas.h5, you may supply an optional Region table (*.rgn.h5) to filter the reads prior to alignment

####Outputs

Aligned Reads File (aligned_subreads.*)

This file contains your aligned output reads. You can specify *.cmp.h5, *.sam or *.bam Note: *.bam output format only works if --algorthim blasr

Modifications CSV File (basemods.csv)

The modifications.csv file contains one row for each (reference position, strand) pair that appeared in the dataset with coverage at least x. x defaults to 3, but is configurable with ‘–minCoverage’ flag to ipdSummary.py. The reference position index is 1-based for compatibility with the gff file the R environment.

The output columns vary depending on whether --control was used to run in case-control mode, or whether the program ran in the default in-silico control mode.

in-silico control mode output columns

case-control mode output columns

Modifications GFF File (basemods.gff)

The modifications.gff is compliant with the GFF Version 3 specification. Each template position / strand pair whose p-value exceeds the pvalue threshold appears as a row. The template position is 1-based, per the GFF spec. The strand column refers to the strand carrying the detected modification, which is the opposite strand from those used to detect the modification. The GFF confidence column is a Phred-transformed pvalue of detection.

The modifications.gff file will not work directly with most genome browsers. You will likely need to make a copy of the GFF file and convert the seqid columns from the generic ‘ref0000x’ names generated by PacBio, to the FASTA headers present in the original reference FASTA file. The mapping table is written in the header of the modifications.gff file in #sequence-header tags. This issue will be resolved in the 1.4 release of kineticsTools.

The auxiliary data column of the GFF file contains other statistics which may be useful downstream analysis or filtering. In particular the coverage level of the reads used to make the call, and +/- 20bp sequence context surrounding the site.

Motifs GFF File (motifs.gff)

This is a reprocessed version of modifications.gff with the following changes. If a detected modification occurs on a motif detected by the motif finder, the modification is annotated with motif data. An attribute ‘motif’ is added containing the motif string, and an attribute ‘id’ is added containing the motif id, which is the motif string for unpaired motifs or ‘motifString1/motifString2’ for paired motifs. If a motif instance exists in the genome, but was not detected in modifications.gff, an entry is added to motifs.gff, indicating the presence of that motif and the kinetics that were observed at that site.

Motifs CSV File (motifs.csv)

This file summarizes the modified motifs discovered by the tool. The CSV contains one row per detected motif, with the following columns.

BLASR is a read mapping program that maps reads to positions in a genome by clustering short exact matches between the read and the genome, and scoring clusters using alignment. The matches are generated by searching all suffixes of a read against the genome using a suffix array. Global chaining methods are used to score clusters of matches. It is exremely useful to have read filtering information, and mapping runtime may decrease substantially when a precomputed suffix array index on the reference sequence is specified.

Studies of the relationship between IPD and sequence context reveal that most of the variation in mean IPD across a genome can be predicted from a 12-base sequence context surrounding the active site of the DNA polymerase. The bounds of the relevant context window correspond to the window of DNA in contact with the polymerase, as seen in DNA/polymerase crystal structures. To simplify the process of finding DNA modifications with PacBio data, the tool includes a pre-trained lookup table mapping 12-mer DNA sequences to mean IPDs observed in C2 chemistry.

kineticsTools uses the Mapping QV generated by BLASR and stored in the BAM file to ignore reads that aren’t confidently mapped. The default minimum Mapping QV required is 10, implying that BLASR has 90% confidence that the read is correctly mapped. Because of the range of read lengths inherent in PacBio data. This can be changed in using the –mapQvThreshold command line argument, or via the SMRT Link configuration dialog for Modification Detection.

There are a few features of PacBio data that require special attention in order to acheive good modification detection performance. kineticsTools inspects the alignment between the observed bases and the reference sequence – in order for an IPD measurement to be included in the analysis, the PacBio read sequence must match the reference sequence for k around the cognate base. In the current module k=1, the IPD distribution at some locus be thought of as a mixture between the ‘normal’ incorporation process IPD, which is sensitive to the local sequence context and DNA modifications and a contaminating ‘pause’ process IPD which have a much longer duration (mean >10x longer than normal), but happen rarely (~1% of IPDs). Note: Our current understanding is that pauses do not carry useful information about the methylation state of the DNA, however a more careful analysis may be warranted. Also note that modifications that drastically increase the Roughly 1% of observed IPDs are generated by pause events. Capping observed IPDs at the global 99th percentile is motivated by theory from robust hypothesis testing. Some sequence contexts may have naturally longer IPDs, to avoid capping too much data at those contexts, the cap threshold is adjusted per context as follows: capThreshold = max(global99, 5*modelPrediction, percentile(ipdObservations, 75))

We test the hypothesis that IPDs observed at a particular locus in the sample have a longer means than IPDs observed at the same locus in unmodified DNA. If we have generated a Whole Genome Amplified dataset, which removes DNA modifications, we use a case-control, two-sample t-test. This tool also provides a pre-calibrated ‘synthetic control’ model which predicts the unmodified IPD, given a 12 base sequence context. In the synthetic control case we use a one-sample t-test, with an adjustment to account for error in the synthetic control model.

Existing motif finding algorithms such as MEME-chip and YMF are sub-optimal for this case for the following reasons:

1. They search for a single motif, rather than attempting to identify a complicated mixture of motifs

2. They generally don't accept the notion of aligned motifs - the input to the tools is a window into the reference sequence which can contain the motif at any offset, rather a single center position that is available with kinetic modification detection.

3. Implementations generally either use a Markov model of the reference (MEME-chip), or do exact counting on the reference, but place restrictions on the size and complexity of the motifs that can be discovered.

Here we give a rough overview of the algorthim used by MotifMaker.

Define a motif as a set of (position relative to methylation, required base) tuples. Positions not listed in the motif are implicitly degenerate. Given a list of modification detections and a genome sequence, we define the following objective function on motifs:

Motif score(motif) = 
(# of detections matching motif) / (# of genome sites matching motif)  *  (Sum of log-pvalue of detections matching motif) = 
(fraction methylated) * (sum of log-pvalues of matches)

We search (close to exhaustively) through the space of all possible motifs, progressively testing longer motifs using a branch-and-bound search. The 'fraction methylated' term must be less than 1, so the maximum achievable score of a child node is the sum of scores of modification hits in the current node, permitting us to prune all search paths whose maximum achievable score is less than the best score discovered so far.

• Inter-Pulse Duration (IPD)
• In the terminology of SMRT sequencing, pulses are single-molecule fluorescence signals interpreted as basecalls. The IPD for a basecall is the time between the current pulse and the next one. Many covalent base modifications cause a temporary pause in the sequencing reaction, which varies depending on the type of modification. For detecting basemods, kineticsTools uses the IPD averaged over many subreads aligned at the same position.

^{For Research Use Only. Not for use in diagnostic procedures. © Copyright 2015, Pacific Biosciences of California, Inc. All rights reserved. Information in this document is subject to change without notice. Pacific Biosciences assumes no responsibility for any errors or omissions in this document. Certain notices, terms, conditions and/or use restrictions may pertain to your use of Pacific Biosciences products and/or third party products. Please refer to the applicable Pacific Biosciences Terms and Conditions of Sale and the applicable license terms at http://www.pacificbiosciences.com/licenses.html.}

^{Visit the PacBio Developer's Network Website for the most up-to-date links to downloads, documentation and more.}

^{Terms of Use & Trademarks | Contact Us}