biod / sambamba

Tools for working with SAM/BAM data
http://thebird.nl/blog/D_Dragon.html
GNU General Public License v2.0
558 stars 104 forks source link

Signed/Unsigned #359

Closed SergejN closed 3 years ago

SergejN commented 6 years ago

Dear sambamba maintainers,

I'm afraid the software handles the coordinates larger than 2^31-1 incorrectly.

original SAM file: D00689:262:CAUR4ANXX:1:1106:1494:1994/1 0 chr2 2551346395 42 88M

sambamba view output: D00689:262:CAUR4ANXX:1:1106:1494:1994/1 0 chr2 -1743620901 42 88M

Is there an easy fix for this issue?

thanks a lot! Sergej

pjotrp commented 6 years ago

Interesting find. @SergejN do you mind pasting a few of those lines fully here? What species is that? I need to check BAM also supports that size. It should.

SergejN commented 6 years ago

Hi Pjotr,

yes, sure. It's the Mexican axolotl (Ambystoma mexicanum). It's genome is 10 times larger than that of a human and split into 14 chromosomes. We are trying to generate chromosome-size scaffolds. Below are a few lines from a SAM file including the header. As far as I know, BAM can't handle sequences longer than 2^31-1 (https://sourceforge.net/p/samtools/mailman/message/28141607/), therefore I was using SAM for output, but maybe there is an easy way to use an unsigned int instead of the signed. Thanks a lot!

