search

amplab / snap

Scalable Nucleotide Alignment Program -- a fast and accurate read aligner for high-throughput sequencing data
https://www.microsoft.com/en-us/research/project/snap/
Apache License 2.0
287 stars 66 forks source link

sorting doesn't work with SAM output? #144

Closed eboyden closed 2 years ago

eboyden commented 2 years ago

As stated, using 2.0.0 for Linux. Sorting works properly with -so -bam to either a file or stdout, but with -so -sam output is not sorted properly.

Command:

/snap/v2.0.0/snap-aligner paired /snap/index/ -pairedInterleavedFastq - -t 16 -mrl 30 -s 0 500 -fs -R "@RG\tID:id\tLB:lb\tPL:ILLUMINA\tPU:pu\tSM:sm" -so -S d -o -sam -

Output:

/samtools/samtools-1.14/samtools view NA12878.bam | head
SN0796:846:HK7WNBCX3:1:1210:16081:78207 163 chr1    484665  0   123M    =   484688  150 CCCCTCCTGGGTTCAAGTGATTCTCCTGCCTCAGCCTCCCAAGTAGCTGGGATTACAGGTCCCAGCCACCACGCCTAGCTAATTTTTGTATGTTTAGTAGAGACAGGGTTTCATCATGTTGGC IIIIIIIIIIIIIHIIIIIIIIHIIIIIIIHIIHIIIIIIIIIIHEHIIIIHHIHIIIGIIIIIHEHHIIIIIIEHHHHHHEHHIIIIIHHIIIIIIIFHEHHIIIIIIIIIIHHGHIHIII  PG:Z:SNAP   NM:i:5  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:5007
SN0796:846:HK7WNBCX3:2:2216:9435:8098   99  chr1    11794284    70  121M    =   11794305    150 CGCAGCCTGGCCTGCAGCTGGGGTCAGGCCAGGGGCAGGGGATGAACCAGGGTCCCCACTCCAGCATCACTCACTTTGTGACCATTCCGGTTTGGTTCTCCCGAGAGGTAAAGAACAAAGA   DHDCCHHH@1@CG1FEGI@GEEC<<1D11<1<G@CH/E<CC<0DHHFHECDHIIHCHIIII?1<<CFCHEE1<CCHHHH1C1<CGGCHEHECHH=CCGHH0C<@C/<G<FGEFHI@?FHH/   PG:Z:SNAP   NM:i:1  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:4641
SN0796:846:HK7WNBCX3:2:2108:13540:94820 163 chr1    135234  1   121M    =   135359  254 ACTGGAGATCAAGTTCTGCGCCTGAAGAGGCTGCCAAAAGTCAAAAGCGGGGCCTGGGAAGGCCGCCGAGAGCCATGAGCTGGGCTGGGCCGAAAGAGGCCACTGGGAGGCAGGAGGAGCT   IIIHIIHIEHHHIHEHHIGHCHHIIHHIIIIIHHHHHCHHHHFIGHIIIIIHIIIIEHHEHHHHIIIIGICHCFHHHHIGIGEGHHHHIGHICDHHHGHDHHIHHEHIIICHDHBGF=GE.   PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:4723
SN0796:846:HK7WNBCX3:2:1113:12783:21047 163 chr1    1007794 0   126M    =   1007821 151 GCTCTGTCGCCCAGGCTGGAGTGCAGTGGCGGGATCTCGGCTCACTGCAAGCTCCGCCTCCCGGGTTCACGCCATTCTCCTGCCTCAGCCTCCCAAGTAGCTGGGACTACAGGCGCCCGCCACTAC  IIIIIIIIIIIIIIIIIIIIIIIIIIIHIIIIIIHIIIIIIIIIIIHIHIIIFHIIIIFHIIIHIHIIIIIIHGGHHDGHHIIIIIHIIIHIIIGHHIIIIIIIIEHHHIIIIIDHIIHHHCHEHH  PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:4732
SN0796:846:HK7WNBCX3:2:1101:14300:95049 83  chr1    10824   3   89M =   10824   -89 AGGCGTGGCGCAGGCGCAGAGAGGCGCGCCGCGCCGGCGCAGGCGCAGAGACACATGCTACCGCGTCCAGGGGTGGAGGCGTGGCGCAG   @HF/HHE<EFE=HGDCC1<<D0HC</<EDEC/<</CC<C1HGC/GCGG@FHEHCD@IHHCC<CC<1C?HD@FFCEIHHEHIHIIHHCHH   PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:3090
SN0796:846:HK7WNBCX3:2:2215:8198:57341  99  chr1    3217926 70  107M    =   3217926 107 TCCCCCTCTTGCCGGCCCTCCCTCCTCTGTCACCGATGTTGGGGTCTGGGAAGGTAGCCAAATGCAGGAGGACACCGGAGCAATAAATCACCCCTCATCCAACATCT IIIIIIIIIIIIIIIHIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIHIIIIIHIIIIIIIFHIHIIIHHHHHHEEHIIIIHIHHFHHIIIIIHHII PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:4257
SN0796:846:HK7WNBCX3:2:2216:17499:68595 163 chr1    1007794 0   126M    =   1007821 151 GCTCTGTCGCCCAGGCTGGAGTGCAGTGGCGGGATCTCGGCTCACTGCAAGCTCCGCCTCCCGGGTTCACGCCAATCTCCTGCCTCAGCCTCCCAAGTAGCTGGGACTACAGGCGCCCGCCACTAC  C@EHIHGHCHHHFEEHHHCH1FGHH?HEEHHD?D/D1D11<1<1C@GG1<1CCEG0=E=DGHEDD0@@GHH?GH/<C0<<<FHH0<0DGG<FEGHH/<<<D/<CEGH?EH??FECEHH?EHHG..B  PG:Z:SNAP   NM:i:1  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:3952
SN0796:846:HK7WNBCX3:1:2103:15237:74721 99  chr1    955157  70  128M    =   955331  296 CCATCTGCCAGCTTCCTTGTCCCCTTCGGACAGCACTGCCTCCATGTGGACATGGCTCTGTTCTATCCGTCCAGCCCCGACCACTACCACCAGCAACAGTGCCCCTGTCTCCTGAGCCGCTGTCGGCG    IIIIIIHIIGHIHEHHHIHIIIIIIEHHHIIIEG1CHHIIIIEIFHHFEFHHIHEHIIHHIHH@GEFHECHHH@GHIIHCEHGHHIIIIIIICHHHIHHIHHHIIHHFGHHII?FEHHHHD<@GHDH/    PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:4808
SN0796:846:HK7WNBCX3:2:2215:21200:80073 99  chr1    11796187    70  121M    =   11796208    150 GCTCTCCTGGGCCCCTCACCTGGATGGGAAAGATCCCGGGGACGATGGGGCAAGTGATGCCCATGTCGGTGCATGCCTTCACAAAGCGGAAGAATGTGTCAGCCTCAAAGAAAAGCTGCGT   FGCEDHHHCDEDHII1EHIICCFFGFHHGHIDG<FHHDCCHCFCE<<<FHDFH<<<DC1DHIIC<GF1<CEC<<<DF<CFHIIII@GCEHE<FHD@GCHHG<GH0DHHI@GHHH<DDCF:/   PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:4496
SN0796:846:HK7WNBCX3:2:2209:10280:67461 99  chr1    11794284    70  121M    =   11794305    150 CGCAGCCTGGCCTGCAGCTGGGGTCAGGCTAGGGGCAGGGGATGAACCAGGGTCCCCACTCCAGCTTCACTCACTTTGTGACCATTCCGGTTTGGTTCTCCCGAGAGGTAAAGAACAAAGA   IIIIIIIIIFIIIIIIIIIIIIHHHIIIG1CGHIIGIIIIIIIIHIIIIIIIIIIII1<FHIIII1<DGHHIFHHIIIIIIIIIIIIHFHHHIHIIIIIIIIIIIHIGIIHIIIIIIIIIH   PG:Z:SNAP   NM:i:3  RG:Z:id LB:Z:lb PL:Z:ILLUMINA   PU:Z:pu SM:Z:sm QS:i:5064

