kmer-model format

[jsthv]

I am trying to simulate signal data by creating a kmer-model for a modified base. My data comes from R10.4.1 flow cell and thus needs a 9-mer kmer model. I have a sequence with a site-specific modification that causes a large decrease in the current. All I want to do is try to simulate the signal by only varying the parameters for the 9-mers that overlap the modification. It would have been very straightforward if I could just take the ont R10.4.1 400 bps 2 column file and simply add 9 lines to the file at the end for only the 9 modified 9mers in my sequence. It seems that I cannot do that or that I am not doing it correctly. The manual says that I have to use a file formatted in the same was as for f5c (r9.4_450bps.nucleotide.6mer.template.model). When I look at that file, it is a 6mer kmer library with a 6 column structure. The file for r9.4_450bps.cpg.6mer.template.model has a 5 column structure. There are no examples of 9mer libraries. Is the problem that I would need a 5^9 line file for A, C, G, T, X which would be unwieldly? Is that why I have to use a 6mer library, which translates to 5^6 lines? Strangely, the errors I get seem to indicate that it could hand 9mer libraries but I don't know how to get it to do so.

For example I get the error:
[read_model::INFO] k-mer size in file /home/jst/models/9mer_XX_CPD_f5c.txt is 9
[read_model::ERROR] File /home/jst/models/9mer_XX_CPD_f5c.txt has too many entries. Expected 262144 kmers in the model, but file had more than that At src/model.c:114

This is what I had in the header:

#model_name R10.4.1_400bps_9mer_custom
#kit R10.4.1_400bps
#strand template
#k 9
#alphabet dnaX

If it did recognize it as a 9mer library with a modified base X, then it should have anticipated 5^9 9-mers, not 4^9 which is 262144.

Is there anyway to simply use the two column ONT 9mer library with a few added modified 9-mers? My case is a bit non standard because my modifications run in tandem, i.e., XX.

[hasindu2008]

This --kmer-model and --meth-model were developed mostly for debugging and development. So it hasn't been thoroughly documented.

Here is an example 9-mer model that will work with --kmer-model. If you only require the level_mean, you can make level_stdv, sd_mean and sd_stdv all 1.

sq0.3.0_r10.4_400bps.nucleotide.9mer.template.model.txt

But seems like what you are after is something with a modified base. So the closest model to use would be the methylation model. Here is a methylation model that works with the --meth-model.
https://unsw-my.sharepoint.com/:u:/g/personal/z5136909_ad_unsw_edu_au/EXr0hnzSpJNIj_8MSq43tiQB_MN94ZRUEw3NMH-SEvUxVA?e=9FlAq0

But squigulator as it is would not be compatible with an 'X' base or the 2-column ONT format though. Supporting an 'X' base or the 2-column model format would need some edits to the code.

[jsthv]

Sounded like a great solution. I just used M as my modification and changed the mean current value for the 9-mers involving them in my sequence. Unfortunately I got back an error saying that your program does not recognize the --meth-model.:

./squigulator BC03_CPD_290mer.fa --meth-model /home/jst/models/r10_400bps.cpg.9mer.template.model -o BC03_290mer_CPD_signal.slow5 -c BC03_290mer_CPD_align.paf -q BC03_290mer_read.fa --ideal --full-contigs
./squigulator: unrecognized option '--meth-model'

Has this option been removed? It might solve my problem if it were reinstated.

When I tried using the --kmer-model option it gave back:

./squigulator BC03_CPD_290mer.fa --kmer-model /home/jst/models/r10_400bps.cpg.9mer.template.model -o BC03_290mer_CPD_signal.slow5 -c BC03_290mer_CPD_align.paf -q BC03_290mer_read.fa --ideal --full-contigs
[INFO] sim_main: Using random seed: 1754604504
[read_model::WARNING] Reading the model from file /home/jst/models/r10_400bps.cpg.9mer.template.model. This is an experimental feature. Use with caution At src/model.c:44
[read_model::INFO] k-mer size in file /home/jst/models/r10_400bps.cpg.9mer.template.model is 9
[read_model::ERROR] File /home/jst/models/r10_400bps.cpg.9mer.template.model has too many entries. Expected 262144 kmers in the model, but file had more than that At src/model.c:114

I assume because it was expecting 4^9 lines, not 5^9 lines as it assumes a canonical set of bases.

[hasindu2008]

May I know if you are using the latest version of squigulator?

[jsthv]

