How can I compare two bed files?

Question

I am trying to translate (lift over) bed files describing genomic regions from hg37 to hg38. I have tried both UCSC's LiftOver tool and CrossMap but saw that they give me different results. I therefore need a way of assessing how correct the results are.

I decided to test by downloading a bed file of the RefSeq genes for hg37 and one for hg38 from UCSC. Then, I run liftOver to translate the hg37 file to hg39 coordinates and now I want to compare the results of the liftover to the bed file I downloaded for hg38.

How can I usefully compare these files? I am hoping for a tool that can report how similar the two are. Ideally, I would like to see three things in the output:

1. Matching regions: those that are shown with identical positions in the liftover results and the UCSC bed file.
2. Non-matching regions: those whose coordinates do not match between the two files.
3. (if possible) fuzzy matches: those that are not identical but pretty close (off by a few nucleotides).

I thought it would be easy to do with a little awk script but it gets quite complicated because certain RefSeq genes are mapped to multiple positions (e.g. MIR4454). And the 3rd requirement (fuzzy) is not trivial to do in awk. I can live without that though, if a tool can give me the other two.

Answer 1

From a set operations viewpoint, consider your hg37 regions as a "reference" set and hg38 as a "map" set.

If you have these sets in BED format, BEDOPS lets you do set operations on them with various overlap criteria.

Specifically, you can use BEDOPS bedmap to map hg38 regions to hg37 regions. The bedmap tool lets you assign overlap criteria from as little as one base to as much as full, mutual overlap.

Once you do mapping, you can do operations, like counting the number of regions that meet the overlap criteria between reference and map sets, and calculate percentages.

For example, use --exact to get exactly-matching regions between hg37 and hg38:

$ bedmap --count --exact hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > exact_counts.txt

You can do fractional ("fuzzy") overlaps, where (for example) at least 50% of an hg38 region has to overlap an hg37 region:

$ bedmap --count --fraction-map 0.5 hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > fuzzy_0p5_counts.txt

Getting a percentage out of a "fuzzy" result can be complicated by overlaps between hg38 regions that overlap a reference hg37 region. A more complicated procedure can be implemented with the --echo-map option and some command-line work to implement logic to deal with that case more stringently.

To get a straight-up count of no overlaps, build the inverse — or mask — of your hg38 regions, using Kent tools fetchChromSizes and BEDOPS bedops --difference:

$ fetchChromSizes hg38 | awk '{print $1"\t0\t"$2}' | sort-bed - > hg38.bounds.bed
$ bedops --difference hg38.bounds.bed hg38.bed > hg38.masked.bed

The regions in hg38.masked.bed make up all of the genomic space in hg38, except for your lift-over regions.

Apply --count --fraction-ref 1.0 to count the number of reference (hg37) elements which completely overlap these masked hg38 regions:

$ bedmap --count --fraction-ref 1.0 hg37.bed hg38.masked.bed | awk '{s+=$1}END{print s;}' > no_counts.txt

The documentation for bedmap describes these overlap criteria and operations in more detail.

Answer 2

I wrote ParsEval to handle these types of comparisons. ParsEval reports a variety of similarity statistics at both the nucleotide level and feature (exon) level. It doesn't explicitly support "fuzzy" matches, but it does report scaled values between 0.0 and 1.0 for similarity, so you could get a similar feel by definig a similarity cutoff for filtering.

Caveat: ParsEval doesn't handle BED files, only GFF3.

Answer 3

I was trying the code from @Alex Reynolds, however, I got retrieved an error:

bedmap --count --exact SRR800856_*.bed | awk '{s+=$1}END{print s;}' > exact_counts.txt
May use bedmap --help for more help.

Error: Option SRR800856_ngs_te_mapper_nonredundant.bed does not start with '--'

But when I tried with just two bed files:

bedmap --count --exact SRR800856_ngs_te_mapper_nonredundant.bed SRR800856_popoolationte_nonredundant.bed 

I don't get any error. Could you tell me what I should add. Thank in advance for your time. P.D. As I am new here I couldn't get to comment on the specific thread.