sorting bam files

[frymor]

Hi,

I have run a bowtie followed by a tophat with the unmapped reads.

I than want to merge the two sam files (output).

I can't do it because the bam file from the bowtie run is not sorted.

Somehow I can't sort this file.
I tried everything I know. This is how I run my script

bowtie -a --phred64-quals -m 1 -n 2 -l 22 -q --un wt_trimmed.unmapped -t -p 6 --chunkmbs 256 --max wt_trimmed.maxHits -S /export/ref_genome  -1 /export/wt_s1_R1 -2 /export/wt_s1_R2 wt_bowtie.sam

samtools view -h -o wttotal_bowtie.bam -Sb wttotal_bowtie.sam
samtools sort wttotal_bowtie.bam wttotal_bowtie.sorted
samtools index wttotal_bowtie.sorted.bam

#topHat for unmapped reads from bowtie run
tophat -o tophat_out -p 6 -r 100 --solexa1.3-quals /export/ref_genome wt_trimmed1.unmapped wt_trimmed2.unmapped
cp tophat_out/accepted_hits.bam wt_tophat.bam

samtools sort wt_tophat.bam wt_tophat.sorted
samtools index wt_tophat.sorted.bam

As you can see I am sorting both files in the same way.
BUT somehow the tophat file is sorted, the bowtie file is not!
tophat:

@HD VN:1.0 SO:coordinate
@SQ SN:2L LN:23011544
@SQ SN:2LHet LN:368872
@SQ SN:2R LN:21146708
@SQ SN:2RHet LN:3288761
@SQ SN:3L LN:24543557
@SQ SN:3LHet LN:2555491
@SQ SN:3R LN:27905053
@SQ SN:3RHet LN:2517507
@SQ SN:4 LN:1351857
@SQ SN:U LN:10049037
@SQ SN:Uextra LN:29004656
@SQ SN:X LN:22422827
@SQ SN:XHet LN:204112
@SQ SN:YHet LN:347038
@SQ SN:dmel_mitochondrion_genome LN:19517
@PG ID:TopHat VN:1.3.0 CL:/usr/bin/tophat -o tophat_out -p 6 -r 100 --solexa1.3-quals /export/ref_genome wt_trimmed1.unmapped wt_trimmed2.unmapped
HWI-EAS225_0031_FC:2:102:14716:18758#0 147 2L 9069 255 54M = 9077 62 CAAGCCCAAAAAGCAGTTTGATACCAGCGATTTTGTTATTGAGAGCGTGCAGAA GD<BGGDD<>DDG>DCACAA?B=?D<<BDDIDIGBC>?A9BED<BBABA?<D@G NM:i:0 NH:i:1
HWI-EAS225_0031_FC:2:102:14716:18758#0 99 2L 9077 255 54M = 9069 62 AAAAGCAGTTTGATACCAGCGATTTTGTTATTGAGAGCGTGCAGAATATACCAC EDCFBDDGGEEGCCBGEGGEGGDGCCB9@>EF<ECAE8E:CC><:76?=DDBGB NM:i:0 NH:i:1
HWI-EAS225_0031_FC:2:56:6678:11287#0 137 2L 9855 255 54M * 0 0 GTTTTATTAATAATCCTAAGCTAAATACTCAATTATATACTTTATATGGTCGGA IHIHIIIIIGIIHIIGIGHFIIHHHIIIIHIIIIIIIIIIIEIHIIIHHIBHII NM:i:0 NH:i:1

...

bowtie:

