bedtools intersect on very large .bam file - -sorted confusion

Question

I need to do a bedtools intersect operation on a very large .bam file. When I use the standard bedtools intersect operation, the process consumes all the memory on my system.

From the bedtools documentation page, I find the -sorted option to bedtools intersect, which sounds like it is the perfect solution to the problem, but the documentation seems incomplete. From the documentation, I need to ensure that the .bed file is sorted and I can use either bedtools sort or the Linux sort command to accomplish this:

bedtools sort -chrThenSizeA -i input.bed > sorted.bed

sort -k1,1 -k2,2n input.bed > sorted.bed

If I then execute the intersect operation on the .bam file, I get an error.

bedtools intersect -sorted -c -a sorted.bed -b input.bam > intersect_linuxsorted.out
ERROR: chromomsome sort ordering for file input.bam is inconsistent with other files.

Tracking this down I find that samtools sort, which was executed on the .bam file, "sorts by the leftmost coordinate, or the read name if -n is used." It would appear there is no way to put the .bam file in order by reference name, which is what bedtools instructs.

For a workaround, bedtools sort does have the -g option, which allows me to specify the ordering of the .bed file. To get this ordering, I execute:

samtools view input.bam | cut -f3 | unique > order.txt
bedtools sort -g order.txt > alt_sorted.bed

This puts the .bed file into the same order as the .bam file, and the bedtools intersect -sorted command now works, but is it working correctly? It seems that bedtools intersect -sorted wants the files in consistent ordering, which can be done, but is that consistent with the intersect algorithm itself? From the documentation's instructions to sort by the first 2 fields, chromosome name and left-most coordinate, it seems inconsistent to order it as the .bam file is ordered, by read name and not reference name.

I'm not sure how to check this assumption.

EDIT: Some extra details regarding sorting of the .bam file.

If the .bam file is sorted with samtools sort input.bam > sorted.bam I get the error, 'ERROR: chromomsome sort ordering for file input.bam is inconsistent with other files`

If the .bam file is sorted with samtools sort -n input.bam > sorted.bam I get this error, Error: Sorted input specified, but the file sorted.bam has the following out of order

It would seem that if I sort without -n, it thinks the chromosomes are out of order, and if I sort with -n option it thinks the coordinate order is wrong. The only fix I have is forcing the order of the .bed file with the -g option, but again, I have no idea if that is consistent with the underlying algorithm's assumptions for the -sorted option.

Answer 1

It sounds like bedtools is behaving properly. The bam file is sorted by read name only if you used the -n option in samtools sort. Otherwise, if you used samtools sort without the name-sort option, it is 'ordered' by reference name (which can be arbitrarily ordered, as I'll explain below) and then by position (which is numerically ordered).

I think you are getting confused about 'ordering' the reference names. There are multiple ways to sort names; Some people might want numeric sorting (chr1, chr2, chr3...) and others might want lexicographic sorting (chr1, chr10, chr11...). Because of this, many tools allow you to specify the way you want to order the reference names.

There isn't necessary only one correct way to sort the reference names; what matters is that the order is the same between the files you're trying to intersect. This is exactly the purpose of bedtools intersect -g.

If you want to check the order of reference names in your bam, you can easily do so by checking the header via samtools view -h input.bam. It will say something like this:

@HD     VN:1.3  SO:coordinate
@SQ     SN:chr1     LN:248956422
@SQ     SN:chr10    LN:133797422
@SQ     SN:chr11    LN:135086622
@SQ     SN:chr12    LN:133275309

The @SQ lines tell you how your bam is sorted. My suspicion is that this order is different from the Linux sort or bedtools sort default, which is why your bedtools intersect was failing.

Answer 2

As per Lynn's answer, sorting is dependent on how chromosomes/contigs are named, how you want them to be sorted, etc.

samtools sort or sambamba sort will sort as per the SAM specification by the order of the contigs in the header of the file. If you've aligned these sequences to a human reference genome, typically it will be chr1, chr2, ... chr10, chr11, ...

To sort the BED file in the same way, you can specify the sort order with the bedtools sort -g option and providing the explicit order.

Alternatively, use sort with the -V, --version-sort option, that provides a "natural sort of (version) numbers within text".

sambamba sort -o output.bam input.bam
sort -k1,1 -V -k2,2n input.bed > sorted.bed
bedtools intersect -sorted -c -a sorted.bed -b input.bam > intersect_linuxsorted.out