@HD VN:1.0 SO:unsorted @SQ SN:chr1 LN:2955898172 @SQ SN:chr2 LN:2923849100 @SQ SN:chr3 LN:2495613285 @SQ SN:chr4 LN:2453368029 @SQ SN:chr5 LN:2630329465 @SQ SN:chr6 LN:3136375085 @SQ SN:chr7 LN:2029504283 @SQ SN:chr8 LN:1710882196 @SQ SN:chr9 LN:1495491496 @SQ SN:chr10 LN:1639492812 @SQ SN:chr11 LN:1436510532 @SQ SN:chr12 LN:1210869585 @SQ SN:chr13 LN:718908455 @SQ SN:chr14 LN:657912494 @PG ID:bowtie2 PN:bowtie2 VN:2.2.9 CL:"bowtie2-align-l --wrapper basic-0 --very-sensitive -x ambMex_genome -p 6 -" 2_D00689:262:CAUR4ANXX:1:1106:1494:1994/1 0 chr2 2551346395 42 88M 0 0 GCATATGAATGCCANTGTCTATGCATAGTACAGTCCAAGCTCCCTGTTAGCACCACTCCCGGGACATGTGTGGAAGATGTGAATGATC BBBBBFFFFFFFFF#BBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:14C73 YT:Z:UU 5_D00689:262:CAUR4ANXX:1:1106:1683:1989/1 16 chr3 811034470 42 125M 0 0 TTGTGTACAAATATTTGCACAGACCGCCCCCCCTACCCCGTCCTTGCTGTCACTATCACTTGGCATCAAATGCACCCTATTCGAGTACACTCCATAATCTCTGCCCCCACNCACCAACCTCCTCT FBB/FFF<FFF/FFB<FFBFFFBFFFFFFBFB<</F<BF//FFFFFFFFFF</<//F/<<<</BFFFB/FB<FF</FF/BFFBFFFFFBFBFFBF<FFBFF<BFFFBB<B#/FFB<FFB/B</BB AS:i:-4 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:75G34G14 YT:Z:UU 7_D00689:262:CAUR4ANXX:1:1106:1622:1994/1 0 chr7 552176181 42 125M 0 0 ATACACAGATTCTTNATTGTGCTAACCAAAAAGAATCAATGTTCAGAAGCTTAGTGTCTATGGGTATCCTAAGCACAAACACATTTTAAATAATGTTTGTTCTGGTAAAAATCAACCCAGATTGA BBBBBFFFFFFFFF#BBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:14C110 YT:Z:UU 9_D00689:262:CAUR4ANXX:1:1106:2059:1993/1 16 chr2 534503978 42 125M 0 0 TATTCCAGGCGAGCATCTGCTTCTGCCTGTAAAAATATGATGCCATGTTGCTTTTCACCAAGTTGCTTTTGTCTAAAACTACACTTAGGAGTGGTGAGATGTTTACTCTANGGTGATAGTGATAT FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBB#FFFFFFFFFBBBBB AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:110G14 YT:Z:UU 12_D00689:262:CAUR4ANXX:1:1106:2380:1992/1 16 chr4 1883888181 1 125M 0 0 AGAATACGATGGTCACTGTTTGCATTTGACTCAATATTGAAGACTTGTATTTTCGAAAATAACCTGTCATCTGCATAAATATAATCAATTGTAGCTTTTAAAGAGCCCTGNTGCCAAGTGTGAGC /BFB/B/F/B<//<FFBFFFFF<FF/</</</FFFFFB/BBBBB<FFB</FF/BFFFFFFB/BB/FB</FFF//F<//FFFBF<FFFFBFFBB<FFBBFFFFFBBFFFB<#FFFB<FFFFBBBBB AS:i:-1 XS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:110T14 YT:Z:UU 11_D00689:262:CAUR4ANXX:1:1106:2442:1985/1 16 chr3 2106236618 42 90M 0 0 GATCTGGTTCCATTTAGTCACAGTCTCGTGGGGGGGTGCTGGTAGTTAGCTAGCGTCCACTTTATGAAATTTGGGNCAAATCCATAGTGG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFB<#FFFFFFFFFBBBBB AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:75C14 YT:Z:UU 1_D00689:262:CAUR4ANXX:1:1106:1247:2000/1 4 0 0 0 0 CAATACACCCCCATNTATGACTGTGCAGTACACTGGGAAAGGAAGACGCNAAACACTTCCTNTCATTGTGCNCATAGGATGAGCTGTTTGGTTGGATCCCCGCCACCTTGGGCTGGAGGAGTCCA BBBBBFFFFFFFFF#<BFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF#<<FFFFFFFFF#<<BFFFFFF#<<FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF YT:Z:UU 8_D00689:262:CAUR4ANXX:1:1106:2103:1985/1 16 chr1 675445417 22 125M 0 0 GGAGATCTCGACCGATATGGAGAACAGCTCACTTGGAGATAATTCTGACAGCCACGAAAATATGGCAAACAATTGACGCTTTAAAATTTCTCGAAAAATTATCTTAAGTTNATAATGTTTTACTG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF/FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFFFFFFFFB<#FFFFFFFFFBBBBB AS:i:-14 XS:i:-64 XN:i:0 XM:i:4 XO:i:0 XG:i:0 NM:i:4 MD:Z:0T1G42C64G14 YT:Z:UU 17_D00689:262:CAUR4ANXX:1:1106:2609:1986/1 0 chr11 807787110 40 8M1I55M1I60M 0 0 CATAGCAATCCAGANCTATAATCTAATATGCAATTTATGTCATTCATTTAATCATGGAATCTCCTTTTTAAACATACAACCTTGTATTATATCTAAAACATACCCTACTAATTAGATTCAAGTGG BBBBBBFFFFFFFF#<<FFFFFFFFFFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFBBB//FFFFFFFF/FFFFFFFFFF<FFFFFFFFFFF/B<F<<FFFFB/BBFBF//F/FBFBFFFF/<B/ AS:i:-17 XN:i:0 XM:i:1 XO:i:2 XG:i:2 NM:i:3 MD:Z:13G109 YT:Z:UU 18_D00689:262:CAUR4ANXX:1:1106:2693:1987/1 0 chr4 1490395106 42 50M1I74M 0 0 CGCTAACCGAAGAANCCGAATTCGATAAGAATCTGTGTCAACATGCCCCTCCCAAACAAGGACATAAAAAGAGAAAGGTACATGACGCCAACAATTGACGATCTGGTCTCGGATCTGCATGGGGC BBBBBFFFFFFFFF#<BFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBB< AS:i:-9 XN:i:0 XM:i:1 XO:i:1 XG:i:1 NM:i:2 MD:Z:14A109 YT:Z:UU 4_D00689:262:CAUR4ANXX:1:1106:1389:1999/1 0 chr6 2031208305 16 38M1I86M 0 0 CCCCAAGTGTATGGNGTGGGTGGAATTTCAAAACACGACCCCACGGCTGTGAGAGGCTGCAGGCACATGCTTTGTCCTGTTGGCGGTCTGAACAGAGCTGTCATCTTTATGATGCGATTTAGGAA BBBBBFFFFFFFFF#BBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:-19 XS:i:-62 XN:i:0 XM:i:3 XO:i:1 XG:i:1 NM:i:4 MD:Z:8T5T81A27 YT:Z:UU 6_D00689:262:CAUR4ANXX:1:1106:1704:1991/1 16 chr8 901520305 6 89M 0 0 GATCAAAGCAAGTTCAAACACAATCAGAAAAAGTATACACACCTCCCTTAATAATAATTTAGACAAGTCTAGGTNAAAACACCACTATT FFBFBFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFBB#FFFFFFFFFBBBBB AS:i:-1 XS:i:-6 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:74C14 YT:Z:UU 20_D00689:262:CAUR4ANXX:1:1106:2841:1992/1 0 chr6 974304881 42 65M1D60M * 0 0 TCACCCCTCCAGATNCTCCTCCCACAGACATGGAACCTCCCCTAGATGGTTTTCCCAAATCCCTCAAAACACCTAAACACACACTGGAGTGCAAACGGTGGGCCGACCCCCCTCCTCAGAATAGA BBBBBFFFFFFFFF#<FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFBFFFFFFFFBFFFFFFFFFFFFFFFFFFFFFF/B AS:i:-14 XN:i:0 XM:i:2 XO:i:1 XG:i:1 NM:i:3 MD:Z:14C50^A45T14 YT:Z:UU

