Add "samtools checksum" command. Needs htslib#1850#2122
Conversation
3986466 to
970fef1
Compare
4d72159 to
2ec7849
Compare
|
Hi James, thank you for this PR and functionality! Can you elaborate more on the logic implemented here? From our testing, there are quite a few cases to consider when comparing matched BAM and CRAM:
|
|
Good questions, and I'm up for looking at test data that exhibits some remaining problems. However generally when doing CRAM checking with This corrects many of the erroneous BAM fields put there by buggy aligners, including MAPQ on unmapped data and CIGAR strings on unmapped data and some BAM flags. It's not likely to change TLEN or other less obvious BAM flags though. It depends on your pipeline I guess, but I think locally we use may use fixmate to correct a lot of the aligner bugs while the data is still name collated, which obscures a lot of the TLEN out by one errors or incorrectly assigned BAM flags. Anyway, if you have cases the round tripping isn't working on them I'm all ears and I'm happy to help add extra rules for known round trip problems. For example there's a new "cigarx" sanitizer option that replaces = and X CIGAR ops with M, and another "cigar" sanitizer which collapsed neighbouring identical operations together (eg 2M3M is 5M) and removes 0 length operations. Bizarrely both of those turn up in the wild - again likely due to buggy upstream processing! It's been quite an eye opener realising just how many oddly encoded files there are. :) |
|
I thought I dealt with NM and MD already, but apparently not. You can do I can change the RG however shouldn't need to be removed. If RG is present, CRAM codes it as a number to reference the Nth RG record (just as BAM encodes RNAME as the Nth SQ record). If it finds RG >= 0 then on decode it converts back. This means it should always round trip provided the RG headers exist. I guess technically we could have a BAM file with an RG tag and no RG header, but I'd argue it's broken data again. Maybe this is something the sanitizer should fix. I'll explore more, but have you seen such files in the wild? Part of me also thinks if there is known broken inputs in a non-common way (unlike MAPQ/CIGAR on unaligned data which is sadly very common due to bwa doing this), then it's not that unfair to require users to explicitly add extra options to filter out those mistakes. It at least brings to their attention a problem with their input data, which could be seen as a benefit. |
4dfd64d to
50bac78
Compare
|
I think I'm ready for review on this new command. There may be a need to add more sanitizer options or recommend fixmate etc first for CRAM round trips, but if there are differences highlighted then perhaps that's valid and correct. |
dkj
left a comment
There was a problem hiding this comment.
From my testing and comparison with bamseqchksum I think this is exactly what we need. Thank you for the -B option!
23f9965 to
d72977f
Compare
The first option collapses duplicate opcodes into a single one, so 2M3M becomes 5M. This is enabled by default in the "on" and "all" categories. The second option replaces "=" and "X" with "M" (and collapsing as per cigdup). This is not technically an error and so is not enabled in any mode unless explicitly asked for. However it is a useful tool when round-tripping to CRAM and back as it does not support the "=" and "X" opcodes. Also correct for entirely pointless CIGAR entries. Eg 10M0D10M becomes 20M with the 0D removed. These are an outcome of replacing = and X with M, but surprisingly they also turn up in some CIGAR strings that don't use = and X.
This computes and reports checksums on the data within SAM, BAM and CRAM files. There are several modes of operation. - By default it computes checksums on the original sequencing data in an order agnostic fashion: primary sequence, quality, read names and some barcode related tags. This can be used to track whether any of the original FASTQish data has been lost at various stages on the pipelines. This aims to be compatible with biobambam2's bamseqchksum command in functionality, albeit not identical in output format. - With -a it performs full analysis of all alignment data including order and orientation and (almost) all aux tags. This is useful for validating file format conversions. - With -m it can merge the checksum output from multiple files to produce a report that matches the sums we would see if we performed a samtools merge operation. This is useful for validating split and merge operations to ensure we have not lost any records.
- Check a couple more memory allocation errors. - Avoid shadowing variables - Renamed some of the single letter variables to help searching - Avoid a memory leak by (inappropriately) mixing the -t and -m options
This serves two purposes.
bamseqchksum
It's a replacement for biobambam's bamseqchksum tool (but about twice as fast), which is heavily used at Sanger. It should give compatible results with the
-qvoption.This checksums sequence, name, quality and barcode tags in an order and orientation agnostic manner, to facilitate validation of no data loss between raw fastq (or unmapped crams) through alignment, markdup, sorting, and other processing operations to get to the final aligned bam/cram.
vs
If we want the QC pass/fail, we can use
-qfor that:Unlike bamseqchksum we report the fail as well as pass, but only if non-zero. (Unless
-vis enabled where zero rows are reported, for easier parsing)File format round trips
We also have a need to validate that conversion between various file types such as SAM to BAM to CRAM doesn't change anything, or at least nothing that we don't already know about (such as CIGAR =/X vs M).
This is the
-aoption, which is a catch-all for all the other necessary options to force secondary/supplementary, more columns, full BAM flags, all tags, file order and orientation validation, etc.TODO
I have a few more things to work on, so labelling it as draft for now, but this is to get feedback on more ideas.
The order of read-groups should probably be alphanumeric sorted rather than "as it comes" out of the hash table.We could consider ways to combine checksum output files. So file1.chk and file2.chk could be combined together to give the same output that we would expect if we were tosamtools mergeorsamtools catthe two BAMs together and checksum the combined data. This is a simple case of calling thesums_updatefunction to merge CRCs again.Test data