I thought that I had posted that I figured out I was using an older version (v0.2.2). I have since been playing with the 400 kHz cpg R10 library you pointed me to in comparison with the standard 500 kHz R10_min library for canonical sequences. I find that the R10_min library gives a better match to the actual data. For some reason which I don't understand, when using -x R10_min I get the expected match between the sequence and the signal with a -6 base shift and the 9-mers shown align properly with the expectation that the highest signals are for T's in the 7th position, but when using the --meth-model option, I have to use a -2 offset and the 9-mers don't correspond to T's in the 7 the 7th position. I have attached a file with what I have observed. It would be great if there was a cpg kmer R10 kmer library that reproduces the R10_min predictions.

Analysis of barcode 01.pdf

[jsthv]

Another problem. I ran squigulator using a modified cpg 9mer file with a new name in which I replaced the level mean value of 89 with 10 for a M base is specific 9-mers around my modification site but I got the same signal profile as I got using the unmodified cpg file. Is squigulator simply grabbing a default cpg kmer file when the option --methyl-model is invoked and ignoring the filename I am giving it? --meth-model option is working because if I replace the M in the sequence with T, I get the expected signal increase.

[jsthv]

On reviewing the screen output (see below) when using squigulator v0.4.0 and the --meth-model option with my edited file, it reports a "builtin DNA R9 nucleotide model loaded" instead which must be is why I can't over-ride it with my own file with the --meth-model option? When I run -x dna-r10-min it says "builtin DNA R10 model loaded"

(base) jst@taylor-desktop:~/squigulator-v0.4.0$ ./squigulator BC01_290mer.fa --meth-model /home/jst/models/r10_400bps.cpg.9mer.template.model -o BC01_290mer_cpg_signal.slow5 -c BC01_290mer_cpg_align.paf -q BC01_cpg_290mer_read.fa --ideal --full-contigs
[INFO] sim_main: Using random seed: 1754683540
[init_core::INFO] builtin DNA R9 nucleotide model loaded
[INFO] load_ref: Loaded 1 reference sequences with total length 0.000290 Mbases
[INFO] print_model_stat: digitisation: 2048.0; sample_rate: 4000.0; range: 748.6; offset_mean: -237.4; offset_std: 14.2; dwell_mean: 9.0; dwell_std: 4.0

How can I get my R10 model loaded?

[jsthv]

I have played around with it some more, and find that if I use the -x dna-r10-min option without the --meth-model and run it on a sequence with M in it, it will calculate a signal for the 5mC that is very different from T. I get the same signal trace whether or not I additionally add the --meth-model option with a modified file.

[hasindu2008]

I thought that I had posted that I figured out I was using an older version (v0.2.2). I have since been playing with the 400 kHz cpg R10 library you pointed me to in comparison with the standard 500 kHz R10_min library for canonical sequences.

• Did you mean 4kHz and 5KHz?

The problem you are facing could be this. Sorry, I forgot to mention this before. The methylation file I gave you is intended to run like this.

 squigulator -x dna-r10-min -o a.slow5  test/nCoV-2019.reference.fasta --meth-model r10_400bps.cpg.9mer.template.model --meth-freq test/mfreq.tsv 

mfreq should be like this https://github.com/hasindu2008/squigulator/blob/dev/test/mfreq.tsv. This tsv file should have three columns, chr, 0-based pos, and methylation frequency.

I have played around with it some more, and find that if I use the -x dna-r10-min option without the --meth-model and run it on a sequence with M in it, it will calculate a signal for the 5mC that is very different from T.