pjotrp commented 6 years ago

That is interesting. We should also check CRAM. If they are also limiting size to 31 bits we ought to raise an issue with samtools. But, yes, there is no reason SAM can not store longer sizes.

SergejN commented 6 years ago

Hi Pjotr, I was just wondering if there is any news or updates.. Thanks!

pjotrp commented 6 years ago

I still need to look at it. Little bit stretched. If you want to have a go fixing that could work too.

pjotrp commented 6 years ago

@SergejN if I can be a co-author I can prioritize. Your project sounds interesting. Even so, it may be a bit tricky depending on what you want to do. E.g. do you want to sort reads? Or is it just view?

SergejN commented 6 years ago

Hi Pjotr,

As of now, I’ve implemented our own “samtools view/sort” replacement in Java. However, it only supports SAM file. This is sufficient to far, although not quite compliant with the SAM format.

BAM is much trickier, since they use int32_t as the type, which must be replaced by LONG and this, in turn, will screw up the whole BAM format. I am not sure I have the time to do that.

Von: Pjotr Prins notifications@github.com Antworten an: biod/sambamba reply@reply.github.com Datum: Mittwoch, 29. August 2018 um 11:52 An: biod/sambamba sambamba@noreply.github.com Cc: SergejN Sergej.Nowoshilow@gmail.com, Mention mention@noreply.github.com Betreff: Re: [biod/sambamba] Signed/Unsigned (#359)

@SergejN if I can be a co-author I can prioritize. Your project sounds interesting. Even so, it may be a bit tricky depending on what you want to do. E.g. do you want to sort reads? Or is it just view?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

pjotrp commented 6 years ago

OK, I hope you understand that it is not a priority at this point. I'll spend some time at possibilities - see if there is an easy fix. D may do better here than samtools' C++.

SergejN commented 6 years ago

I know though that quite some other groups working with large genomes run into similar limitations of the existing tools. Therefore, at some point we want to raise awareness of the limitations of the software in the field.

SergejN commented 6 years ago

Yes, definitely. No problem, sorry if I made an impression I’m pushing. As I said, I was able to solve the issue in our particular case, but there are others that have similar problems. So, I was just curious in the interest of the community.

pjotrp commented 6 years ago

No problem. It is on record now :)

pjotrp commented 3 years ago

No activity