Feature request: only check first column of genome file for -g option of bedtools intersect

[bentyeh]

Note: this could be considered a sub-issue of issue #1116; i.e., if the feature requested in issue #1116 is implemented, then this issue would become obsolete.

---

Issue

For bedtools intersect -sorted operations using non-lexicographically sorted chromosomes, users are instructed to supply a genome file using the -g <genome_file> option.

Currently, bedtools intersect performs at least 2 checks on the genome file:

1. it must have 2 columns (tab-delimited)
2. the second column (representing chromosome length) must be non-negative

For the purpose of bedtools intersect, however, the second column is irrelevant. Would it be possible to remove these 2 checks (and any additional irrelevant checks)?

Minimum working example

Consider the following files:

• a.bed

  chr1	1	10
  

• b.bed

  chr10	1	10
  chr1	5	15
  

• genome1.bed

  chr10
  chr1
  

• genome2.bed

  chr10	0
  chr1	0
  

• genome3.bed

  chr10	1
  chr1	1
  

Let's try running bedtools intersect with the -sorted option:

bedtools intersect -a a.bed -b b.bed -sorted

# or an equivalent one-liner without creating a.bed and b.bed
# bedtools intersect -a <(echo -e "chr1\t1\t10\n") -b <(echo -e "chr10\t1\t10\nchr1\t5\t15\n") -sorted

Output:

ERROR: Database file b.bed contains chromosome chr10, but the query file does not.
       Please rerun with the -g option for a genome file.
       See documentation for details.

Now, we try to provide a 1-column file specifying our custom chromosome ordering:

bedtools intersect -a a.bed -b b.bed -sorted -g genome1.bed

# bedtools intersect -a <(echo -e "chr1\t1\t10\n") -b <(echo -e "chr10\t1\t10\nchr1\t5\t15\n") -sorted -g <(echo -e "chr10\nchr1\n")

Output

Error: The genome file genome1.bed has no valid entries (are you sure it's a 2-column bedtools genome file). Exiting.

Next, we try a 2-column genome file, with chromosome sizes set to 0

bedtools intersect -a a.bed -b b.bed -sorted -g genome2.bed

# bedtools intersect -a <(echo -e "chr1\t1\t10\n") -b <(echo -e "chr10\t1\t10\nchr1\t5\t15\n") -sorted -g <(echo -e "chr10\t0\nchr1\t0\n")

Output

Error: chr10 has length equal to 0. Each chromosome must have non-zero length. Exiting.

Finally, we try a 2-column genome file, with chromosome sizes set to 1

bedtools intersect -a a.bed -b b.bed -sorted -g genome3.bed

# bedtools intersect -a <(echo -e "chr1\t1\t10\n") -b <(echo -e "chr10\t1\t10\nchr1\t5\t15\n") -sorted -g <(echo -e "chr10\t1\nchr1\t1\n")

Output

chr1	5	10

Note that the genome file used chromosome sizes incompatible with the actual BED output, but this is not checked by bedtools: the genome file said chromosome chr1 had length 1, yet the output contains a region from positions 5 to 10. Therefore, even in the current (v.2.31.1) implementation of bedtools, the 2nd column of the genome file is clearly irrelevant to the actual intersect operation.

Real-world motivation

My motivating example was trying to filter out reads in a BAM file that overlapped specific regions. Both my BAM file and my mask BED file were pre-sorted according to an alphanumeric chromosome ordering (e.g., chr1, chr2, ..., chr10, chr11, ...).

Consequently, bedtools intersect -sorted requires that I provide a genome file for non-lexicographically sorted chromosomes. Ideally, given that both inputs were already grouped by chromosome (and coordinate sorted within each chromosome group), there would be no need to provide a genome file (see issue #1116). But if that would be harder to implement, at least the requirement for a 2-column genome file could be eliminated; a simple 1-column text file specifying the chromosome order should suffice.