ValueError: quality and sequence mismatch: 98 != 97

[mkankain-hruh]

Bamsurgeon downloaded 20/05/2022. Using pysam 0.19.0 py39h5030a8b_0 from bioconda. Below last log lines before crash. Inputing variants to ultradeep BAM file with mean x500

INFO 2022-05-20 02:16:51,517 secondary reads count:0
INFO 2022-05-20 02:16:51,517 supplementary reads count:0
INFO 2022-05-20 02:16:51,517 loaded 12686 reads, (0 excluded, 0 null or secondary or supplementary--> ignored)

Traceback (most recent call last):
File "/usr/local/source/bamsurgeon/bin/addsnv.py", line 479, in
run()
File "/usr/local/source/bamsurgeon/bin/addsnv.py", line 476, in run
main(args)
File "/usr/local/source/bamsurgeon/bin/addsnv.py", line 422, in main
rr.replace_reads(args.bamFileName, outbam_mutsfile, args.outBamFile, keepqual=True, seed=args.seed)
File "/usr/local/source/bamsurgeon/bin/bamsurgeon/replace_reads.py", line 165, in replace_reads
rdict[extqname].qual = read.qual
File "pysam/libcalignedsegment.pyx", line 2772, in pysam.libcalignedsegment.AlignedSegment.qual.set
File "pysam/libcalignedsegment.pyx", line 1514, in pysam.libcalignedsegment.AlignedSegment.query_qualities.set
ValueError: quality and sequence mismatch: 98 != 97

[drmrgd]

I am getting a very similar error too, which I think may be related to an issue matching the quality string and sequence string for mate pairs:

INFO 2022-08-11 21:18:20,977 loading donor reads into dictionary...

INFO 2022-08-11 21:18:21,042 secondary reads count:0
INFO 2022-08-11 21:18:21,042 supplementary reads count:0
INFO 2022-08-11 21:18:21,042 loaded 28448 reads, (0 excluded, 0 null or secondary or supplementary--> ignored)

Traceback (most recent call last):
  File "/Users/jg7758/Dropbox/biofx_tools/bamsurgeon/bin/addindel.py", line 389, in <module>
    run()
  File "/Users/jg7758/Dropbox/biofx_tools/bamsurgeon/bin/addindel.py", line 386, in run
    main(args)
  File "/Users/jg7758/Dropbox/biofx_tools/bamsurgeon/bin/addindel.py", line 334, in main
    rr.replace_reads(args.bamFileName, outbam_mutsfile, args.outBamFile, keepqual=True, seed=args.seed)
  File "/Users/jg7758/OneDrive - Exact Sciences/biofx_tools/bamsurgeon/bin/bamsurgeon/replace_reads.py", line 165, in replace_reads
    rdict[extqname].qual = read.qual
  File "pysam/libcalignedsegment.pyx", line 2772, in pysam.libcalignedsegment.AlignedSegment.qual.__set__
  File "pysam/libcalignedsegment.pyx", line 1514, in pysam.libcalignedsegment.AlignedSegment.query_qualities.__set__
ValueError: quality and sequence mismatch: 90 != 74

When I tried to look into this, I found it interesting that the lengths correspond to each read in a mate pair in a trans way. For example, in my original BAM, I have this read pair:

220211_A00995_0203_BHW2YGDRXY-Crick:2:2101:4182:26083	163	chr10	89692878	60	90M	=	89692894	90	CAATTCACTGTAAAGCTGGAAAGGGACGAACTGGTGTAATGATATGTGCATATTTATTACATCGGGGCAAATTTTTAAAGGCACAAGAGG	@DFEEHDDCEAAFEF?7>DEG+GFEBE9DECIGDAGB?EEGDCAF>@EFCE?DEFACFABC9F8DEDEC9)9FFFE)77G:FCD8EECEF	ZA:i:0	MB:Z:TCTTGATCGTACAC	ZB:i:11	MC:Z:74M	ZC:f:0.108911	MD:Z:90	PG:Z:bwa	NM:i:0	MQ:i:60OQ:Z:FFFFFFFF:FFFFFF:,:FFF,FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFF:FFFFFFF:,:FFFF,::F:FFF:FFFFFFFF:F:,:FF:	UQ:i:0	AS:i:90	CS:Z:GCCGTCGTTTTAT	XS:i:75	QX:Z:FFFFFFFFFFFFFF-FFFFFFFFFFFFF	RX:Z:TCTTGATCGTACAC-GCCGTCGTTTTAT
220211_A00995_0203_BHW2YGDRXY-Crick:2:2101:4182:26083	83	chr10	89692894	60	74M	=	89692878	-90	TGGAAAGGGACGAACTGGTGTAATGATATGTGCATATTTATTACATCGGGGCAAATTTTTAAAGGCACAAGAGG	CDGFGHCDFB;GFBGDD:EDAHDEGEAEDDEGH:BCFFAEEBBHCE2DDD@H:G;EFGE:HH7CGHAHF8EHB9	MB:Z:TCTTGATCGTACAC	MC:Z:90M11S	MD:Z:74	PG:Z:bwa   	NM:i:0	MQ:i:60	OQ:Z:FFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFF:F:F:FFFF:FF,FFFFFF,FFF:	UQ:i:0	AS:i:74	CS:Z:GCCGTCGTTTTAT	XS:i:59	QX:Z:FFFFFFFFFFFFFF-FFFFFFFFFFFFF	RX:Z:TCTTGATCGTACAC-GCCGTCGTTTTAT