@HD VN:1.0 SO:unsorted
@SQ SN:YHet LN:347038
@SQ SN:dmel_mitochondrion_genome LN:19517
@SQ SN:2L LN:23011544
@SQ SN:X LN:22422827
@SQ SN:3L LN:24543557
@SQ SN:4 LN:1351857
@SQ SN:2R LN:21146708
@SQ SN:3R LN:27905053
@SQ SN:Uextra LN:29004656
@SQ SN:2RHet LN:3288761
@SQ SN:2LHet LN:368872
@SQ SN:3LHet LN:2555491
@SQ SN:3RHet LN:2517507
@SQ SN:U LN:10049037
@SQ SN:XHet LN:204112
@PG ID:Bowtie VN:0.12.7 CL:"bowtie -a --phred64-quals -m 1 -n 2 -l 22 -q --un wt_trimmed.unmapped -t -p 6 --chunkmbs 256 --max wt_trimmed.maxHits -S /export/ref_genome -1 /export/wt_s1_R1 -2 /export/wt_s1_R2 wt_bowtie.sam"
HWI-EAS225_0031_FC:2:1:10172:6768#0 163 dmel_mitochondrion_genome 1832 255 54M = 1850 72 TTGGAACAGGATGAACTGTTTATCCACCTTTATCCGCTGGAATTGCTCATGGTG BIIGIIIIIIIIIIIIIEIIIDIIGIEII<GGGGHHIIII>IIIGIIIGGIIEH XA:i:1 MD:Z:0C53 NM:i:1
HWI-EAS225_0031_FC:2:1:10172:6768#0 83 dmel_mitochondrion_genome 1850 255 54M = 1832 -72 TTTATCCACCTTTATCCGCTGGAATTGCTCATGGTGGAGCTTCAGTTGATTTAG HIIFIGGCDGI8IGHHFFIIIIIIIHIIIIIHIIIIIDIIIIHIIIIIIIIIII XA:i:0 MD:Z:54 NM:i:0
HWI-EAS225_0031_FC:2:4:13364:12037#0 163 2L 6805 255 54M = 6819 68 AGAGGTGAAAATATATTAAAATTGCCGCTCATTTTCTTCGCGCTAGAATTAGGA HIGIIDIIIIIIIFIIIIIIIIEIIIIHIIFIIIIIIIIIHH+HGGFFGIIIII XA:i:0 MD:Z:42G11 NM:i:1
...

I just can't figure out what I do wrong.
Does samtools sort has some kind of a threshold for file size, number of lines or anything?

I would appreciate any elp or ideas as to how to try and fix it.

Thanks

Assa

[gaffa]

Does the sorting command run to completion, or do you get any kind of error message? Samtools actually does not modify the SO tag after sorting, so you shouldn't use this tag to determine if the file is sorted. If samtools index runs without complains then the file is probably sorted; you could also just look at the first few reads to see if they are placed in the right order.

[gringer]

You can see the first few reads in what has been posted. Tophat reads are 2L-9069,2L-9077,2L-9855; Bowtie reads are dmel_mitochondrion_genome-1832, dmel_mitochondrion_genome-1852, 2L-6805 (am I picking out the right columns?). From those 3 results, they appear sorted, but it would be a good idea to check through a few more of the reads.

[frymor]

Does the sorting command run to completion, or do you get any kind of error message? Samtools actually does not modify the SO tag after sorting, so you shouldn't use this tag to determine if the file is sorted. If samtools index runs without complains then the file is probably sorted; you could also just look at the first few reads to see if they are placed in the right order.

yes it does.

I compared both outputs. They are both sorted, but the other way around.

tophat sorted my Yhet and dmel at the end of the files and begins with the numbered chromosoems (2L, 2R etc.), while bowtie starts with the numbers.

Is there a way to change it or to change the output file?
Is there a way to sort the files according some specific order?

You can see the first few reads in what has been posted. Tophat reads are 2L-9069,2L-9077,2L-9855; Bowtie reads are dmel_mitochondrion_genome-1832, dmel_mitochondrion_genome-1852, 2L-6805 (am I picking out the right columns?). From those 3 results, they appear sorted, but it would be a good idea to check through a few more of the reads.

yes they are but in a different one. If I want to merge them, I need them to be in the same order. But I can't set it.

Is there any other way to do it?

[Richard Finney]

try "samtools reheader" command

samtools is probably not using the name label, rather an index pointing to the order of the lines in the header.

give both files same header (make sure it matches what you have in the sambam files).

[frymor]

try "samtools reheader" command

samtools is probably not using the name label, rather an index pointing to the order of the lines in the header.

give both files same header (make sure it matches what you have in the sambam files).

this would have bin a good idea, if the samtools reheader function will work.

I get gibberish when running the command