Not sure if this intentional, in which case the documentation should make this clearer (as it is for indexing), or not. This is possibly related to https://github.com/amplab/snap/issues/136, but the fix for that seems to be OSX-specific? Was this also fixed in the Linux dev branch?

bolosky commented 2 years ago

Yeah, that seems totally broken. And no, there aren't any fixes for it in the dev branch.

I suspect no one has tried it for ages and so we missed that it broke. I'll take a look at it.

--Bill

From: eboyden @.> Sent: Wednesday, January 5, 2022 6:29 PM To: amplab/snap @.> Cc: Subscribed @.***> Subject: [amplab/snap] sorting doesn't work with SAM output? (Issue #144)

As stated, using 2.0.0 for Linux. Sorting works properly with -so -bam to either a file or stdout, but with -so -sam output is not sorted properly.

Command:

/snap/v2.0.0/snap-aligner paired /snap/index/ -pairedInterleavedFastq - -t 16 -mrl 30 -s 0 500 -fs -R @.***\tID:id\tLB:lb\tPL:ILLUMINA\tPU:pu\tSM:sm" -so -S d -o -sam -

Output:

/samtools/samtools-1.14/samtools view NA12878.bam | head

SN0796:846:HK7WNBCX3:1:1210:16081:78207 163 chr1 484665 0 123M = 484688 150 CCCCTCCTGGGTTCAAGTGATTCTCCTGCCTCAGCCTCCCAAGTAGCTGGGATTACAGGTCCCAGCCACCACGCCTAGCTAATTTTTGTATGTTTAGTAGAGACAGGGTTTCATCATGTTGGC IIIIIIIIIIIIIHIIIIIIIIHIIIIIIIHIIHIIIIIIIIIIHEHIIIIHHIHIIIGIIIIIHEHHIIIIIIEHHHHHHEHHIIIIIHHIIIIIIIFHEHHIIIIIIIIIIHHGHIHIII PG:Z:SNAP NM:i:5 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:5007

