1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ As a workaround, you could just split your bed file into one file per chromosome, run the intersect for each of those, and then merge the resulting bam files. Clearly not optimal, but if all else fails, it might give you a way out. $\endgroup$ – terdon Sep 25 at 20:36
  • $\begingroup$ Thanks for the suggestion. Unfortunately, that's not a workable option given the workflow that I currently have. $\endgroup$ – KirkD-CO Sep 25 at 20:52
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – KirkD-CO Sep 25 at 21:04
  • 1
    $\begingroup$ samtools sort -n does not sort by reference name (the 3rd column of a bam file). It sorts by read name (ie. the first column of a bam file). $\endgroup$ – Lynn Sep 25 at 21:12
  • $\begingroup$ Inadvertently deleted my first comment. while updating my original post. Oops. Here's what was in my clipboard buffer: Yes, I've done sorting with samtools sort without the -n option, giving the error above. If I sort with the -n option, I get a different error: Error: Sorted input specified, but the file sorted.bam has the following out of order If I look at the header of the .bam file when the -n option wasn't used, the first line is '@HD VN:1.0 SO:coordinate` . I interpret SO:coordinate to mean that the Sort Order is by coordinate, which is consistent with what the documentation says. $\endgroup$ – KirkD-CO Sep 25 at 23:33
1
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Thanks for the response. As I mention above, I have used the -g option that you suggest, which seems the most useful solution for my workflow. The concern I had is that this specified ordering is not consistent with how bedtools would normally sort the .bed file and question is whether this atypical sorting would affect the intersect operation. $\endgroup$ – KirkD-CO Sep 27 at 17:22
  • $\begingroup$ As far as I'm aware, bedtools doesn't care what positional order you sort them in, just that all the files are in the same order. So bedtools sort -g is probably your best option to force a specific chromosomal ordering $\endgroup$ – James Hawley Sep 30 at 13:39
  • $\begingroup$ Thanks, James! That is consistent with what I see, but not knowing what's going on under the hood makes me question things a bit. $\endgroup$ – KirkD-CO Oct 1 at 0:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.