Dependencies¶

Modification Detection¶

The basic mode of kineticsTools does an independent comparison of IPDs at each position on the genome, for each strand, and emits various statistics to CSV and GFF (after applying a significance filter).

Modifications Identification¶

Synthetic Control¶

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.

Filtering and Trimming¶

kineticsTools uses the Mapping QV generated by BLASR or pbmm2 and stored in the alignments or BAM file (or AlignmentSet) 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 readlengths inherent in PacBio dataThis can be changed in using the --mapQvThreshold command line argument, or via the SMRTPortal configuration dialog for Modification Detection.

There are a few features of PacBio data that require special attention in order to achieve 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))

Statistical Testing¶

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.

Aligned Reads¶

A standard PacBio alignment file - either AlignmentSet XML or BAM - containing alignments and IPD information supplies the kinetic data used to perform modification detection.

Reference Sequence¶

The tool requires the reference sequence used to perform alignments. This can be either a FASTA file or a ReferenceSet XML.

modifications.gff¶

The modifications.gff is compliant with the GFF Version 3 specification ( <http://www.sequenceontology.org/gff3.shtml> ). 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.

Note on genome browser compatibility

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.

modifications.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. Note that this output type scales poorly and is not recommended for large genomes; the HDF5 output should perform much better in these cases. We have preserved the CSV option to support legacy applications but this is no longer produce by the pipelines in SMRT Link/pbsmrtpipe.

Output columns¶

in-silico control mode

case-control mode