samtools view -H dilp_tophat.bam > header.sam
samtools reheader header.sam bowtie_small.bam

what I get looks like that:

                                                                                                                                ?&?^&
K??Y?6ǴͺB?yI?2j[?&<?;??xm3?SA??Z
?x?<{0=??y???P  ?S:??`#M?&?E?0\,?<??@/??n?Q'????O?|??H#?????S8:???BCO[??
? !$??                                                                   ?m?Yv??^s??ܪ?ۧ??c??]{?B^?~???
      ?e?e^??@0???.Î
                   9b? @????@H??? $֊	ؐ??!?$0????<}zg??3?????5??W????뭷n?z||?/??-??޿}?o?5kk_	????S??GG?_xpp|?????K?^r???ӓ?????\????蹷n9Z??8ELE)c?ּʩ???Zg:K?%??1?]ư_gB??L??e"y?<v?3??;=????|??[???;???zꩧ
_                                                                                                                                                                                                        59%?P??q@Lk?Hn?e:???
t??OU
???+F?i???P???j^?VAo??L?]?5)գ???????\ä6\&WN?Jd"×?F?$???B?U<O?Z>?ܷ???Ygl΅???<??R??0??->??M?
                                                                                         x? F???n??D?Ř???????NN???)??Z?<??J*???ڊ,?s???J?(͍???Ӳ?????SJ??*r???????S???????|????HT/?\/׫?J????]??й??:]ɜ?T??B?exj???e?L?g?âS??"P?D?k?3?spVKZ?Z1VVyay.????vj?V???9??%??{?(?
                                                                                                                                                                                                                                                                    p?T&o???@?8???B?f?I?3K???Ѕ*x??q?g2M#|t?W,?X?$N??#?%?7<??J;^?–??\?%q?J?*??K??mmNi?q???C?
?#=?Q???#???SѬ?????VL=֋?C|?'?J??(+?l?;?d^?Ns?????1?IQj:?S???$?a?F(??C
                                                                     o? ޫ֣??>?AE??3e=6?^EhKa?Q?????X???a????߉^?	??6"?n???P??o?????}??dR???CC'?Z?ԕ??v?T???0??j??b87ޙ?_J?1??~?4V P??#???WtP??AE??jR?c???U	?`	?v-?hE?H?? ?h?'???8?W?؋Q??W?>>???5Ї???Rŵ???'O??ӃVy#T?Qڕ?5?o?9F)P?F??????t????]???uA}?V???*??ȄT7?^?*?Cw?ݻ??????wo!?{;ˇ???/?{?]A8?%;ߡ???????9F??o?=O.<-;??r??sIrƫ?????Y?,&??vp?????P?I??&V7מp??KY??΀??#(U??_?U?M?6&Mb?@????#mL?ge??||?????v$??X???y??A?Qe?Z?"?T?8?rW??.??@(?$I???I?U$J?,R	ph?)????7vP??'w???yQ?????XTPR???_d??h,?7??E??HZ?ފ????ǔ?^4?"z?˳ش%?Iڮu????c?W???`ϕ& U?*??ڊm??y?x\?????
m??????wvz6T?]?e??@??V?p#???S                                        $w~|Ԫ???enJ?T??OJ-?ʔ???p??J< ??2?'???g ?8`?y??????PP5,
                             ?O?R?x/^	<??{?4???!{???=?[)??{???w

A lot of it.

What am I doing wrong here?

Thanks
Assa

[gringer]

by default,

samtools reheader

outputs a bam file. Try

samtools reheader header.sam bowtie_small.bam > output.bam

[frymor]

Thanks, now it works

I should have thought about it.

I still need to try and merge the files. I will post a second response as soon as I figure it out.

Assa

[Robby]

But check your result after using the reheader function. I tried that as well, but after using reheader the position of my reads were different. So I received not just a new header but a complete different mapping

[sommeh]

Dear all,

I have tried all the above recomendations in order to sort my accepted_hits.bam output files from tophat. I still get all the errors mentioned above including:

"Error: sort order of reads in BAMs must be the same"

What could be the problem?

Many thanks!

[tboothby]

I am also getting the same error as sommeh. I have tried sorting, indexing, and changing the header as described above. No luck. Any ideas anyone?

[gringer]

please post the exact commands that you ran, and the headers of output files:

samtools view <file> | head # for BAM files
head <file> # for other text files
grep -v '^@' <file > | head # may also be useful for SAM files

[tboothby]

OK, I have two conditions I want to test. I ran them through tophat and then did the following for each:

samtools view -H accepted_hits.bam > header.sam
samtools reheader header.sam accepted_hits.bam > output.bam
samtools sort output.bam sorted_output

I thought this would sort my two files into the same order. Below are the first few lines of each file.


Condition 1

HWI-1KL118:19:C099JACXX:7:2104:6000:72219	129	lcl|comp0_c0_seq1_FPKM_all255	100M	=	9	108	CGGCCGTCTAGGCCGCGGAAGTTACCTGGTTGATCCTGCCAGTAGTCATATGCTTGTCTCAAAGATTAAGCCATGCATGTCTAAGTATAAACTCTTTGTA	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:1
HWI-1KL118:19:C099JACXX:7:2106:9821:36467	129	lcl|comp0_c0_seq1_FPKM_all255	100M	=	44	143	CGGCCGTCTAGGCCGCGGAAGTTACCTGGTTGATCCTGCCAGTAGTCATATGCTTGTCTCAAAGATTAAGCCATGCATGTCTAAGTATAAACTCTTTGTA	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:1
HWI-1KL118:19:C099JACXX:7:2306:10812:41121	73	lcl|comp0_c0_seq1_FPKM_all255	25M1469N75M	*	0	0	CGGCCGTCTAGGCCGCGGAAGTCGCTACACTGATGAATTCAACGAGTTTTCTACCTGGGCCGATAGGCTCGGGTAATCTTTTGAAATTTCATCGTGATGG	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:2	XS:A:-	NH:i:1
HWI-1KL118:19:C099JACXX:7:1206:6846:6176	113	lcl|comp0_c0_seq1_FPKM_all255	100M	=	17	116	CGGCCGTCTAGGCCGCGGAAGTTACCTGGTTGATCCTGCCAGTAGTCATATGCTTGTCTCAAAGATTAAGCCATGCATGTCTAAGTATAAACTCTTTGTA	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:1
HWI-1KL118:19:C099JACXX:7:1307:18878:68878	113	lcl|comp0_c0_seq1_FPKM_all255	100M	=	108	207	CGGCCGTCTAGGCCGCGGAAGTTACCTGGTTGATCCTGCCAGTAGTCATATGCTTGTCTCAAAGATTAAGCCATGCATGTCTAAGTATAAACTCTTTGTA	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:1
HWI-1KL118:19:C099JACXX:7:2302:20038:198362	113	lcl|comp0_c0_seq1_FPKM_all255	100M	=	96	193	GCCGTCTAGGCCGCGGAAGTTACCTGGTTGATCCTGCCAGTAGTCATATGCTTGTCTCAAAGATTAAGCCATGCATGTCTAAGTATAAACTCTTTGTACT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:1
HWI-1KL118:19:C099JACXX:7:2307:13318:135884	89	lcl|comp0_c0_seq1_FPKM_all255	100M	*	0	0	CGTCTAGGCCGCGGAAGTTACCTGGTTGATCCTGCCAGTAGTCATATGCTTGTCTCAAAGATTAAGCCTTTCATGTCTAAGTATAAACTCTTTGTACTGT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:2	NH:i:1

Condition 2

HWI-1KL118:19:C099JACXX:7:1202:18064:79537	385	comp0_c0_seq1	1	3100M	=	65	164	CGGGCGGTGTGTACAAAGCAATGATCCTTCCGCAGGTTCACCTACGGAAACCTTGTTACGACTTCTCCTTCCTCTAAATGATAAGGCTTAATGAACTTCT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:2	CC:Z:comp1_c0_seq1	CP:i:1750	HI:i:0
HWI-1KL118:19:C099JACXX:7:2304:19725:82467	129	comp0_c0_seq1	1	3100M	=	67	166	CGGGCGGTGTGTACAAAGCAATGATCCTTCCGCAGGTTCACCTACGGAAACCTTGTTACGACTTCTCCTTCCTCTAAATGATAAGGCTTAATGAACTTCT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:2	CC:Z:comp1_c0_seq1	CP:i:1750	HI:i:0
HWI-1KL118:19:C099JACXX:7:1101:17072:30058	113	comp0_c0_seq1	1	3100M	=	76	175	CGGGCGGTGTGTACAAAGCAATGATCCTTCCGCAGGTTCACCTACGGAAACCTTGTTACGACTTCTCCTTCCTCTAAATGATAAGGCTTAATGAACTTCT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:2	CC:Z:comp1_c0_seq1	CP:i:1750	HI:i:0
HWI-1KL118:19:C099JACXX:7:1103:4285:164358	89	comp0_c0_seq1	1	3100M	*	0	0	CGGGCGGTGTGTACAAAGCAATGATCCTTCCGCAGGTTCACCTACGGAAACCTTGTTACGACTTCTCCTTCCTCTAAATGATAAGGCTTAATGAACTTCT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:2	CC:Z:comp1_c0_seq1	CP:i:1750	HI:i:0
HWI-1KL118:19:C099JACXX:7:1103:19190:197593	385	comp0_c0_seq1	1	3100M	=	65	164	CGGGCGGTGTGTACAAAGCAATGATCCTTCCGCAGGTTCACCTACGGAAACCTTGTTACGACTTCTCCTTCCTCTAAATGATAAGGCTTAATGAACTTCT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:2	CC:Z:comp1_c0_seq1	CP:i:1750	HI:i:0
HWI-1KL118:19:C099JACXX:7:1104:8932:26899	369	comp0_c0_seq1	1	3100M	=	70	169	CGGGCGGTGTGTACAAAGCAATGATCCTTCCGCAGGTTCACCTACGGAAACCTTGTTACGACTTCTCCTTCCTCTAAATGATAAGGCTTAATGAACTTCT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:2	CC:Z:comp1_c0_seq1	CP:i:1750	HI:i:0
HWI-1KL118:19:C099JACXX:7:1104:13134:90864	89	comp0_c0_seq1	1	3100M	*	0	0	CGGGCGGTGTGTACAAAGCAATGATCCTTCCGCAGGTTCACCTACGGAAACCTTGTTACGACTTCTCCTTCCTCTAAATGATAAGGCTTAATGAACTTCT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:0	NH:i:2	CC:Z:comp1_c0_seq1	CP:i:1750

[gringer]

Both files seem to be sorted correctly from what I've seen (i.e. contig first, then location, then read name). What is the problem you're having?

I suspect there is some oddity in the input data. I notice that those two files have what seems to be different contig names: lcl|comp0_c0_seq1_FPKM_all in the first file (I think, it would be useful if you could wrap your file output in code rather than quote), and comp0_c0_seq1 in the second file.

If you're trying to merge these files, it is possible that samtools is complaining because the two files have different contigs in the headers, but without the exact warning message and command sequence, it's difficult to say.

[tboothby]

I am trying the following command

cuffdiff /Users/wolniaklab/Desktop/mv_test_data/merged_asm/merged.gtf /Users/wolniaklab/Desktop/mv_test_data/3-5hr_2/3-5hr_2_output/sorted_output.bam /Users/wolniaklab/Desktop/mv_test_data/6-8hr_2/6-8hr_2_output/sorted_output.bam

and I get the following:

wolniaklab:3-5_2vs6-8_2 wolniaklab$ cuffdiff /Users/wolniaklab/Desktop/mv_test_data/merged_asm/merged.gtf /Users/wolniaklab/Desktop/mv_test_data/3-5hr_2/3-5hr_2_output/sorted_output.bam /Users/wolniaklab/Desktop/mv_test_data/6-8hr_2/6-8hr_2_output/sorted_output.bam
You are using Cufflinks v1.3.0, which is the most recent release.
Error: sort order of reads in BAMs must be the same

Edit: I find it strange that contig names are different... I am going rerun tophat and see if the this difference in naming persists.