High 6mA threshold in epi2me chromatin accessibility dataset

[billytcl]

I’m taking a look at the chromatin accessibility dataset from ONT’s epi2me labs. I downloaded the raw pod5 files and called it with dorado using sip and 5mCG_5hmCG,6mA.

Using modkit 0.6.4’s modkit summary as an overall QC, I saw the following thresholds:

  # bases                     C,A                         
  # total_reads_used          24902062                    
  # count_reads_A             24902062                    
  # count_reads_C             24902062                    
  # pass_threshold_A          0.9980469                   
  # pass_threshold_C          0.7910156
  # modification_codes_for_A  a
  # modification_codes_for_C  h,m
base  code  pass_count   pass_frac    all_count
   all_frac
   A     -     61054750491  0.9983692    63991911696
   0.99128604
   A     a     99731455     0.001630812  562521620
   0.008713911 
   C     -     809455209    0.38108593   873978594
   0.37053645
   C     h     30801412     0.014501093  89012849
   0.03773834
   C     m     1283818670   0.604413     1395693163
   0.59172523

Is the A threshold supposed to be that high? I seem to recall that was a bug in older modkit versions for 5mC all context.

[ArtRand]

@billytcl sorry about the delay, which sample is this? Here's a summary from "PBA15156", one of the CS runs.

❯ modkit --version
modkit 0.6.4
❯ modkit mb summary PBA15156/basecalls.bam -t 32
> not subsampling, using all reads
> calculating threshold that removes lowest 10(th) percentile of calls
> collecting base modification calls from unmapped records
# bases                     C,A
# total_reads_used          11428573
# count_reads_A             11428573
# count_reads_C             11428573
# pass_threshold_A          0.9628906
# pass_threshold_C          0.7519531
# modification_codes_for_A  a
# modification_codes_for_C  h,m
 base  code  pass_count   pass_frac     all_count    all_frac
 A     -     28546627881  0.99594784    31254250146  0.9867772
 A     a     116145325    0.0040521314  418807467    0.01322283
 C     -     373728036    0.36350924    407770235    0.35716933
 C     h     24156399     0.023495896   57899376     0.05071454
 C     m     630227023    0.6129948     676002501    0.5921162

[billytcl]

No problem! I rebasecalled and merged the two flowcells from epi2me labs open data - PAY22766 and PBA15156. As far as I know, the dataset on AWS was done with hac v5.2.0 and modkit v0.5 which is pretty outdated given the various improvements. We subsequently used dorado v1.4.0 to rebasecall on sup v5.2.0 and modkit v0.6.4. This is the relevant header of the bam:

@RG     ID:0c6d9c31-dd4a-4ec8-b367-b67bad30254e_dna_r10.4.1_e8.2_400bps_sup@v5.2.0      PU:PAY22766     PM:PC48B260     DT:2025-03-11T12:23:30.820000+00:00     PL:ONT  DS:runid=0c6d9c31-dd4a-4ec8-b367-b67bad30254e basecall_model=dna_r10.4.1_e8.2_400bps_sup@v5.2.0 modbase_models=dna_r10.4.1_e8.2_400bps_sup@v5.2.0_5mCG_5hmCG@v2,dna_r10.4.1_e8.2_400bps_sup@v5.2.0_6mA@v1 experiment_id=20250311_PS_HG002_PLusMinus_ChromStenc acquisition_start_time=2025-03-11T12:25:23.543000+00:00 model_stride=6   LB:ONLL09557_Plus_CS_C
@RG     ID:e690665b-d836-49f0-bd83-f90610ee90f0_dna_r10.4.1_e8.2_400bps_sup@v5.2.0      PU:PBA15156     PM:PC48B260     DT:2025-03-11T12:23:30.876000+00:00     PL:ONT  DS:runid=e690665b-d836-49f0-bd83-f90610ee90f0 basecall_model=dna_r10.4.1_e8.2_400bps_sup@v5.2.0 modbase_models=dna_r10.4.1_e8.2_400bps_sup@v5.2.0_5mCG_5hmCG@v2,dna_r10.4.1_e8.2_400bps_sup@v5.2.0_6mA@v1 experiment_id=20250311_PS_HG002_PLusMinus_ChromStenc acquisition_start_time=2025-03-11T12:25:40.717000+00:00 model_stride=6   LB:ONLL09555_Plus_CS_A
@PG     ID:basecaller   PN:dorado       VN:1.4.0+ba44a013       CL:dorado "basecaller" "sup,5mCG_5hmCG,6mA" "chunk_input" "--no-trim"   DS:gpu:NVIDIA H100 80GB HBM3

Modkit summary:

