#!/usr/bin/perl
# perl script NewMapSort.pl
#Uses some BioPerl.
# a test for bio::restriction
#For given enzyme, shows sorted fragment lengths.
#
#
#Outputs to ?excel table?
#
use strict;
use Bio::DB::Fasta;
use Bio::Seq;
use Bio::Restriction::Analysis;
use Bio::PrimarySeq;
use Data::Dumper;
use Bio::Restriction::Enzyme;
my $filename;
my $fragmentstart;
my $sequenceID;
if (defined(@ARGV)){
$filename = @ARGV[0]; } 
else {print "enter filename:";
$filename = <STDIN>;
chomp $filename;}
print "enter fragment start:";
my $fragmentstart =  <STDIN>;
#Open file and see what sequence names are
my $db = Bio::DB::Fasta->new("/home/staffa/clients/colaneria/Mouse_chromosomes/$filename",
				-makeid => \&make_my_id);
my @ids = $db->get_all_ids; 
#print @ids;
print "looking for $fragmentstart\n";
my $ID;
foreach $ID (@ids){
	print "ID= $ID\n";
	if ($ID =~ /W7000001/){
	$sequenceID = $ID;
	print "found $sequenceID\n";
	last;}
	}
#Do Digestion
print "for sequence $sequenceID\n";
my $obj = $db->get_Seq_by_id("$sequenceID");
my $enz = 'HpaII';
my $ra=Bio::Restriction::Analysis->new(-seq=>$obj);
my @frags=$ra->fragments($enz);
print "All sizes, sorted ", join (" ", $ra->sizes($enz, 0, 1)), "\n";
sub make_my_id {
    my $line = shift;
    $line =~ /^>(\w+)\s+/;
    $1;
      }