SN0796:846:HK7WNBCX3:2:2216:9435:8098 99 chr1 11794284 70 121M = 11794305 150 CGCAGCCTGGCCTGCAGCTGGGGTCAGGCCAGGGGCAGGGGATGAACCAGGGTCCCCACTCCAGCATCACTCACTTTGTGACCATTCCGGTTTGGTTCTCCCGAGAGGTAAAGAACAAAGA @.@@*.**@*.**@./<G<FGEFHI@?FHH/ PG:Z:SNAP NM:i:1 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:4641

SN0796:846:HK7WNBCX3:2:2108:13540:94820 163 chr1 135234 1 121M = 135359 254 ACTGGAGATCAAGTTCTGCGCCTGAAGAGGCTGCCAAAAGTCAAAAGCGGGGCCTGGGAAGGCCGCCGAGAGCCATGAGCTGGGCTGGGCCGAAAGAGGCCACTGGGAGGCAGGAGGAGCT IIIHIIHIEHHHIHEHHIGHCHHIIHHIIIIIHHHHHCHHHHFIGHIIIIIHIIIIEHHEHHHHIIIIGICHCFHHHHIGIGEGHHHHIGHICDHHHGHDHHIHHEHIIICHDHBGF=GE. PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:4723

SN0796:846:HK7WNBCX3:2:1113:12783:21047 163 chr1 1007794 0 126M = 1007821 151 GCTCTGTCGCCCAGGCTGGAGTGCAGTGGCGGGATCTCGGCTCACTGCAAGCTCCGCCTCCCGGGTTCACGCCATTCTCCTGCCTCAGCCTCCCAAGTAGCTGGGACTACAGGCGCCCGCCACTAC IIIIIIIIIIIIIIIIIIIIIIIIIIIHIIIIIIHIIIIIIIIIIIHIHIIIFHIIIIFHIIIHIHIIIIIIHGGHHDGHHIIIIIHIIIHIIIGHHIIIIIIIIEHHHIIIIIDHIIHHHCHEHH PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:4732

SN0796:846:HK7WNBCX3:2:1101:14300:95049 83 chr1 10824 3 89M = 10824 -89 AGGCGTGGCGCAGGCGCAGAGAGGCGCGCCGCGCCGGCGCAGGCGCAGAGACACATGCTACCGCGTCCAGGGGTGGAGGCGTGGCGCAG @@.@@. PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:3090

SN0796:846:HK7WNBCX3:2:2215:8198:57341 99 chr1 3217926 70 107M = 3217926 107 TCCCCCTCTTGCCGGCCCTCCCTCCTCTGTCACCGATGTTGGGGTCTGGGAAGGTAGCCAAATGCAGGAGGACACCGGAGCAATAAATCACCCCTCATCCAACATCT IIIIIIIIIIIIIIIHIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIHIIIIIHIIIIIIIFHIHIIIHHHHHHEEHIIIIHIHHFHHIIIIIHHII PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:4257

