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

<p>If you have these sets in BED format, <a href="http://bedops.readthedocs.io/en/latest/index.html" rel="nofollow noreferrer">BEDOPS</a> lets you do set operations on them with various overlap criteria.</p>

<p>Specifically, you can use <a href="http://bedops.readthedocs.io/en/latest/index.html" rel="nofollow noreferrer">BEDOPS <code>bedmap</code></a> to map hg38 regions to hg37 regions. The <code>bedmap</code> tool lets you assign overlap criteria from as little as one base to as much as full, mutual overlap. </p>

<p>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.</p>

<p>For example, use <code>--exact</code> to get exactly-matching regions between hg37 and hg38:</p>

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

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

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

<p>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 <code>--echo-map</code> option and some command-line work to implement logic to deal with that case more stringently.</p>

<p>To get a straight-up count of no overlaps, build the inverse — or mask — of your hg38 regions, using <a href="https://github.com/ENCODE-DCC/kentUtils" rel="nofollow noreferrer">Kent tools</a> <code>fetchChromSizes</code> and BEDOPS <code>bedops --difference</code>:</p>

<pre><code>$ fetchChromSizes hg38 | awk '{print $1"\t0\t"$2}' | sort-bed - &gt; hg38.bounds.bed
$ bedops --difference hg38.bounds.bed hg38.bed &gt; hg38.masked.bed
</code></pre>

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

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

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

<p>The <a href="http://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#overlap-criteria" rel="nofollow noreferrer">documentation for <code>bedmap</code></a> describes these overlap criteria and operations in more detail.</p>

CC BY-SA 3.0
2017-07-17T19:06:47.940Z
2017-07-17T19:37:29.100Z
2017-07-17T19:37:29.100Z
Alex Reynolds
<p>I wrote <a href="http://aegean.readthedocs.io/en/stable/parseval.html" rel="nofollow noreferrer">ParsEval</a> 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.</p>

<p>Caveat: ParsEval doesn't handle BED files, only GFF3.</p>

CC BY-SA 3.0
2017-07-18T09:47:04.510Z
Ah, if ParsEval also handled BED it would surely be the holy grail of annotation comparison. ;-)
Konrad Rudolph
CC BY-SA 3.0
2017-07-18T18:58:28.037Z
@KonradRudolph Ah, if everyone would use BED for visualization (like it was designed) and GFF3 for generic feature annotation (like it was designed), the world would be a better place. ;-)
Daniel Standage
CC BY-SA 3.0
2017-07-18T18:58:53.257Z
@KonradRudolph I joking...partially. But yes, your point is well taken. :-)
Daniel Standage
CC BY-SA 3.0
2017-07-18T19:13:21.800Z
I was mainly just trying to make a cheap pun based on the tool’s name. But [it’s no secret](https://bioinformatics.stackexchange.com/a/586/29) that I find the G[TF]F format horrible.
Konrad Rudolph
CC BY-SA 3.0
2017-07-17T18:09:10.553Z
2017-07-17T18:09:10.553Z
Daniel Standage
<p>I was trying the code from @Alex Reynolds, however, I got retrieved an error: </p>

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

Error: Option SRR800856_ngs_te_mapper_nonredundant.bed does not start with '--'
</code></pre>

<p>But when I tried with just two bed files: </p>

<pre><code>bedmap --count --exact SRR800856_ngs_te_mapper_nonredundant.bed SRR800856_popoolationte_nonredundant.bed 
</code></pre>

<p>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. </p>

CC BY-SA 4.0
2020-04-01T23:44:08.060Z
The `bedmap` application only takes two arguments: a **reference** file, and a **map** file. 

In a certain use case, it can take one file as an argument, which acts as both a **reference** and **map** dataset (the reference is mapped against itself). However, that's probably not relevant here.

Basically, you can't pass a wildcard argument to `bedmap` — for `bedmap` you must specify either 1 or 2 files. 

Please ask a new question to separate what you're trying to do from this parent question (and tag it with `bedops`) and I should be able to find it, take a look and see if I can help.
Alex Reynolds
CC BY-SA 4.0
2020-04-02T05:24:17.863Z
Thank you! I will post a new question.
Silvia VC
CC BY-SA 4.0
2020-04-01T22:29:24.910Z
2020-04-01T22:29:24.910Z
Silvia VC
<p>I am trying to translate (lift over) bed files describing genomic regions from hg37 to hg38. I have tried both UCSC's <a href="http://genome.ucsc.edu/cgi-bin/hgLiftOver" rel="noreferrer">LiftOver</a> tool and <a href="http://crossmap.sourceforge.net/" rel="noreferrer">CrossMap</a> but saw that they give me different results. I therefore need a way of assessing how correct the results are. </p>

<p>I decided to test by downloading a bed file of the RefSeq genes for hg37 and one for hg38 from UCSC. Then, I run <code>liftOver</code> 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.</p>

<p>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:</p>

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

<p>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. <a href="http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=41553" rel="noreferrer">MIR4454</a>). 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. </p>

<bed><awk><liftover>
How can I compare two bed files?
# How can I compare two bed files?

(Asked by: username_0 on 2017-07-17)

I am trying to translate (lift over) bed files describing genomic regions from hg37 to hg38. I have tried both UCSC's [LiftOver](http://genome.ucsc.edu/cgi-bin/hgLiftOver) tool and [CrossMap](http://crossmap.sourceforge.net/) 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](http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=41553)). 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.

**Tags:** bed, awk, liftover

---

## 3 answers

### Answer by username_1 on 2017-07-17. Score: 3

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 by username_2 on 2017-07-17. Score: 2

I wrote [ParsEval](http://aegean.readthedocs.io/en/stable/parseval.html) 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.

> **Comments:**
> **username_3**: Ah, if ParsEval also handled BED it would surely be the holy grail of annotation comparison. ;-)  
> **username_2**: @KonradRudolph Ah, if everyone would use BED for visualization (like it was designed) and GFF3 for generic feature annotation (like it was designed), the world would be a better place. ;-)  
> **username_2**: @KonradRudolph I joking...partially. But yes, your point is well taken. :-)  
> **username_3**: I was mainly just trying to make a cheap pun based on the tool’s name. But [it’s no secret](https://bioinformatics.stackexchange.com/a/586/29) that I find the G[TF]F format horrible.  

---

### Answer by username_4 on 2020-04-01. Score: 0

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.

> **Comments:**
> **username_1**: The `bedmap` application only takes two arguments: a **reference** file, and a **map** file.   In a certain use case, it can take one file as an argument, which acts as both a **reference** and **map** dataset (the reference is mapped against itself). However, that's probably not relevant here.  Basically, you can't pass a wildcard argument to `bedmap` — for `bedmap` you must specify either 1 or 2 files.   Please ask a new question to separate what you're trying to do from this parent question (and tag it with `bedops`) and I should be able to find it, take a look and see if I can help.  
> **username_4**: Thank you! I will post a new question.  

---

