<p>I want to find the coordinates of all occurrences of the sequence recognized by a restriction enzyme. I know that using EMBOSS I may do this, but this task seems perfectly fitted for short-read sequence alignment software. However, I didn't find any reference for this.</p>

<p>I used bwa for the task and quickly obtained some results. However, to be on the safe side I will like to ask is someone has done something similar or has some advice, perhaps I am stretching the use of bwa. </p>

<p>I tried the following:</p>

<pre><code>echo -e "&gt;DpnII\bGATC" &gt; DnpII.fa
bwa aln -N -n 0 -o 0 -e 0 -l 4 -k 0 dm3 DpnII.fa &gt; DpnII.sai
bwa samse -n 100000000 dm3 DpnII.sai DpnII.fa &gt; DpnII.sam
</code></pre>

<p>The results seem to match the right sites, also the number of sites (489570 for Drosophila) that I obtain are close to what I expect. </p>

<p>Thanks.</p>

Here is an update for people landing on this question: The bwa method that I proposed does not work properly. I don't remember exactly why it was wrong but I stop using it. I think it was returning regions with NNNs.

My solution was to use Biopython to search for a pattern. This is quite fast and accurate. Here is the code that finds and sorts the results:

https://github.com/maxplanck-ie/HiCExplorer/blob/master/hicexplorer/findRestSite.py

[update]
Based on the comments by [Michael Dondrup][1] I fixed an error in the code.


  [1]: https://www.biostars.org/u/55/