SN0796:846:HK7WNBCX3:2:2216:17499:68595 163 chr1 1007794 0 126M = 1007821 151 GCTCTGTCGCCCAGGCTGGAGTGCAGTGGCGGGATCTCGGCTCACTGCAAGCTCCGCCTCCCGGGTTCACGCCAATCTCCTGCCTCAGCCTCCCAAGTAGCTGGGACTACAGGCGCCCGCCACTAC @.**@.1CCEG0=E=DGHEDD0@@***@***.******@***.***%3c1CCEG0=E=DGHEDD0@@GHH?GH/%3cC0%3c%3c%3cFHH0%3c0DGG%3cFEGHH/%3c%3c%3cD/%3cCEGH?EH??FECEHH?EHHG..B PG:Z:SNAP NM:i:1 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:3952

SN0796:846:HK7WNBCX3:1:2103:15237:74721 99 chr1 955157 70 128M = 955331 296 CCATCTGCCAGCTTCCTTGTCCCCTTCGGACAGCACTGCCTCCATGTGGACATGGCTCTGTTCTATCCGTCCAGCCCCGACCACTACCACCAGCAACAGTGCCCCTGTCTCCTGAGCCGCTGTCGGCG @.@@./ PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:4808

SN0796:846:HK7WNBCX3:2:2215:21200:80073 99 chr1 11796187 70 121M = 11796208 150 GCTCTCCTGGGCCCCTCACCTGGATGGGAAAGATCCCGGGGACGATGGGGCAAGTGATGCCCATGTCGGTGCATGCCTTCACAAAGCGGAAGAATGTGTCAGCCTCAAAGAAAAGCTGCGT @.**@*.**@*.***<DDCF:/ PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:4496

SN0796:846:HK7WNBCX3:2:2209:10280:67461 99 chr1 11794284 70 121M = 11794305 150 CGCAGCCTGGCCTGCAGCTGGGGTCAGGCTAGGGGCAGGGGATGAACCAGGGTCCCCACTCCAGCTTCACTCACTTTGTGACCATTCCGGTTTGGTTCTCCCGAGAGGTAAAGAACAAAGA IIIIIIIIIFIIIIIIIIIIIIHHHIIIG1CGHIIGIIIIIIIIHIIIIIIIIIIII1<FHIIII1<DGHHIFHHIIIIIIIIIIIIHFHHHIHIIIIIIIIIIIHIGIIHIIIIIIIIIH PG:Z:SNAP NM:i:3 RG:Z:id LB:Z:lb PL:Z:ILLUMINA PU:Z:pu SM:Z:sm QS:i:5064

Not sure if this intentional, in which case the documentation should make this clearer (as it is for indexing), or not. This is possibly related to #136https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Famplab%2Fsnap%2Fissues%2F136&data=04%7C01%7Cbolosky%40microsoft.com%7C64f8ba91f13442a367de08d9d0bc3b9b%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637770329141487720%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=jPpKMqIfzSgE7ha033ntDb%2FlJ4kRrwRsZ2Fg9JN1GEg%3D&reserved=0, but the fix for that seems to be OSX-specific? Was this also fixed in the Linux dev branch?