The first is 90 bases long and the second is 74, and the sequence length and quality length are correct. In the temporary mutation BAM that is being compared to the original BAM, I can see that the sequence strings look correct (i.e. the 5bp deletion I requested is in the correct location), the CIGAR is updated correctly, and the string lengths look correct:

samtools view addindel.48284eaa-12e3-47fd-b813-277a15003e8d.muts.bam | grep 220211_A00995_0203_BHW2YGDRXY-Crick:2:2101:4182:26083
220211_A00995_0203_BHW2YGDRXY-Crick:2:2101:4182:26083	163	chr10	89692878	60	48M5D42M	=	89692894	95	CAATTCACTGTAAAGCTGGAAAGGGACGAACTGGTGTAATGATATGTGTTATTACATCGGGGCAAATTTTTAAAGGCACAAGAGGCCCTA	@DFEEHDDCEAAFEF?7>DEG+GFEBE9DECIGDAGB?EEGDCAF>@EFCE?DEFACFABC9F8DEDEC9)9FFFE)77G:FCD8EECEF	NM:i:5	MD:Z:48^CATAT42	MC:Z:32M5D42M	MQ:i:60	AS:i:79	XS:i:64	XA:Z:chr9,-33676194,42M5D48M,8;BS:i:1
220211_A00995_0203_BHW2YGDRXY-Crick:2:2101:4182:26083	83	chr10	89692894	60	32M5D42M	=	89692878	-95	TGGAAAGGGACGAACTGGTGTAATGATATGTGTTATTACATCGGGGCAAATTTTTAAAGGCACAAGAGGCCCTA	CDGFGHCDFB;GFBGDD:EDAHDEGEAEDDEGH:BCFFAEEBBHCE2DDD@H:G;EFGE:HH7CGHAHF8EHB9	NM:i:5	MD:Z:32^CATAT42	MC:Z:48M5D42M	MQ:i:60AS:i:63	XS:i:48	BS:i:1

So, I think that part is correct.

I dug into replace_reads.py, which is where the error is occurring I think, and it seems like it's choking on the very first read. With some debugging cod in there, I can see that the lengths within this record are correct:

~/Dropbox/biofx_tools/bamsurgeon/bin/bamsurgeon/replace_reads.py -b nosc.bam -r addindel.48284eaa-12e3-47fd-b813-277a15003e8d.muts.bam -o sim.bam --keepqual --progress
[E::idx_find_and_load] Could not retrieve index file for 'nosc.bam'
INFO 2022-08-12 08:57:39,936 loading donor reads into dictionary...

INFO 2022-08-12 08:57:40,017 secondary reads count:0
INFO 2022-08-12 08:57:40,017 supplementary reads count:0
INFO 2022-08-12 08:57:40,018 loaded 28448 reads, (0 excluded, 0 null or secondary or supplementary--> ignored)

new Reads -> []

Curr Read:
220211_A00995_0203_BHW2YGDRXY-Crick:2:2101:4182:26083	163	#10	89692878	60	90M	#10	89692894	90	CAATTCACTGTAAAGCTGGAAAGGGACGAACTGGTGTAATGATATGTGCATATTTATTACATCGGGGCAAATTTTTAAAGGCACAAGAGG	array('B', [31, 35, 37, 36, 36, 39, 35, 35, 34, 36, 32, 32, 37, 36, 37, 30, 22, 29, 35, 36, 38, 10, 38, 37, 36, 33, 36, 24, 35, 36, 34, 40, 38, 35, 32, 38, 33, 30, 36, 36, 38, 35, 34, 32, 37, 29, 31, 36, 37, 34, 36, 30, 35, 36, 37, 32, 34, 37, 32, 33, 34, 24, 37, 23, 35, 36, 35, 36, 34, 24, 8, 24, 37, 37, 37, 36, 8, 22, 22, 38, 25, 37, 34, 35, 23, 36, 36, 34, 36, 37])	[('ZA', 0), ('MB', 'TCTTGATCGTACAC'), ('ZB', 11), ('MC', '74M'), ('ZC', 0.10891088843345642), ('MD', '90'), ('PG', 'bwa'), ('NM', 0), ('MQ', 60), ('OQ', 'FFFFFFFF:FFFFFF:,:FFF,FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFF:FFFFFFF:,:FFFF,::F:FFF:FFFFFFFF:F:,:FF:'), ('UQ', 0), ('AS', 90), ('CS', 'GCCGTCGTTTTAT'), ('XS', 75), ('QX', 'FFFFFFFFFFFFFF-FFFFFFFFFFFFF'), ('RX', 'TCTTGATCGTACAC-GCCGTCGTTTTAT')]

Seq Len: 90
qual len: 90

So, I'm wondering if the issue comes from this record being compared to the mate pair and since the lengths of the two mates are different, this is causing the issue. Does that seem possible or am I on the wrong track? Is this something that you've seen before and have a suggestion on how to fix?

[drmrgd]

One more thing of note: it seems to be related to newer code in this repo. When I downgraded and pulled bamsurgeon==1.3 just now to test, it seems to work. So, is there a parity check that's either causing an issue here, or detecting an issue with my input BAM?

[adamewing]

Hmm, will try to make a mixed-length .bam for testing and report back.

Thanks fot the details and debugging.

[adamewing]

This may have been related to #212, if still trying to use bamsurgeon let me know how you go with the latest updates.

[apraga]

Hi, I'm have the same issue with both 1.4.1 and latest on the master branch. It does not happen with version 1.3, as mentioned above.