Number of reads processed by modkit

[ZhihanNUS]

Hi, it is my first time to use Modkit to extract modified base signals but have some problems.

I use Dorado v0.7.1 to do SUP + modification calling (m6A,pseU).
I use minimap2 to align the unaligned_mod_bam file
Number of read aligned to 'Mus_musculus.GRCm39.dna.primary_assembly.fa' is 19,446,790
However, when I use modkit to extract m6A, it says Processed ~7505886 reads and skipped zero reads

I wonder why modkit filtered the other reads?
Here is what I did,
modkit pileup input.bam output.bed --ref Mus_musculus.GRCm39.dna.primary_assembly.fa --motif A 0 -t 36

I notice it says Threshold of 0.6464844 for base A is low. Consider increasing the filter-percentile or specifying a higher threshold.
Is it related to my problem?

Thank you very much!

[ArtRand]

Hello @ZhihanNUS,

The pileup command will only use primary alignments to generate the bedMethyl table. If there are errors with records in the modBAM file you can inspect the reasons by using the --log-filepath $log_file option. The threshold is not related to the number of reads that are used, skipped, or error. Let me take a look at what threshold values I get on similar (mouse) data and get back to you.

[ZhihanNUS]

Hi @ArtRand ,

I have checked my aligned bam file and the number of reads of primary alignment is 19M, I tried to realign the bam file and run modkit, but I still find the mismatch. Since there is no skipped reads, does it mean modkit pocessed all reads in bamfile, or there is any other rules to filter some reads?
This mismatch between number of reads from aligned bam file and number of reads processed by modkit happened in all samples. Do you have any other suggestions to check the reason for this?

Thank you very much!

[ArtRand]

Hello @ZhihanNUS,

Since there is no skipped reads, does it mean modkit pocessed all reads in bamfile, or there is any other rules to filter some reads?

Can you check the output of the --log-filepath? It will tell you if a read isn't used (usually because there isn't any modified base information).

The counter in modkit pileup that you see on the command line isn't completely perfect (thus why it says "Processed approx X reads"). The information in the log file will be more accurate. If you're not seeing any errors in the log file, all of the primary alignments are being used. If you have a log file I'd be happy to review it.

[ZhihanNUS]

m6Alog.txt
Hello @ArtRand ,
Here it is the log file. And may I know if you have a recommend threshold value for mouse data?

[hoppo]

Hi,
I would also like to piggy-back on this query.
When I run dorado aligner, I get statistics which tell me the following:-
total/primary/unmapped 10256383/4888094/542564

I then run modkit pileup and I get the following numbers (also in the output in --log-filepath output):-
[src/pileup/subcommand.rs::819][2024-07-03 17:12:19][INFO] Done, processed 527371612 rows. Processed ~5681665 reads and skipped ~2651 reads.

I'm wondering where the discrepancy between the '4,888,094' which dorado aligner reports as 'primary' and the 5,684,316 (5,681,665 + 2651) which modkit pileup reports processing and skipping (which should just be the 'primary' alignments).

[ArtRand]

Hello @hoppo

The pileup function doesn't count every read that is skipped. The skipped count are primary alignments that are not used in the aggregation, you could call this "fail count" but some of the reasons aren't exactly failure cases. The read-ids for these cases as well as the skip reason will end up in the log file. Non-primary alignments or records with zero-length sequence are also ignored, but they aren't added to the count.

The reason the counts are approximate is due to how Modkit parallelizes over the genome during pileup. The genome is divided into intervals (100kb by default) and multiple intervals are computed on in parallel. The results are coalesced and emitted. The problem is, a read can span between two intervals:

|---A---|---B---|
    <-----> Read in two intervals

In this case the read may be "processed" or "skipped" during both calculations, resulting in a double counting. You could keep track and not double count - but the extra bookkeeping will incur a decent amount of overhead and I'm not sure it's worth it.

The extract suite of commands don't have this problem, because they work on the read level. The commands will still process intervals of the genome in parallel, but it's clever not to double count (i.e. emit records for a read twice). So instead, what I think is a better solution is to provide a "tag-check" command that tabulates:

1. Number of primary alignments
   1.1 Number of primary alignments that have valid tags and modification calls
   1.2 Number of primary alignments with invalid tags (corrupt or similar)
   1.3 Number of primary alignments that don't have modification information (zero modification calls, things like that)
2. Number of non-primary alignments
   2.1, 2.2, 2.3 (same things but for non-primary alignments).

[ArtRand]

Hello @ZhihanNUS,

The log you sent doesn't have any messages with reads that have errored. The estimated pass threshold for A-modifications (0.6347656) is lower than I might expect. Could you run modkit sample-probs with --hist and see what the distribution looks like?