- Reply to this email directly, view it on GitHubhttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Famplab%2Fsnap%2Fissues%2F144&data=04%7C01%7Cbolosky%40microsoft.com%7C64f8ba91f13442a367de08d9d0bc3b9b%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637770329141487720%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=6BQvRykdq9CKQCQe9t2rP2UlfWT7LQmsnjBPzeqIx1Y%3D&reserved=0, or unsubscribehttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAHPTWIJBJ7NX4N5H67IRJLUUT447ANCNFSM5LLJ7ILQ&data=04%7C01%7Cbolosky%40microsoft.com%7C64f8ba91f13442a367de08d9d0bc3b9b%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637770329141487720%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=aKq3w9uYnYthytJf3Ww6Wj2gLDFpuRIJXmT0jvwtE94%3D&reserved=0. Triage notifications on the go with GitHub Mobile for iOShttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fapps.apple.com%2Fapp%2Fapple-store%2Fid1477376905%3Fct%3Dnotification-email%26mt%3D8%26pt%3D524675&data=04%7C01%7Cbolosky%40microsoft.com%7C64f8ba91f13442a367de08d9d0bc3b9b%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637770329141487720%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=eA%2BVVU8Whtq5X255gExzfChvG2GZ2jIs46qcIZ9rBes%3D&reserved=0 or Androidhttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fplay.google.com%2Fstore%2Fapps%2Fdetails%3Fid%3Dcom.github.android%26referrer%3Dutm_campaign%253Dnotification-email%2526utm_medium%253Demail%2526utm_source%253Dgithub&data=04%7C01%7Cbolosky%40microsoft.com%7C64f8ba91f13442a367de08d9d0bc3b9b%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637770329141487720%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=%2B0gQJbd8LAlm8yGd1Fg11DeBkrt04fD81x1qu%2FbnfDQ%3D&reserved=0. You are receiving this because you are subscribed to this thread.Message ID: @.**@.>>

bolosky commented 2 years ago

I put a fix for this in the dev branch. Give it a try and let me know if you think it's fixed before I put it into master.

We were sorting SAM reads by chromosome and uninitialized memory rather than chromosome and pos. This got messed up when we switched around to make the contigs come out in the same order as in the index. We fixed BAM properly but not SAM.

I think it's fine now.

eboyden commented 2 years ago

Hi, thanks for working on this. There's progress, but it still doesn't seem to be quite right. Piping sorted sam output to samtools view/index got about halfway through (based on bam file size upon exit) when it threw this error:

[E::hts_idx_push] Chromosome blocks not continuous
[E::bam_write_idx1] Read 'SN0796:846:HK7WNBCX3:1:1101:10234:13106' with ref_name='chr9', ref_length=138394717, flags=83, pos=130500851 cannot be indexed

Tail of bam output:

/tools/samtools/samtools-1.14/samtools view NA12878.bam | tail -5
SN0796:846:HK7WNBCX3:2:2216:20002:65969 99  chr22   50627680    70  127M    =   50627707    150 CGAGGTCGTCGGCAAAGATCAGCACGATGTTGGGCGGACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAG @FCFH0CDHHHHCEHHHH@@GCHIHHIIIGF1@HHHDCCCD<CCHC<CCEH<C1@FEGHH?C1CGH=HHIHCHE=@CFC<E?/<<E<C0<D?CHHH?D0==@?0D0:E:D<C=:DFH=.8.8..8/: PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:NA12878    PL:Z:ILLUMINA   PU:Z:pu SM:Z:NA12878    QS:i:3578
SN0796:846:HK7WNBCX3:2:2216:20002:65969 147 chr22   50627707    70  123M    =   50627680    -150    TGTTGGGCGTACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGGGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCAG -G@HE7<:7-HHHD</GC<//</D</C?HCC<0<<D@IHGG<<11HD1D<</DH?C/<//<HG@1D1<1F<F<@DHEHE@EG<1EEGCC<<1CGF@CEHGHG<C11EHGHEEGG<1CHHHHHH PG:Z:SNAP   NM:i:2  RG:Z:id LB:Z:NA12878    PL:Z:ILLUMINA   PU:Z:pu SM:Z:NA12878    QS:i:3982
SN0796:846:HK7WNBCX3:2:2215:20009:45230 163 chr22   50627708    70  121M    =   50627729    150 GTTGGGCGGACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCA   IIIIICHHIIIIIDHIIIIIIDHIHIHIIHHIHHIIHIHHIHHIIIGEHHHIIIIICHHDHG1EDHIGIH?@<CG<CGHGIIGDGHIIHH?@EHHGHCEHGEEHH.BHHHIEHH?EEEHIC   PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:NA12878    PL:Z:ILLUMINA   PU:Z:pu SM:Z:NA12878    QS:i:4919
SN0796:846:HK7WNBCX3:2:2215:20009:45230 83  chr22   50627729    70  129M    =   50627708    -150    CAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCAGGCTCTTTCCGATACCGCAGACCCAGGCG   CIIIIIHIIIIHHHHHDEHGHHHHHHHHHHHHHHIIHIIIIIIHHHFHDHIHFIHHHFIHHFHHHGIIIGIIHHHIIHHHHHIIHEHHHHHHGFG?IIIIIIIIIIHEIHED/HHH@@1@GHECHHHHD   PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:NA12878    PL:Z:ILLUMINA   PU:Z:pu SM:Z:NA12878    QS:i:4535
SN0796:846:HK7WNBCX3:1:1101:10234:13106 83  chr9    130500851   70  125M    =   130500826   -150    GCAACTTCACACACATCTTAACAAACAGCTGCCCCAGCCACCCCAGCTCTGCCTGAATTAATTGAACCCAGTGTGTGTTGTTATTGTTAATTTACATTTTTCTTTGTTTTGAATCTGGTTTACAG   HH@FG<1<1EGF@F?HEECG@HHHFEHEIHHIHFHEIGDCDG?EHHCHGHFHIHFCC1HHEEHFD<HH@EGHIHIIHIIIIHGIHEHFCHHHIHIHHHEHHHHHHIHFGHEFGCEIHCIGHHFHH   PG:Z:SNAP   NM:i:0  RG:Z:id LB:Z:NA12878    PL:Z:ILLUMINA   PU:Z:pu SM:Z:NA12878    QS:i:4477

