IPDSUMMARY(1)

📦
version 0.6.1+git20220223.1326a4d · 3 versions
, in debian bookworm, ubuntu jammy

✨ IPDSUMMARY — the ipdSummary(1) glow-up ✨

okay bestie, sit down, hydrate, because we're about to read a whole man page together and honestly? it's a bop. this is ipdSummary, section 1 of your manual, and it's all about sniffing out DNA base-modifications from kinetic signatures. 🧬💅

🌸 NAME 🌸

ipdSummary — Detect DNA base-modifications from kinetic signatures. that's it, that's the whole job, and it does it gorgeously.

🎀 DESCRIPTION 🎀

so kineticsTool (yes, the tool is called kineticsTool, keep up 💖) loads up all the IPDs observed at every position in the genome, and then it's like "hmm, is this giving unmodified, or is this giving modified?" It compares those observed IPDs against the value you'd expect for unmodified DNA, and then it outputs the result of that statistical test. slay.

and where does the expected IPD value for unmodified DNA even come from? two places, bestie:

  • an in-silico control, or
  • an amplified control.

The in silico control is trained by PacBio and shipped right there with the package. It predicts the IPD using the local sequence context around the current position. no sequencing required, it just knows. 👑

An amplified control dataset is generated by sequencing unmodified DNA with the same sequence as the test sample. And an amplified control sample is usually generated by whole-genome amplification of the original sample. very much "we made a clean copy so we'd have a baseline" energy.

🦋 Dependencies 🦋

kineticsTools depends on:

that's the whole friend group. no cap.

🌈 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). independent comparison = every position gets its own little vibe check. ✨

💜 Modifications Identification 💜

kineticsTools also has a Modification Identification mode that can decode multi-site IPD 'fingerprints' into a reduced set of calls of specific modifications. and this feature? it's giving main character. benefits:

  • Different modifications occurring on the same base can be distinguished (for example m5C and m4C)
  • The signal from one modification is combined into one statistic, improving sensitivity, removing extra peaks, and correctly centering the call

one statistic, no messy duplicates, perfectly centered. periodt. 💅

✨ OPTIONS ✨

Please call this program with --help to see the available options. that's literally the whole section, go ask it yourself, it'll tell you everything. 🫶

🔥 ALGORITHM 🔥

🧚 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. it's a cheat sheet and it's iconic. 📖⭐

🌸 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 data, this 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 can 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))
okay so the pauses? they're loud and dramatic and they mean nothing about methylation. cap them. respectfully. 💀 but cap them per context so you don't accidentally mute a context that's just naturally chatty.

⭐ 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. two-sample when you've got a control sample, one-sample + an adjustment when you're rolling with the synthetic model. it's giving rigorous. 📊

🌈 EXAMPLE USAGE 🌈

Basic use with BAM input, GFF output:

ipdSummary aligned.bam --reference ref.fasta --identify m6A,m4C --gff basemods.gff

With dataset input, methyl fraction calculation and GFF+CSV output:

ipdSummary mapped.alignmentset.xml --reference ref.fasta --identify m6A,m4C --methylFraction --gff basemods.gff --csv kinetics.csv

copy, paste, go off. 🫶

🎀 INPUTS 🎀

🦋 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.

💜 OUTPUTS 💜

The modification detection tool provides results in a variety of formats suitable for in-depth statistical analysis, quick reference, and consumption by visualization tools such as PacBio SMRTView. Results are generally indexed by reference position and reference strand. In all cases the strand value refers to the strand carrying the modification in DNA sample. Remember that the kinetic effect of the modification is observed in read sequences aligning to the opposite strand. So reads aligning to the positive strand carry information about modification on the negative strand and vice versa, but in this toolkit we alway report the strand containing the putative modification. 🧬

this strand thing is the #1 place people get confused, so say it with me: the reads that see the modification are on the opposite strand from the one we report. the toolkit always reports the strand with the putative modification. iykyk. 🙈

The following output options are available:

  • --gff FILENAME: GFF format
  • --csv FILENAME: comma-separated value format
  • --bigwig FILENAME: BigWig file (mostly only useful for SMRTView)

