Pysam error trying to uplift bam/sam file coordinates between reference and alt genome.

Question

[thomashartwig80]

Hi, I have allele-specific ChIP-seq data from an F1 hybrid of two maize parents. I am mapping the reads against a referenceand a personal genome of the reference with SNPs and INDELs exchanged. However, when I try to uplift the coordinates from the bam mapped against the personal genome back to the reference,I get the error below. The authors did not respond. The error seems to come from pysam trying the exchange the CIGAR.

Hey after some debugging I get this:

[Hi I tried to debug it and I get the following:

[g2gtools debug] Converting A00742:39:HKM32DSXX:4:1662:30481:35806 10 1288218 41H35M

[g2gtools debug] CIGAR CONVERSION : PHASE 6 : Testing length and conversion

[g2gtools debug] CIGAR SEQ LENGTH=76 != SEQ_LEN=35

[g2gtools debug] old cigar != new cigar

[g2gtools debug] CIGAR CONVERSION : 41H35M ==>

**[g2gtools debug] [(4, 41), (0, 35), (4, -41)]**

Traceback (most recent call last):

  File "/netscratch/dep_psl/grp_frommer/Thomas/bin/miniconda3/envs/g2gtools/bin/g2gtools", line 4, in <module>
    __import__('pkg_resources').run_script('g2gtools==0.2.9', 'g2gtools')
  File "/home/thartwig/.local/lib/python2.7/site-packages/pkg_resources/__init__.py", line 666, in run_script
    self.require(requires)[0].run_script(script_name, ns)
  File "/home/thartwig/.local/lib/python2.7/site-packages/pkg_resources/__init__.py", line 1462, in run_script
    exec(code, namespace, namespace)
  File "/netscratch/dep_psl/grp_frommer/Thomas/bin/miniconda3/envs/g2gtools/lib/python2.7/site-packages/g2gtools-0.2.9-py2.7.egg/EGG-INFO/scripts/g2gtools", line 132, in <module>
    G2GToolsApp()
  File "/netscratch/dep_psl/grp_frommer/Thomas/bin/miniconda3/envs/g2gtools/lib/python2.7/site-packages/g2gtools-0.2.9-py2.7.egg/EGG-INFO/scripts/g2gtools", line 99, in __init__
    getattr(self, args.command)()
  File "/netscratch/dep_psl/grp_frommer/Thomas/bin/miniconda3/envs/g2gtools/lib/python2.7/site-packages/g2gtools-0.2.9-py2.7.egg/EGG-INFO/scripts/g2gtools", line 102, in convert
    g2gtools.g2g_commands.command_convert(sys.argv[2:], self.script_name + ' convert')
  File "/netscratch/dep_psl/grp_frommer/Thomas/bin/miniconda3/envs/g2gtools/lib/python2.7/site-packages/g2gtools-0.2.9-py2.7.egg/g2gtools/g2g_commands.py", line 121, in command_convert
    bsam.convert_bam_file(vci_file=args.vci, file_in=args.input, file_out=args.output, reverse=args.reverse)
  File "/netscratch/dep_psl/grp_frommer/Thomas/bin/miniconda3/envs/g2gtools/lib/python2.7/site-packages/g2gtools-0.2.9-py2.7.egg/g2gtools/bsam.py", line 486, in convert_bam_file
    alignment_new.cigar = convert_cigar(alignment.cigar, read_chr, vci_file, alignment.seq, read1_strand, alignment.pos)
  File "pysam/libcalignedsegment.pyx", line 2651, in pysam.libcalignedsegment.AlignedSegment.cigar.__set__
  **File "pysam/libcalignedsegment.pyx", line 2220, in pysam.libcalignedsegment.AlignedSegment.cigartuples.__set__
OverflowError: can't convert negative value to uint32_t**

[Asaf]

Looks like the CIGAR length is the problem, can you verify that in you BAM file?

[thomashartwig80]

Hi Asaf,

yes it looks like an issue with hard clipped secondary mapped reads. For some reason, g2gtools set their length as a negative value. When we deactivate hard clipping in BWA the tool runs through. Unfortunately, there seem to be additional errors with the CIGAR length.

[Asaf]

So your problem is not in pysam, it's the input