Chromosome counts in order; something seems to go awry at chr10:

/tools/samtools/samtools-1.14/samtools view NA12878.bam | cut -f3 | uniq -c
2267638 chr1
 580208 chr2
 505434 chr3
  57344 chr4
 275218 chr5
1149756 chr6
 817274 chr7
 188890 chr8
 333774 chr9
      2 chr10
     44 chr11
     24 chr12
     10 chr13
     16 chr14
     26 chr15
     28 chr16
     38 chr17
      6 chr18
     12 chr19
      2 chr20
      4 chr21
      4 chr22
      1 chr9

Note that I did not rebuild the index; it was built using v2.0.0. Please advise if this would be a factor, and rebuilding the index would fix the problem.

bolosky commented 2 years ago

Funny. Just as I was reading this the following assert just went off in the sorter (the highlighted part is that the reads are being emitted in sort order):

_ASSERT(b->length <= readBytes && b->contigAndPos >= previous);

Obviously needs more work

From: eboyden @.> Sent: Thursday, January 6, 2022 8:58 PM To: amplab/snap @.> Cc: Bill Bolosky @.>; Comment @.> Subject: Re: [amplab/snap] sorting doesn't work with SAM output? (Issue #144)

Hi, thanks for working on this. There's progress, but it still doesn't seem to be quite right. Piping sorted sam output to samtools view/index got about halfway through (based on bam file size upon exit) when it threw this error:

[E::hts_idx_push] Chromosome blocks not continuous

[E::bam_write_idx1] Read 'SN0796:846:HK7WNBCX3:1:1101:10234:13106' with ref_name='chr9', ref_length=138394717, flags=83, pos=130500851 cannot be indexed

Tail of bam output:

/tools/samtools/samtools-1.14/samtools view NA12878.bam | tail -5

SN0796:846:HK7WNBCX3:2:2216:20002:65969 99 chr22 50627680 70 127M = 50627707 150 CGAGGTCGTCGGCAAAGATCAGCACGATGTTGGGCGGACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAG @FCFH0CDHHHHCEHHHH@@@.**@*.**@*.***<E?/<<E<C0<D?CHHH?D0==@?0D0:E:D<C=:DFH=.8.8..8/: PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:3578

SN0796:846:HK7WNBCX3:2:2216:20002:65969 147 chr22 50627707 70 123M = 50627680 -150 TGTTGGGCGTACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGGGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCAG @.**@*.**@*.**@.@@.***<C11EHGHEEGG<1CHHHHHH PG:Z:SNAP NM:i:2 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:3982

SN0796:846:HK7WNBCX3:2:2215:20009:45230 163 chr22 50627708 70 121M = 50627729 150 GTTGGGCGGACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCA IIIIICHHIIIIIDHIIIIIIDHIHIHIIHHIHHIIHIHHIHHIIIGEHHHIIIIICHHDHG1EDHIGIH?@@.?EEEHIC<mailto:IIIIICHHIIIIIDHIIIIIIDHIHIHIIHHIHHIIHIHHIHHIIIGEHHHIIIIICHHDHG1EDHIGIH?@@.?EEEHIC> PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:4919

