Modkit extract calls: 5hmC ("h") calls correspond mostly to <50% Mod_qual from modkit extract full

[ahong017]

Versions
modkit: 0.5.1
Dorado: 1.0.2+c758d2f6

modkit commands
modkit extract calls
modkit extract full

Background

Hello,
I first used Dorado to basecall ONT data with 5mC/5hmC modified-base detection. I then aligned the resulting BAM and ran both modkit extract calls and modkit extract full.
I used modkit extract calls to obtain the discrete base-modification call (call_code) and modkit extract full to obtain the corresponding modification probability (Mod_qual).

Dorado modified-base basecalling
./dorado basecaller \ dna_r10.4.1_e8.2_400bps_hac@v5.2.0 \ --modified-bases-models dna_r10.4.1_e8.2_400bps_hac@v5.0.0_5mC_5hmC@v3 \ R10.4.1_ONT_data.pod5 \ --emit-moves > dorado_modified_bases_basecalled_data.bam

Alignment
./dorado aligner \ GCF_000001405.26_GRCh38_genomic.mmi \ dorado_modified_bases_basecalled_data.bam \ -o dorado_modified_bases_basecalled_data_aligned.bam

modkit extract calls
./modkit extract calls \ --mapped \ dorado_modified_bases_basecalled_data_aligned.bam \ modkit_extract_calls_output.tsv

modkit extract full
./modkit extract full \ --mapped \ dorado_modified_bases_basecalled_data_aligned.bam \ modkit_extract_full_output.tsv

For modkit extract full, I multiplied the modification probability value by 100 so that it was represented on a 0–100 scale. I then restricted both output files to reads overlapping CpG sites that were classified as either 100% methylated or 0% methylated in both my EM-seq and WGBS ground-truth datasets. Next, I matched corresponding rows between the modkit extract calls and modkit extract full outputs using:

• read ID
• Forward_read_position
• Ref_position

For each call_code in the modkit extract calls output (h, m, or -), I counted the corresponding modification probability values from modkit extract full in 10% bins.

Issue

I then observed the following pattern:

• "-" calls were primarily associated with values below 50%.
• "m" calls were primarily associated with values above 50%.
  - However, "h" calls, representing 5hmC, were associated almost entirely with values below 50%.
  This was unexpected because I had assumed that an h call would correspond to a relatively high probability of 5hmC. The distribution is shown in the screenshot below.

My understanding is that modkit extract calls determines the per-read discrete modification call using thresholds derived from the input data rather than necessarily using a fixed 50% cutoff.

Could you clarify how the call_code is determined when both 5mC and 5hmC are present?

In particular:

1. Is Mod_qual in modkit extract full the probability associated specifically with the modification represented by call_code, or does its interpretation depend on the modification probabilities encoded in the MM/ML tags?
2. Is it expected for an h call to have a Mod_qual corresponding to <50% on this scale?
3. Is there a particular threshold or method you would recommend if I want to classify individual reads as modified versus unmodified when using a Dorado model that predicts both 5mC and 5hmC?

Here is the command I used to get the table shown in the screenshot above:

awk -F'\t' -v OFS='\t' -v gt_label="100%" '
function get_bin(value) {
    if (value < 0 || value > 100) return 0
    if (value <= 10) return 1
    if (value <= 20) return 2
    if (value <= 30) return 3
    if (value <= 40) return 4
    if (value <= 50) return 5
    if (value <= 60) return 6
    if (value <= 70) return 7
    if (value <= 80) return 8
    if (value <= 90) return 9
    return 10
}
BEGIN {
    bin_label[1]="0-10"; bin_label[2]="11-20"; bin_label[3]="21-30"
    bin_label[4]="31-40"; bin_label[5]="41-50"; bin_label[6]="51-60"
    bin_label[7]="61-70"; bin_label[8]="71-80"; bin_label[9]="81-90"
    bin_label[10]="91-100"
    call_order[1]="h"; call_order[2]="m"; call_order[3]="-"
}
NR == FNR {
    key=$1 SUBSEP $2 SUBSEP $3
    value=$13
    if (value == "" || value !~ /^[-+]?[0-9]*\.?[0-9]+$/) {
        invalid_percent_rows++
        next
    }
    percent_by_key[key]=value+0
    next
}
{
    key=$1 SUBSEP $2 SUBSEP $3
    call=$15
    if (!(key in percent_by_key)) {
        unmatched_yesno_rows++
        next
    }
    if (call != "h" && call != "m" && call != "-") {
        invalid_call_rows++
        next
    }
    value=percent_by_key[key]
    bin=get_bin(value)
    if (bin == 0) {
        out_of_range_rows++
        next
    }
    count[call,bin]++
    matched_rows++
}
END {
    for (c=1; c<=3; c++)
        for (b=1; b<=10; b++)
            print "Dorado",gt_label,call_order[c],bin_label[b],count[call_order[c],b]+0

    print "Dataset: " gt_label > "/dev/stderr"
    print "Matched rows: " matched_rows+0 > "/dev/stderr"
    print "Unmatched yes/no rows: " unmatched_yesno_rows+0 > "/dev/stderr"
    print "Invalid yes/no calls: " invalid_call_rows+0 > "/dev/stderr"
    print "Invalid percentage rows: " invalid_percent_rows+0 > "/dev/stderr"
    print "Out-of-range percentage rows: " out_of_range_rows+0 > "/dev/stderr"
    print "" > "/dev/stderr"
}' \
TR_dorado_per_read_actual_methy_%_predictions_from_modkit_extract_full_100%_ground_truth_merged.tsv \
TR_dorado_per_read_yes_no_predictions_from_modkit_extract_calls_100%_ground_truth_merged.tsv \
>> TR_dorado_h_m_dash_counts_by_methylation_percent_bin.tsv

awk -F'\t' -v OFS='\t' -v gt_label="0%" '
function get_bin(value) {
    if (value < 0 || value > 100) return 0
    if (value <= 10) return 1
    if (value <= 20) return 2
    if (value <= 30) return 3
    if (value <= 40) return 4
    if (value <= 50) return 5
    if (value <= 60) return 6
    if (value <= 70) return 7
    if (value <= 80) return 8
    if (value <= 90) return 9
    return 10
}
BEGIN {
    bin_label[1]="0-10"; bin_label[2]="11-20"; bin_label[3]="21-30"
    bin_label[4]="31-40"; bin_label[5]="41-50"; bin_label[6]="51-60"
    bin_label[7]="61-70"; bin_label[8]="71-80"; bin_label[9]="81-90"
    bin_label[10]="91-100"
    call_order[1]="h"; call_order[2]="m"; call_order[3]="-"
}
NR == FNR {
    key=$1 SUBSEP $2 SUBSEP $3
    value=$13
    if (value == "" || value !~ /^[-+]?[0-9]*\.?[0-9]+$/) {
        invalid_percent_rows++
        next
    }
    percent_by_key[key]=value+0
    next
}
{
    key=$1 SUBSEP $2 SUBSEP $3
    call=$15
    if (!(key in percent_by_key)) {
        unmatched_yesno_rows++
        next
    }
    if (call != "h" && call != "m" && call != "-") {
        invalid_call_rows++
        next
    }
    value=percent_by_key[key]
    bin=get_bin(value)
    if (bin == 0) {
        out_of_range_rows++
        next
    }
    count[call,bin]++
    matched_rows++
}
END {
    for (c=1; c<=3; c++)
        for (b=1; b<=10; b++)
            print "Dorado",gt_label,call_order[c],bin_label[b],count[call_order[c],b]+0

    print "Dataset: " gt_label > "/dev/stderr"
    print "Matched rows: " matched_rows+0 > "/dev/stderr"
    print "Unmatched yes/no rows: " unmatched_yesno_rows+0 > "/dev/stderr"
    print "Invalid yes/no calls: " invalid_call_rows+0 > "/dev/stderr"
    print "Invalid percentage rows: " invalid_percent_rows+0 > "/dev/stderr"
    print "Out-of-range percentage rows: " out_of_range_rows+0 > "/dev/stderr"
    print "" > "/dev/stderr"
}' \
TR_dorado_per_read_actual_methy_%_predictions_from_modkit_extract_full_0%_ground_truth_merged.tsv \
TR_dorado_per_read_yes_no_predictions_from_modkit_extract_calls_0%_ground_truth_merged.tsv \
>> TR_dorado_h_m_dash_counts_by_methylation_percent_bin.tsv

echo "Done. Results written to:"
echo "TR_dorado_h_m_dash_counts_by_methylation_percent_bin.tsv"
echo
column -t -s $'\t' TR_dorado_h_m_dash_counts_by_methylation_percent_bin.tsv