# bases                     C,A 
# total_reads_used          24902062 
# count_reads_A             24902062 
# count_reads_C             24902062 
# pass_threshold_A          0.9980469 
# pass_threshold_C          0.7910156 
# modification_codes_for_A  a 
# modification_codes_for_C  h,m 
 base  code  pass_count   pass_frac    all_count    all_frac 
 A     -     61054750491  0.9983692    63991911696  0.99128604 
 A     a     99731455     0.001630812  562521620    0.008713911 
 C     -     809455209    0.38108593   873978594    0.37053645 
 C     h     30801412     0.014501093  89012849     0.03773834 
 C     m     1283818670   0.604413     1395693163   0.59172523 

As you can see here, it looks like the pass_threshold for A is 0.998. It seems really high!

[ArtRand]

Hello @billytcl,

Thanks for bringing this up. I'm in the process of reproducing the basecalls like you have. In the meantime I've checked on a few thousand reads and looked at both modbam summary and the first few log lines of pileup.

❯ modkit --version
modkit 0.6.4

Summary:

❯ modkit mb summary test_basecalls_sort.bam --matched-only -r ${HG38} --threads 64
> not subsampling, using all reads
> calculating threshold that removes lowest 10(th) percentile of calls
> calculating summary at matches to primary sequence bases
# bases                     A,C
# total_reads_used          4628
# count_reads_A             4628
# count_reads_C             4628
# pass_threshold_A          0.9980469
# pass_threshold_C          0.8105469
# modification_codes_for_A  a
# modification_codes_for_C  h,m
 base  code  pass_count  pass_frac     all_count  all_frac
 A     -     10412446    0.9985748     10870189   0.99252224
 A     a     14861       0.0014252001  81897      0.0074777533
 C     -     127824      0.37254307    137772     0.36202723
 C     h     4073        0.01187076    12295      0.032307908
 C     m     211215      0.61558616    230490     0.60566485

I'm seeing the 10th percentile cutoff be calculated as the highest quality bin. With default settings, dorado basecaller will elide modification calls when $P(\text{modified}) \gt 0.95$ and Modkit interprets these calls as having $\text{mod qual} = 0.0 == P(\text{unmodified}) = 0.998$. You can disable this behavior by adding --modified-bases-threshold 0.0 to the dorado basecaller command.

The CS samples we uploaded to Epi2Me intentionally have very low background 6mA at heterochromatic regions. What is happening here is that all of the 6mA sites are in the lowest 10th percentile. The pileup function has a check for this:

❯ modkit pileup \
∙ test_basecalls_sort.bam \
∙ test_pileup.bed \
∙ --modified-bases 6mA 5mC 5hmC \
∙ --ref ${HG38}
> discarded 121 contigs with zero aligned reads
> parsed 3 base modification(s). Base modifications other than 'A:a,C:h,C:m' will be counted as 'N_other'.
> adding single-base motif: 'A 0'
> adding single-base motif: 'C 0'
> attempting to sample at least 10042 reads
> calculating threshold value with probabilities aligned to reference base(s) A,C
> estimated threshold for A too high (0.9980469, qual: 255), setting threshold 0.9511719 (qual: 243), highest discovered explicit canonical value

The last line shows how the threshold will be adjusted.

I think it would be good for modkit mb summary to have a similar log and/or adjust the threshold value.

For what it's worth, the chromatin stenciling functionality in Modkit is perfectly happy with HAC basecalls. Whilst developing it we focused on HAC. SUP should work fine as well, and I'll take a closer look to make sure that everything is working with the setup you have here.

[billytcl]

Awesome! So an interim workaround would be to use the pileup correction rather than reading the summary?

…