SN0796:846:HK7WNBCX3:2:2215:20009:45230 83 chr22 50627729 70 129M = 50627708 -150 CAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCAGGCTCTTTCCGATACCGCAGACCCAGGCG CIIIIIHIIIIHHHHHDEHGHHHHHHHHHHHHHHIIHIIIIIIHHHFHDHIHFIHHHFIHHFHHHGIIIGIIHHHIIHHHHHIIHEHHHHHHGFG?IIIIIIIIIIHEIHED/HHH@@@.*** PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:4535

SN0796:846:HK7WNBCX3:1:1101:10234:13106 83 chr9 130500851 70 125M = 130500826 -150 GCAACTTCACACACATCTTAACAAACAGCTGCCCCAGCCACCCCAGCTCTGCCTGAATTAATTGAACCCAGTGTGTGTTGTTATTGTTAATTTACATTTTTCTTTGTTTTGAATCTGGTTTACAG @.**@*.**@*.**@. PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:4477

- Reply to this email directly, view it on GitHubhttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Famplab%2Fsnap%2Fissues%2F144%23issuecomment-1007141145&data=04%7C01%7Cbolosky%40microsoft.com%7C4b61bacd47fd478948af08d9d19a55e9%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637771283067053659%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=T5kFfRTMvd2KX6PXTFXdcKtyrK2yPyt35wlauAEKKcw%3D&reserved=0, or unsubscribehttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAHPTWMCR7XI4L2SL4QSQMDUUZXG7ANCNFSM5LLJ7ILQ&data=04%7C01%7Cbolosky%40microsoft.com%7C4b61bacd47fd478948af08d9d19a55e9%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637771283067053659%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=tcH0fDtxyrGFM8IXhHX8loY2z5X%2BWh0lxcXc65gDl0I%3D&reserved=0. Triage notifications on the go with GitHub Mobile for iOShttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fapps.apple.com%2Fapp%2Fapple-store%2Fid1477376905%3Fct%3Dnotification-email%26mt%3D8%26pt%3D524675&data=04%7C01%7Cbolosky%40microsoft.com%7C4b61bacd47fd478948af08d9d19a55e9%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637771283067053659%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=EYgcgMKqDktfErOYXPoFIebkGBHBD%2FPNmrqMGQMh200%3D&reserved=0 or Androidhttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fplay.google.com%2Fstore%2Fapps%2Fdetails%3Fid%3Dcom.github.android%26referrer%3Dutm_campaign%253Dnotification-email%2526utm_medium%253Demail%2526utm_source%253Dgithub&data=04%7C01%7Cbolosky%40microsoft.com%7C4b61bacd47fd478948af08d9d19a55e9%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637771283067053659%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=lnZBpmi4NLHlz5MOODGymFl6ojUVAOT%2BRj%2BJFG4NT0A%3D&reserved=0. You are receiving this because you commented.Message ID: @.**@.>>

bolosky commented 2 years ago

OK, try it now.

2.0.1.dev.3.

--Bill

From: eboyden @.> Sent: Thursday, January 6, 2022 8:58 PM To: amplab/snap @.> Cc: Bill Bolosky @.>; Comment @.> Subject: Re: [amplab/snap] sorting doesn't work with SAM output? (Issue #144)

Hi, thanks for working on this. There's progress, but it still doesn't seem to be quite right. Piping sorted sam output to samtools view/index got about halfway through (based on bam file size upon exit) when it threw this error:

[E::hts_idx_push] Chromosome blocks not continuous

[E::bam_write_idx1] Read 'SN0796:846:HK7WNBCX3:1:1101:10234:13106' with ref_name='chr9', ref_length=138394717, flags=83, pos=130500851 cannot be indexed

Tail of bam output:

/tools/samtools/samtools-1.14/samtools view NA12878.bam | tail -5

SN0796:846:HK7WNBCX3:2:2216:20002:65969 99 chr22 50627680 70 127M = 50627707 150 CGAGGTCGTCGGCAAAGATCAGCACGATGTTGGGCGGACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAG @FCFH0CDHHHHCEHHHH@@@.**@*.**@*.***<E?/<<E<C0<D?CHHH?D0==@?0D0:E:D<C=:DFH=.8.8..8/: PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:3578