If you mean change the fasta file to have an M? That is not suitable, as M in fasta is an IUPAC ambiguous base that is A or C [https://genome.ucsc.edu/goldenPath/help/iupac.html], which if I remember right, squigulator will just replace with an 'A'.

[jsthv]

Thanks for all your help but it still doesn't seem to solve the problem. It still seems that it is ignoring my modified cpg file though it acknowledges opening it. As per your directions (though I also added the --sample-rate option to match the cpg file rate) I have put C's into the 145 mer sequence at 70 and 71 which I duplicated to make a 290mer. At the second identical site (216/217) I either inserted two T's or two C's. In the tsv file I only identified the 70 and 71 C's as being methylated 100%. When I ran squigulator v0.4.0 with the command line (below) using my file where I changed the signal level for the two 9mers as follows: CAAGAAMMA and AAGAAMMAG corresponding to 9mers containing C70 and C71 which positions the modifications at the 7th position of the 9mer where it should make the greatest contribution to the signal to 19.9999 each vs original values of 109.52 for CAAGAACCA and 91.8 for AAGAACCAG and 108.3 for CAAGAAMMA and 89.6 for AAGAAMMAG. What I find is that the methylated CC site at 70/71 gives the identical signal as the canonical CC site at 216/217 as I get using the unmodified cpg model file, indicating that is is not using the modified cpg file and only using the dna-r10-min file. I know that the tsv file must be working because if I replace the two C's with T's in the reference sequence at 70/71 it reports an error that it didn't find a C at 70.

(base) jst@taylor-desktop:~/squigulator-v0.4.0$ ./squigulator BC03_CPD_CC_290mer.fa -x dna-r10-min -o BC03_CPD_CC_290mer_signal5.slow5 -c BC03_CPD_CC_290mer_align5.paf -q BC03_CPD_CC_290mer_read5.fa --sample-rate 4000 --meth-model /home/jst/models/r10_400bps.cpg.9mer.CPD_20.model --meth-freq BC03_CPD_CC_290mer.tsv --ideal --full-contigs
[set_profile::INFO] dna-r10-min is 5kHz from squigulator v0.3.0 onwards. Specify --sample-rate 4000 for old 4kHz.
[set_profile::WARNING] Parameters and models for dna-r10-min 5khz are still crude. If you have good modification-free data, please share! At src/sim.c:170
[sim_main::WARNING] Option --meth-freq is experimental. Please report any issues. At src/sim.c:995
[INFO] sim_main: Using random seed: 1754838710
[INFO] sim_main: bps=400.0,sample_rate=4000.0,set dwell_mean=sample_rate/bps=10.0
[init_core::INFO] builtin DNA R10 nucleotide model loaded
[read_model::WARNING] Reading the model from file /home/jst/models/r10_400bps.cpg.9mer.CPD_20.model. This is an experimental feature. Use with caution At src/model.c:46
[read_model::INFO] k-mer size in file /home/jst/models/r10_400bps.cpg.9mer.CPD_20.model is 9
[INFO] load_ref: Loaded 1 reference sequences with total length 0.000290 Mbases
[INFO] print_model_stat: digitisation: 8192.0; sample_rate: 4000.0; range: 1536.6; offset_mean: 13.4; offset_std: 16.3; dwell_mean: 10.0; dwell_std: 4.0
[INFO] sim_main: 1/1 reads done
[main] Version: 0.4.0
[main] CMD: ./squigulator -x dna-r10-min -o BC03_CPD_CC_290mer_signal5.slow5 -c BC03_CPD_CC_290mer_align5.paf -q BC03_CPD_CC_290mer_read5.fa --sample-rate 4000 --meth-model /home/jst/models/r10_400bps.cpg.9mer.CPD_20.model --meth-freq BC03_CPD_CC_290mer.tsv --ideal --full-contigs BC03_CPD_CC_290mer.fa
[main] Real time: 1.167 sec; CPU time: 1.168 sec; Peak RAM: 0.604 GB

[hasindu2008]

Ohh! Sorry for the trouble. This seems to be a bug!
Please give me a bit of time to fix it.

[hasindu2008]

@jsthv
Would you be able to compile the branch https://github.com/hasindu2008/squigulator/tree/noncpg and test if you now see that the signal values are actually coming from the model file?
I did not thoroughly test, but a simple test with CCCCCCCCC and MMMMMMMMM was working.

[jsthv]

Compiled and working! Thank you! Very helpful!! The only minor visualization problem is that you used black to identify the modification site, but that obscures the signal so one cannot see where it goes: Perhaps no color or light grey or orange? Unless there is someway for me to change it. ***@***.*** From: Hasindu Gamaarachchi ***@***.***> Sent: Tuesday, August 26, 2025 12:52 AM To: hasindu2008/squigulator ***@***.***> Cc: Taylor, John-Stephen ***@***.***>; Mention ***@***.***> Subject: Re: [hasindu2008/squigulator] kmer-model format (Issue #29) [https://avatars.githubusercontent.com/u/12987163?s=20&v=4]hasindu2008 left a comment (hasindu2008/squigulator#29)<#29 (comment)> @jsthv<https://github.com/jsthv> Would you be able to compile the branch https://github.com/hasindu2008/squigulator/tree/noncpg and test if you now see that the signal values are actually coming from the model file? I did not thoroughly test, but a simple test with CCCCCCCCC and MMMMMMMMM was working. - Reply to this email directly, view it on GitHub<#29 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/BEVIT5M2RU6DHRZRJEG3DIL3PPYWXAVCNFSM6AAAAACDGKE52SVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTEMRSG4YTMMJUGQ>. You are receiving this because you were mentioned.Message ID: ***@***.******@***.***>>

[hasindu2008]

Hi,

Squigulator doesn't do any visualization - could you please explain a bit more?

[jsthv]

That was my mistake. The issue was with squigualiser which I used to display the signal alignments.

One final question. How do I cite the methylation 9mer model that you suggested I use and provided a link to?