On Mon, Jul 13, 2026 at 10:10 AM Arthur Rand ***@***.***> wrote: *ArtRand* left a comment (nanoporetech/modkit#639) <#639 (comment)> Hello @billytcl <https://github.com/billytcl>, Thanks for bringing this up. I'm in the process of reproducing the basecalls like you have. In the meantime I've checked on a few thousand reads and looked at both modbam summary and the first few log lines of pileup. ❯ modkit --version modkit 0.6.4 Summary: ❯ modkit mb summary test_basecalls_sort.bam --matched-only -r ${HG38} --threads 64> not subsampling, using all reads> calculating threshold that removes lowest 10(th) percentile of calls> calculating summary at matches to primary sequence bases# bases A,C# total_reads_used 4628# count_reads_A 4628# count_reads_C 4628# pass_threshold_A 0.9980469# pass_threshold_C 0.8105469# modification_codes_for_A a# modification_codes_for_C h,m base code pass_count pass_frac all_count all_frac A - 10412446 0.9985748 10870189 0.99252224 A a 14861 0.0014252001 81897 0.0074777533 C - 127824 0.37254307 137772 0.36202723 C h 4073 0.01187076 12295 0.032307908 C m 211215 0.61558616 230490 0.60566485 I'm seeing the 10th percentile cutoff be calculated as the highest quality bin. With default settings, dorado basecaller will elide modification calls when $P(\text{modified}) \gt 0.95$ and Modkit interprets these calls as having $\text{mod qual} = 0.0 == P(\text{unmodified}) = 0.998$. You can disable this behavior by adding --modified-bases-threshold 0.0 to the dorado basecaller command. The CS samples we uploaded to Epi2Me intentionally have very low background 6mA at heterochromatic regions. What is happening here is that *all* of the 6mA sites are in the lowest 10th percentile. The pileup function has a check for this: ❯ modkit pileup \ ∙ test_basecalls_sort.bam \ ∙ test_pileup.bed \ ∙ --modified-bases 6mA 5mC 5hmC \ ∙ --ref ${HG38}> discarded 121 contigs with zero aligned reads> parsed 3 base modification(s). Base modifications other than 'A:a,C:h,C:m' will be counted as 'N_other'.> adding single-base motif: 'A 0'> adding single-base motif: 'C 0'> attempting to sample at least 10042 reads> calculating threshold value with probabilities aligned to reference base(s) A,C> estimated threshold for A too high (0.9980469, qual: 255), setting threshold 0.9511719 (qual: 243), highest discovered explicit canonical value The last line shows how the threshold will be adjusted. I think it would be good for modkit mb summary to have a similar log and/or adjust the threshold value. For what it's worth, the chromatin stenciling functionality in Modkit is perfectly happy with HAC basecalls. Whilst developing it we focused on HAC. SUP should work fine as well, and I'll take a closer look to make sure that everything is working with the setup you have here. — Reply to this email directly, view it on GitHub <#639?email_source=notifications&email_token=ACPHYT4YWBUULUW72CBOTGL5EUJX7A5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTIOJWGA2TQNBSGM4KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-4960584238>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACPHYT4YBIR3357ZJ4FYKLT5EUJX7AVCNFSNUABFKJSXA33TNF2G64TZHM3DAOJSGI2DKMRWHNEXG43VMU5TINZYGI4TOMRVGQ22C5QC> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[ArtRand]

Something like that. What are you trying to do? Maybe I can come up with something more streamlined until the next version is out.

[billytcl]

We are trying to model the relationship between CG methylation and the associated open chromatin marked with 6mA. Originally our tool just called modkit summary to parse the 6mA threshold, and then parse the reads in the bams (we can also use call-mods to make our lives easier but it's not that fast and takes up storage space...). But I guess we can just parse out the pileup log! It's a little more painful to read it out but we only need to do it once for the HG002 public dataset...

…

On Mon, Jul 13, 2026 at 3:05 PM Arthur Rand ***@***.***> wrote: *ArtRand* left a comment (nanoporetech/modkit#639) <#639 (comment)> Something like that. What are you trying to do? Maybe I can come up with something more streamlined until the next version is out. — Reply to this email directly, view it on GitHub <#639?email_source=notifications&email_token=ACPHYT57IFEHRMNABYH5RWL5EVMKHA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTIOJWGMZDAMBQGE3KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-4963200016>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACPHYT7DFRE7FGA4IFZHWMD5EVMKHAVCNFSNUABFKJSXA33TNF2G64TZHM3DAOJSGI2DKMRWHNEXG43VMU5TINZYGI4TOMRVGQ22C5QC> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[ArtRand]

Here's a build where modkit modbam summary has the same threshold behavior as pileup. There is a branch with the changes as well. I'll get these merged into the next release - but if you want to try this in your pipeline the CLI won't change in the next version.

There is also a new line in the output of modkit modbam check-tags that lists the modified bases so you can paste this line directly into pileup - but that's not directly related to what you're doing here.

modkit_u16_x86_64.tar.gz

[billytcl]

Thank you! Will give it a shot.

…

On Tue, Jul 14, 2026 at 5:49 PM Arthur Rand ***@***.***> wrote: *ArtRand* left a comment (nanoporetech/modkit#639) <#639 (comment)> Here's a build where modkit modbam summary has the same threshold behavior as pileup. There is a branch <https://github.com/nanoporetech/modkit/tree/ar/qol-list-mods-and-threshold> with the changes as well. I'll get these merged into the next release - but if you want to try this in your pipeline the CLI won't change in the next version. There is also a new line in modkit modbam check-tags that lists the modified bases so you can paste this line directly into pileup - but that's not directly related to what you're doing here. modkit_u16_x86_64.tar.gz <https://github.com/user-attachments/files/30028885/modkit_u16_x86_64.tar.gz> — Reply to this email directly, view it on GitHub <#639?email_source=notifications&email_token=ACPHYT6KPUUAENEHBUAZSRD5E3IIJA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTIOJXGU3DANRVGUYKM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-4975606550>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACPHYT4PIJC6WD33CKOIAFD5E3IIJAVCNFSNUABFKJSXA33TNF2G64TZHM3DAOJSGI2DKMRWHNEXG43VMU5TINZYGI4TOMRVGQ22C5QC> . You are receiving this because you were mentioned.Message ID: ***@***.***>