How to interpret the SMRT Link base modification algorithm output?

Question

I have got PacBio data from a Sequel machine at a coverage > 50X. I have run it in SMRT Link with their Base Modification Detection algorithm. I am now trying to figure out how to interpret the output and where to find the information about the modified bases. One of the outputs is a GFF3 file. The top of this file looks like this:

> ##gff-version 3

...

> I kinModCall  modified_base   1168    1168    31  -   .   coverage=48;context=TTAAACTAATTTCCAAAAATTTTGAACAGTTTTTAATACAA;IPDRatio=2.96
> I kinModCall  modified_base   1497    1497    32  -   .   coverage=50;context=CTTGAAAACACGAATATCTGTTTAGACAATAGTTTTGCCAA;IPDRatio=3.19
> I kinModCall  modified_base   1509    1509    33  -   .   coverage=48;context=CATGAAAAATTTCTTGAAAACACGAATATCTGTTTAGACAA;IPDRatio=2.84
> I kinModCall  modified_base   1529    1529    39  -   .   coverage=46;context=AAGAAAAGACCAAGAAAAACCATGAAAAATTTCTTGAAAAC;IPDRatio=3.84
> I kinModCall  modified_base   1904    1904    36  -   .   coverage=54;context=TAAAACCTGACAGAAAGAATTGGTCCCGTCATCTACTTTCT;IPDRatio=2.86
> I kinModCall  modified_base   2152    2152    31  -   .   coverage=59;context=TTAAAAATTAATATAAAATCAGGGAAATTTTTTTAAATTTT;IPDRatio=2.98
> I kinModCall  modified_base   2226    2226    36  -   .   coverage=63;context=TGGGCGGGGCGTCTCTGCCGGGAAAATTTTGCGTTTAAAAA;IPDRatio=2.59
> I kinModCall  modified_base   2228    2228    34  -   .   coverage=52;context=CGTGGGCGGGGCGTCTCTGCCGGGAAAATTTTGCGTTTAAA;IPDRatio=2.98
> I kinModCall  modified_base   2237    2237    33  -   .   coverage=60;context=GCACGGTTTCGTGGGCGGGGCGTCTCTGCCGGGAAAATTTT;IPDRatio=3.07
> I kinModCall  modified_base   2240    2240    44  -   .   coverage=63;context=GCGGCACGGTTTCGTGGGCGGGGCGTCTCTGCCGGGAAAAT;IPDRatio=2.89
> I kinModCall  modified_base   2311    2311    31  -   .   coverage=63;context=AAAATCGACGAAAATCGGTATAAAATCAAATAAAAATAGAA;IPDRatio=4.35
> I kinModCall  modified_base   2313    2313    39  -   .   coverage=62;context=GAAAAATCGACGAAAATCGGTATAAAATCAAATAAAAATAG;IPDRatio=4.00
> I kinModCall  modified_base   2314    2314    34  -   .   coverage=63;context=AGAAAAATCGACGAAAATCGGTATAAAATCAAATAAAAATA;IPDRatio=2.63
> I kinModCall  modified_base   2333    2333    35  -   .   coverage=62;context=AAAAAAAAGAGAAAAAATGAGAAAAATCGACGAAAATCGGT;IPDRatio=3.77
> I kinModCall  modified_base   2334    2334    35  -   .   coverage=63;context=CAAAAAAAAGAGAAAAAATGAGAAAAATCGACGAAAATCGG;IPDRatio=2.77

I expected something along the lines of 'chromosome I position 1497 methylcytosine', etc.

So it seems that the algorithm has detected modified bases at the following positions but it cannot call them, maybe because of a lack of coverage?

Update: the screenshot below shows there is no CSV file in the Data tab of the completed analysis.

Answer 1

It sounds like you're looking for the modifications.csv file, rather than the modifications.gff file. The .gff file is for viewing in SMRT View (or other viewers). You can see schema for both files at PacBio's Methylome Analysis Technical Note. Specifically, modified_base means the algorithm could not determine the exact base and modification.

I also found a video tutorial for performing the analysis in SMRTLink on YouTube

EDIT: Per PacBio's Methylome Analysis Technical Note documentation, the .csv "...is also produced when motif analysis is not active, such as with the RS_Modification_Detection analysis protocol." Not sure why it's missing, perhaps because the algorithm couldn't clearly identify any of the modifications. You'll probably need to contact tech support on this issue.

In either case, the .gff you posted includes the significant hits (p < 0.01), but the algorithm was unable to classify them, likely because of insufficient coverage. PacBio's Methylome Analysis Technical Note documentation gives specific coverage recommendations for each event type, but 100x is the lower bound. Lots of good information there.

Probably not what you wanted to hear, but I hope that gets you going.