how to compute all restriction enzyme sites in the human genome?

[lunacab]

Dear colleagues,
I have a very simple question to ask but I am struggling with it...
I have a restriction enzyme of 6 nucleotides and i want to find ALL sites within the human genome (hg19 for instance) where the restriction enzyme matches the sequence.
I was trying to use blast but it seems that I am using a too short sequence so it never returns a list.
Any recommendations on how to compute that?
thanks a lot in advance

[malachig]

This should be fun. A real classic bioinformatics task for beginners.

There are some good books out there for learning how to solve these problems.
Beginning Perl for Bioinformatics
Bioinformatics Programming in Python: A Practical Course for Beginners

For working environments you could try:
DNA Linux

This kind of task is also an excellent starting point for learning simple scripting tasks on your own. In other words, you could use this as an excuse to learn some Python, Perl, Regex, Awk, etc.

There are also packages/libraries of code that will have already solved many of these types of basic bioinformatics tasks. To name just a few of these: BioPerl, BioPython, EMBOSS, etc.

[malachig]

In case you feel that my previous post was dodging your question ... attached is an example Perl script that you could use as a starting point. It uses regex to identify occurrences of one string (an RE sequence) within another string (a chromosome).

In this example if you want to get all the EcoRI sites on chromosome 22 you would do this (from a linux prompt):
./findRestrictionSites.pl --genome_version=hg19 --chromosome=22 --re_site=GAATTC

The output will be one site per line in the format: chr:start-end

There is also a list of online RE analysis tools here.

[obig]

If you prefer to use R/Bioconductor, you might investigate the BSgenome and Biostrings packages. Here's a document walking your through the process:

http://www.bioconductor.org/packages/2.3/bioc/vignettes/BSgenome/inst/doc/GenomeSearching.pdf

[lunacab]

Thanks a lot! Very very useful!

[ParthavJailwala]

I have used BioStrings and BSgenome to find restriction sites in the mouse genome...it works great. The only caveat is that you have to use 'matchPattern()' on a per chromosome basis, and then append all the output files if a single per genome file is desired.

[Vandelnokk]

HiCUP

Hi,

check out HiCUP digester in its pipeline:

http://www.bioinformatics.babraham.ac.uk/projects/hicup/scripts_description/#Digester

Best

[craigdj]

Hi lunacab,

Would you be willing to share your data regarding the restriction site coordinates in the human genome? It would be incredibly helpful!

[dariober]

It uses regex to identify occurrences of one string (an RE sequence) within another string (a chromosome).

Just as a comment, if I'm not mistaken your scripts reverse-complements the regular expression, which is something that cannot be done. I'd rather reverse complement the reference sequence even if it is more "expensive".

[blancha]

EMBOSS is an old program, but it works remarkably well for this type of task.
Don't be fooled by the dated website.
It is a very efficient program.

http://emboss.sourceforge.net/apps/cvs/emboss/apps/restrict.html

[malachig]

Just as a comment, if I'm not mistaken your scripts reverse-complements the regular expression, which is something that cannot be done. I'd rather reverse complement the reference sequence even if it is more "expensive".

I'm not sure I follow. It doesn't reverse-complement the "regular expression" it reverse complements the restriction enzyme sequence (a string) that is used in the regular expression. We can either (A) search for our string of interest in the reference sequence and its reverse complement, or (B) search for our string of interest and its reverse complement in the reference sequence.

These two approaches should be equivalent. The script uses option B.