If you are running base modification analysis through SMRT Link or a pbsmrtpipe pipeline, the GFF, HDF5, and BigWig outputs are automatically generated. (yes, HDF5, even though it's not in that little bullet list — SMRT Link just does it for you 💅)

🌸 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.

lowkey a jump scare the first time you drag it into a browser and everything is named ref0000x. the mapping table is right there in the header though, bestie. you got this. 🍓

That auxiliary data column? It's not just decoration bestie, it's carrying extra stats you can totally use for downstream analysis or filtering. Specifically: the coverage level of the reads that made the call, plus the +/- 20bp sequence context chilling around the site. Obsessed.

Here's the whole GFF column breakdown, one by one, no skips:

  • seqid — the Fasta contig name 🧬
  • source — name of the tool that made it, which is 'kinModCall'
  • type — the modification type. In identification mode this is m6A, m4C, or m5C for identified bases, OR the generic tag 'modified_base' if a kinetic event got detected that does not match a known modification signature
  • start — modification position on the contig
  • end — modification position on the contig
  • score — the Phred transformed p-value of detection. This is the single-site detection p-value, periodt.
  • strand — the sample strand containing the modification
  • phase — not applicable, it's giving nothing, moving on
  • attributes — the extra fields relevant to base mods. Let me break this down because there's a lot in here:
    • IPDRatio — the traditional IPDRatio
    • context — the reference sequence, -20bp to +20bp around the modification
    • coverage level — the number of IPD observations used after Mapping QV filtering and accuracy filtering
    • if the row results from an identified modification, there's also an identificationQv tag with the from the modification identification procedure. identificationQv is the phred-transformed probability of an incorrect identification, for bases that were identified as having a particular modification.
    • frac, fracLow, fracUp — the estimated fraction of molecules carrying the modification, and the 5% confidence intervals of the estimate.
Heads up bestie: the methylated fraction estimation is a beta-level feature, and should only be used for exploratory purposes. It's a vibe check, not a verdict. 💅

🎀 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. Now, x defaults to 3, but you can configure it with the '--minCoverage' flag to ipdSummary.py. Slay, customizable.

The reference position index is 1-based for compatibility with the gff file the R environment. Keep that in mind so you don't off-by-one yourself into a whole situation. 😭

Real talk: this output type scales poorly and is not recommended for large genomes. The HDF5 output should perform much better in these cases. We've preserved the CSV option to support legacy applications, but this is no longer produce by the pipelines in SMRT Link/pbsmrtpipe. So like… it's here, it's cute, but it's legacy behavior. 🫶

✨ Output columns in-silico control mode ✨

  • refId — reference sequence ID of this observation
  • tpl — 1-based template position
  • strand — native sample strand where kinetics were generated. '0' is the strand of the original FASTA, '1' is opposite strand from FASTA
  • base — the cognate base at this position in the reference
  • score — Phred-transformed pvalue that a kinetic deviation exists at this position
  • tMean — capped mean of normalized IPDs observed at this position
  • tErr — capped standard error of normalized IPDs observed at this position (standard deviation / sqrt(coverage)
  • modelPrediction — normalized mean IPD predicted by the synthetic control model for this sequence context
  • ipdRatio — tMean / modelPrediction
  • coverage — count of valid IPDs at this position (see Filtering section for details)
  • frac — estimate of the fraction of molecules that carry the modification
  • fracLow — 2.5% confidence bound of frac estimate
  • fracUpp — 97.5% confidence bound of frac estimate

That's the full column roster for in-silico control mode, and not a single one is getting skipped. 💖

case-control mode 🎀

okay bestie, when you're running in case-control mode the output is a whole other vibe — every row is one observation and the columns are giving full kinetic tea. here's the whole roster, no skips, no cap:

  • refId — the reference sequence ID of this observation. it's literally telling you which ref sequence we're looking at 🧚
  • tpl — the 1-based template position. (1-based, bestie, not 0-based, don't come for me)
  • strand — the native sample strand where the kinetics were generated. '0' is the strand of the original FASTA, '1' is the opposite strand from FASTA ✨
  • base — the cognate base at this position in the reference. the real one, the matching one 💅
  • score — the Phred-transformed pvalue that a kinetic deviation exists at this position 👀
  • caseMean — the mean of the normalized case IPDs observed at this position
  • controlMean — the mean of the normalized control IPDs observed at this position
  • caseStd — the standard deviation of the case IPDs observed at this position
  • controlStd — the standard deviation of the control IPDs observed at this position
  • ipdRatio — tMean / modelPrediction
  • testStatistic — the t-test statistic 🍓
  • coverage — the mean of the case and control coverage
  • controlCoverage — the count of valid control IPDs at this position (peep the Filtering section for the deets 👑)
  • caseCoverage — the count of valid case IPDs at this position (again, the Filtering section has the details, iykyk)
no seriously, if you wanna know what counts as "valid," the Filtering section is right there. go read it, bestie, I'm not gatekeeping, I'm just not repeating myself 💖

👀 SEE ALSO 👀

and that's the whole page, we ate and left no crumbs 💅✨ this little masterpiece is IPDSUMMARY(1), dated December 2015 — yeah she's an elder but she still slays. mwah, go summarize some modifications, ily 🫶🌸

📦 where this page comes from