wgatools call error

[rejo27]

Hi, Weijie
This is a very useful tool !
I get an error when running the command wgatools call a.maf -s -l0 -o a.vcf.

2024-09-03T12:06:49.010920177+08:00 WARN maf index not found, will not generate contig info
thread 'main' panicked at src/tools/caller.rs:351:50:
byte index 780 is out of bounds of `TTTCTATTATTTATGTGTTTTTGACGAGCAACTGATTTTCTTTTTTAATAGGAAATGGTACGCTTGATGAAGGAAAAACGTCTTCCAGCCGCATTCAAGTGTTACCATAATTTTCATAGAGCTGAGGACATATCACCTGATAACCTTTTTTATAAATTGGTTGTCCATGTGCACAGCGATTCAGGATTTAGGCGTTACCACAAAGAGATGAGGTTATCTTTAAACACGTCTTAGTTACTTATGTCTTTCAGTTTGA`[...]

I wrote the sequence in maf format into a separate file, and the same error occurred. Do you have any good suggestions?
This is my maf file

##maf version=1


a
s	A_a.Chr1	2449602	779	+	30427671	TTTCTATTATTTATGTGTTTTTGACGAGCAACTGATTTTCTTTTTTAATAGGAAATGGTACGCTTGATGAAGGAAAAACGTCTTCCAGCCGCATTCAAGTGTTACCATAATTTTCATAGAGCTGAGGACATATCACCTGATAACCTTTTTTATAAATTGGTTGTCCATGTGCACAGCGATTCAGGATTTAGGCGTTACCACAAAGAGATGAGGTTATCTTTAAACACGTCT-TAGTTACTTATGTCTTTCAGTTTGACCAAACTACGAAAGTATTGATTTCTGTTTACTGAACCATATTAATTACCGTGTTTCAGGCACATGCCTAGTTTGTGGCCTTTATATCGAGGTAAATGCTTTTACTTGTTTCTTA--TTCTATAACTTTGCTTTCAAATTCCATGTTAATGCACCATTTCTGATTCATTTGTTAGAGTTCAATAATTTTGTCTCTACTTTGCACAGTTCCTAGTCTTTGTGAAATAACTTTGTTAGGAAATCAAATTTTCTTGGATGTTTCAGATTTGGTTGTTGTGGTAGCATGCTATGGACCATTGCCATCTGTTTTT-C-GTCTAACTATATTAGTTAGTGATCTTCTTTTTTCTTTTCAGGTTTCTTTGTTGATATTAATTTGTTCAAGTCAAACAAAGGGAGGGATCTGATGGCTCTGAAAAGCATTGATA------------ATGCCAGTGAAAATGATGGTCGAGGAGAAAAGGATGGTTTGGCTGATGATGATGCTAACTTAATGATTAAAATGAAGTTTCTCACATACAAGGTAATTTTAG
s	B_b.scaffold_6	87649	793	-	25389057	TTTCCATTATTTATGTGTTTTTGATGAGCAACCGAATTTTTTTTTTAATAGGAAATGGTACGCTTGATGAAGGAAAAACGTCTTCCAGCCGCATTCAAATGTTACCATAATTTTCATAGAGCTGACGACATATCACCTGATAACCTATTTTATAAATTGGTTGTTCATGTGCACAGCGATTTAGGATTTAGGCGTTACCAAAAAGAGATGAGGTTATCTTTAAAGACGTCTCTAGTTACTTATGTCTTTCAGTTTGGTCAAACTAAGAAAATGTTGATTTCTGTATACTGAACCATATTTTTTACCGTATTTCAGGCACATGCCTAGTTTGTGGCCTTTATATCGAGGTAAATATTTTTGTTTGTTTCTTA--TTCTTTAGCTTTGCTCTAAAA-TCCATGATTAGGCACCATTTCTGATTCATTTGCTAGAGTTCAATAGTTTTGTCTCTACTTTCCACAGTTCCTAGTCTTTGTGAAATAACTTTGTTAGGAAATCAAATTTTCTTGGACGTTTCAGATTTGGTTGTTGTGGTAGCATACTATGGATCATTGCCATATGTTTTTNTCGTTTAACTATATTAGTTAGTGATCTTCTTTTTTCTTTTCAGGTTTCTTTGTTGATATTAATTTGTTCAAGTCAAACAAAGGGAGGGATCTGATGGCTCTGAAAAGCATTGATAATGCGGTCAAAGATGCCAGTGAAAATGATGGTCAACGAGAAAAGGATGGTTTGGCTGATGGTGATGCTAACTTAATGATCAAATTAAAGTTTCTCACATACAAAGTAATTTTTT

I look forward to your reply.

[wjwei-handsome]

Hi @rejo27 !

Thank you for your use!

It took me a while to find the problem. The reason is that pairs of - symbols appear in the MAF file you provided, such as:

TTA--TTC
TTA--TTC

Actually, such a pair - symbol has no meaning.

But I think this may not be your own reason. Please forgive me for attributing it to the alignment software, so I adjusted the code for this problem and tested it, the latest version of the code should resolve this error!

Please try it, if you have any further questions, please let me know :)

Best regards,
Wenjie

[rejo27]

Hi, Weijie
The new version solved my problem, thank you.

[jiadong324]

Hi Wenjie,

Thanks for developing this nice tool!

Here is the my MAF file
msa_13048_chr10.txt. I created from mafft output.

I call variants with command wgatools call msa_13048_chr10.maf -n test -l 0 and here are the outputs:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  test
chr10:119953412-119953812       170     .       CCTTCCTCTTCCT   T       .       .       SVTYPE=DEL;SVLEN=12;END=182     GT:QI   1|1:13048_fa_hap1@170@170@P
chr10:119953412-119953812       375     .       TTTTCTTTT       T       .       .       SVTYPE=DEL;SVLEN=8;END=383      GT:QI   1|1:13048_fa_hap1@363@363@P

I found some issues by looking at the MAF file with wgatools tview

1. The SV length is correct, but the position looks strange. The first one is detected at 170, but it should be around 150. The paf conversion gives cigar 170=12D193=8D128=. I am wondering what causes the position difference.
2. It looks like that you only call variants from one query sequence, even though there are other queries. Is that possible to call variants from all queries in the MAF file?

Thanks!

[wjwei-handsome]

Hi @jiadong324 ,

Sorry for my late reply, approaching the Lunar New Year means busy, you know!

Back to the problem itself, I think this is still due to such a case:

AA---CC
AA---CC
AA---CC

But this is unavoidable in msa. So you are right, At present, I only target the first query genome with the command of call. I think it's a great suggestion to expand it to multiple query genomes and provide options.

Maybe I'll finish it during the New Year holidays!

[jiadong324]

Totally, enjoy the New Year holiday!

It's great if you can add the function to call multiple query genomes or a user specified genome. I guess you also use the first genome in MAF as the reference genome and all variants are detected on that reference genome accordingly.

[wjwei-handsome]

Hi @jiadong324 Dr.Lin,

I implemented this function! You can use the latest commit and use the parameter example like this:

 wgatools call -l0 msa_13048_chr10.maf --query-name 13048_fa_hap2

Please let me know if this can solve your problem!
Happy Chinese New Year 🎆

Best,
Wenjie