Question: Pysam error trying to uplift bam/sam file coordinates between reference and alt genome.
0
gravatar for thomashartwig80
12 months ago by
thomashartwig8010 wrote:

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**
ADD COMMENTlink modified 11 months ago • written 12 months ago by thomashartwig8010

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

ADD REPLYlink written 12 months ago by Asaf8.4k

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.

ADD REPLYlink written 11 months ago by thomashartwig8010

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

ADD REPLYlink written 11 months ago by Asaf8.4k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2012 users visited in the last hour