SN0796:846:HK7WNBCX3:2:2216:20002:65969 147 chr22 50627707 70 123M = 50627680 -150 TGTTGGGCGTACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGGGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCAG @.**@*.**@*.**@.@@.***<C11EHGHEEGG<1CHHHHHH PG:Z:SNAP NM:i:2 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:3982

SN0796:846:HK7WNBCX3:2:2215:20009:45230 163 chr22 50627708 70 121M = 50627729 150 GTTGGGCGGACGGGCAACGGCCAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCA IIIIICHHIIIIIDHIIIIIIDHIHIHIIHHIHHIIHIHHIHHIIIGEHHHIIIIICHHDHG1EDHIGIH?@@.?EEEHIC<mailto:IIIIICHHIIIIIDHIIIIIIDHIHIHIIHHIHHIIHIHHIHHIIIGEHHHIIIIICHHDHG1EDHIGIH?@@.?EEEHIC> PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:4919

SN0796:846:HK7WNBCX3:2:2215:20009:45230 83 chr22 50627729 70 129M = 50627708 -150 CAGGCCAGCAGCCAGGGCCAGGAGGAGGGACCGCGGTGCCCCCATGGACATGGGACCGAGGGGTCTGTCCCAAGAGAGGGAGGGCTACTTGGCTCCAGCAGGCTCTTTCCGATACCGCAGACCCAGGCG CIIIIIHIIIIHHHHHDEHGHHHHHHHHHHHHHHIIHIIIIIIHHHFHDHIHFIHHHFIHHFHHHGIIIGIIHHHIIHHHHHIIHEHHHHHHGFG?IIIIIIIIIIHEIHED/HHH@@@.*** PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:4535

SN0796:846:HK7WNBCX3:1:1101:10234:13106 83 chr9 130500851 70 125M = 130500826 -150 GCAACTTCACACACATCTTAACAAACAGCTGCCCCAGCCACCCCAGCTCTGCCTGAATTAATTGAACCCAGTGTGTGTTGTTATTGTTAATTTACATTTTTCTTTGTTTTGAATCTGGTTTACAG @.**@*.**@*.**@. PG:Z:SNAP NM:i:0 RG:Z:id LB:Z:NA12878 PL:Z:ILLUMINA PU:Z:pu SM:Z:NA12878 QS:i:4477

- Reply to this email directly, view it on GitHubhttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Famplab%2Fsnap%2Fissues%2F144%23issuecomment-1007141145&data=04%7C01%7Cbolosky%40microsoft.com%7C4b61bacd47fd478948af08d9d19a55e9%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637771283067053659%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=T5kFfRTMvd2KX6PXTFXdcKtyrK2yPyt35wlauAEKKcw%3D&reserved=0, or unsubscribehttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAHPTWMCR7XI4L2SL4QSQMDUUZXG7ANCNFSM5LLJ7ILQ&data=04%7C01%7Cbolosky%40microsoft.com%7C4b61bacd47fd478948af08d9d19a55e9%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637771283067053659%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=tcH0fDtxyrGFM8IXhHX8loY2z5X%2BWh0lxcXc65gDl0I%3D&reserved=0. Triage notifications on the go with GitHub Mobile for iOShttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fapps.apple.com%2Fapp%2Fapple-store%2Fid1477376905%3Fct%3Dnotification-email%26mt%3D8%26pt%3D524675&data=04%7C01%7Cbolosky%40microsoft.com%7C4b61bacd47fd478948af08d9d19a55e9%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637771283067053659%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=EYgcgMKqDktfErOYXPoFIebkGBHBD%2FPNmrqMGQMh200%3D&reserved=0 or Androidhttps://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fplay.google.com%2Fstore%2Fapps%2Fdetails%3Fid%3Dcom.github.android%26referrer%3Dutm_campaign%253Dnotification-email%2526utm_medium%253Demail%2526utm_source%253Dgithub&data=04%7C01%7Cbolosky%40microsoft.com%7C4b61bacd47fd478948af08d9d19a55e9%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637771283067053659%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=lnZBpmi4NLHlz5MOODGymFl6ojUVAOT%2BRj%2BJFG4NT0A%3D&reserved=0. You are receiving this because you commented.Message ID: @.**@.>>

eboyden commented 2 years ago

Excellent, works now, thanks.

bolosky commented 2 years ago

Fixed in 2.0.1