tanghaibao / jcvi

Python library to facilitate genome assembly, annotation, and comparative genomics
BSD 2-Clause "Simplified" License
712 stars 185 forks source link

jcvi.assembly.allmaps path error #590

Closed faye-romero closed 10 months ago

faye-romero commented 10 months ago

Hello @tanghaibao, I am getting the following error when trying to run jcvi.assembly.allmaps path:

#Command: python3 -m jcvi.assembly.allmaps path mybedfile.bed mygenome.fasta -w weights.txt

[12:08:03] DEBUG    Load file `mybedfile.bed`                                                                                                                                                     base.py:34
Traceback (most recent call last):
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/site-packages/jcvi/assembly/allmaps.py", line 2032, in <module>
    main()
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/site-packages/jcvi/assembly/allmaps.py", line 853, in main
    p.dispatch(globals())
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/site-packages/jcvi/apps/base.py", line 136, in dispatch
    globals[action](sys.argv[2:])
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/site-packages/jcvi/assembly/allmaps.py", line 1474, in path
    cc = Map(
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/site-packages/jcvi/assembly/allmaps.py", line 522, in __init__
    bed = Bed(filename)
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/site-packages/jcvi/formats/bed.py", line 156, in __init__
    b = BedLine(line)
  File "/home/fromero3/.conda/envs/allmaps/lib/python3.10/site-packages/jcvi/formats/bed.py", line 61, in __init__
    assert self.start <= self.end, "start={0} end={1}".format(self.start, self.end)
AssertionError: start=1 end=-1

I have used ALLMAPS extensively with success in the past and have only come across this error today. Currently, jcvi.assembly.allmaps merge works just fine, and I have been able to generate the weights file and the input .bed file with no problem. I have also set jcvi and liftover in my path. Any help would be much appreciated trying to debug this error.

tanghaibao commented 10 months ago

@faye-romero

The error suggests that there is a record in the BED file that has start > end.

Would you please check this in your mybedfile.bed? start is the 2nd column and end would be the 3rd column. The faulty line would have 0 in the 2nd column and 0 in the 3rd column, change the end to 1 like below.

-seq    0    0
+seq    0    1
faye-romero commented 10 months ago

@tanghaibao ,

Thank you for the quick reply. I have many lines in mybedfile.bed that look like this:

*   0   -1  mylinkagemap-39:24.154000   *:-1
*   0   -1  mylinkagemap-39:24.179000   *:-1
*   0   -1  mylinkagemap-39:24.179000   *:-1
*   0   -1  mylinkagemap-39:24.179000   *:-1
....

I am assuming these entries are results of markers that were not able to be mapped to my draft genome (i.e. see these matching entries in mylinkagemap.csv):

#Scaffold ID, scaffold position, LG, genetic position
*,-1,39,24.154
*,-1,39,24.179
*,-1,39,24.179
*,-1,39,24.179
...

Should I just create dummy scaffold names for all of these entries, and replace any -1 values with 0? Or just remove the entries entirely?

Please note that I am attempting to use ALLMAPS with several draft genomes generated via different methods. This particular draft was generated using hifiasm trio-binning, so it's a maternal-haplotype-resolved assembly; only 97% of markers were successfully mapped. When I use ALLMAPS on the hifiasm primary assembly, I do not run into this issue, as 99.9% of markers successfully map.

tanghaibao commented 10 months ago

@faye-romero

If these markers don't have scaffold position then they won't contribute to scaffolding. Just remove those for now.

faye-romero commented 10 months ago

@tanghaibao , that did the trick. Thank you for your help.