A tool to search a genome of your choice for a restriction endonuclease recognition site and output the co-ordinates of all cut sites.

genome_RE_sites.pl

#!/usr/bin/perl
use warnings;
use strict;
#############################################################
# Name: Genome RE Sites                                     #
# Author: Phil Ewels                                        #
# Version 1.0 – 25/05/2011                                  #
# ————————————————————————————————————————————————————————— #
# Outputs a file with locations of restriction endonuclease #
# sites (not resulting fragments)                           #
# ————————————————————————————————————————————————————————— #
# Genome RE Sites is licensed under a Creative Commons      #
# Attribution-ShareAlike 3.0 Unported License.              #
# Based on a work at www.tallphil.co.uk.                    #
#############################################################

#===============================
#==== CONFIGURATION OPTIONS ====
#===============================

# Set by command line (optional)
my ($output,$re_search) = @ARGV;
if (!defined $output) {
	die "Usage is genome_RE_sites.pl [output file] [search string]\nLeave blank to use defaults\n";
} elsif(!defined $re_search) {
	# Restriction site to use
	$re_search = ‘AAGCTT’; # HindIII as default
	warn "Using file defaults: search string = $re_search, output = $output\n";
} else {
	warn "Using command line variables: search string = $re_search, output = $output\n";
}

# Path to chromosome fasta files. Replace chromosome number with %s
my $fn_base = ‘D:\Genome Sequences\Mouse\chr%s.fa’;
warn "Looking for Genome Sequences in $fn_base\n\n";

# Chromosomes to use. Default is (1..19,’X',’Y') – other options are ‘MT’ etc.
my @chromosomes = (1..19,’X',’Y');

#==================================
#== END OF CONFIGURATION OPTIONS ==
#==================================

open (OUT,’>’,$output) or die $!;

# go through each chromosome
foreach my $chromosome (@chromosomes) {
	my $filename = sprintf($fn_base, $chromosome);
	open (IN,$filename) or die "Can’t read file: $!";
	warn "Starting Chromosome $chromosome ($filename)\n";
	my $sequence = ”;
	$_ = ; # Remove fasta header
	while (my $line = ) {
		chomp ($line);
		$sequence .= uc($line); # Make everything upper case
	}

	my $offset = 0;
	my $pos = index($sequence, $re_search, $offset);
	while ($pos != -1) {
		print OUT $chromosome."\t". # Chromosome Name
		$pos."\t". # Position Start
		( $pos + length($re_search) )."\n"; # Position Finish
		$offset = $pos + length($re_search);
		$pos = index($sequence, $re_search, $offset);
	}

	close IN;
}