#!/usr/bin/perl -w
use strict;
use SequenceFunctionsDev;
use Getopt::Std;

my %options = ();
getopt("iobmsrdlf", \%options);
#$options{i} = 'D:/primers/probe_sets/skeletal5x.txt';
#$options{o} = 'C:/Documents and Settings/Russ/Desktop/sequence_test';
#$options{b} = 'D:/human_genome_hg16';
#$options{m} = 'rc'; #Regions to mask r = repeat mask c =  conserved
#$options{s} = '10'; #Number of random sets
#$options{r} = 'u'; # regions to search
#$options{d} = '500'; #distance to search
#$options{l} = '20'; #limit size
#$options{f} = 5; #order of model
#my $outdir = 'C:/Documents and Settings/Russ/Desktop/primers/primers/skeletal5xdownstream5kb';
#my $seq2 = 'C:/Documents and Settings/Russ/Desktop/wassermanMasked1000bpUpstream.seq';

if ($options{i}){}else{print"\n\nInput sequence file not defined\n"; err()};
if ($options{o}){}else{print"\n\nOutput directory not defined\n"; err()};
if ($options{b}){}else{print"\n\nBackground directory not defined\n"; err()};
if ($options{r}){}else{print"\n\nSearch region not defined\n"; err()};
if ($options{d}){}else{$options{u} = '1000'};
if ($options{l}){}else{$options{d} = '20'};
if ($options{s}){}else{$options{s} = '0'};
if ($options{f}){}else{$options{f} = '5'};
if ($options{m} =~ /c/){$options{con} = 1}else{$options{con} = 0};
if ($options{m} =~ /r/){$options{rep} = 1}else{$options{rep} = 0};


my $root = undef;
if ($options{i} =~ /\/(\w+)\.\w+$/){$root = $1};

my $genesets = new SequenceFunctionsDev;
$genesets->set_all_directories($options{b}, '', '', '');
$genesets->set_info($options{d}, '0', $options{con}, $options{rep}, $options{l});
#$genesets->load_unique_unigene();
$genesets->load_refgene_full();
$genesets->load_conserved();
$genesets->find_conserved();
$genesets->find_nonconserved();
$genesets->load_gene_set($root, $options{i});
$genesets->find_nonrepeat($root);

if ($options{r} =~ /u/){
    my $region = "$options{d}"."bpUpstream";
    if ($options{m} =~ /c/){$region = "Conserved"."$region"};
    if ($options{m} =~ /r/){$region = "Masked"."$region"};
    my %seq = ();
    $seq{test} = "$options{o}\/$root$region\.seq";
    $seq{svg} = "$options{o}\/$root$region\.svg";
    $seq{background} = "$options{o}\/Background$region\.seq";
    $seq{markov} = "$options{o}\/Background$region\.bfile";
    $seq{scan} = "$options{o}\/Scan$region\.seq";
    $seq{random} = "$options{o}\/$root$region";
    $genesets->list_upstream($root, $seq{test});
    #$genesets->list_upstream_primers($root, '800', '100', '100', "$outdir");
    #$genesets->upstream_sequence_view($root, $seq{svg});
    $genesets->list_genome_upstream($seq{background});
    $genesets->list_scanace_upstream($seq{scan});
    $genesets->calculate_markov($seq{background}, $options{f});
    $genesets->list_markov($seq{markov}, $options{f});
    if ($options{s}){$genesets->make_upstream_randoms($seq{test}, $options{s})};
    print "Upstream regions done\n";
};

if ($options{r} =~ /i/){
    my $region = "$options{d}"."bpIntron";
    if ($options{m} =~ /c/){$region = "Conserved"."$region"};
    if ($options{m} =~ /r/){$region = "Masked"."$region"};
    my %seq = ();
    $seq{test} = "$options{o}\/$root$region\.seq";
    $seq{background}  = "$options{o}\/Background$region\.seq";
    $seq{markov}  = "$options{o}\/Background$region\.bfile";
    $seq{scan} = "$options{o}\/Scan$region\.seq";
    $seq{random} = "$options{o}\/$root$region";
    $genesets->list_first_intron($root, $seq{test});
    #$genesets->list_first_intron_primers($root, '800', '100', '100', "$outdir");
    $genesets->list_genome_first_intron($seq{background});
    $genesets->list_scanace_first_intron($seq{scan});
    $genesets->calculate_markov($seq{background}, $options{f});
    $genesets->list_markov($seq{markov}, $options{f});
    if ($options{s}){$genesets->make_first_intron_randoms($seq{test}, $options{s})};
    print "First Introns done\n";
};

if ($options{r} =~ /d/){
    my $region = "$options{d}"."bpDownstream";
    if ($options{m} =~ /c/){$region = "Conserved"."$region"};
    if ($options{m} =~ /r/){$region = "Masked"."$region"};
    my %seq = ();
    $seq{test} = "$options{o}\/$root$region\.seq";
    $seq{background} = "$options{o}\/Background$region\.seq";
    $seq{markov} = "$options{o}\/Background$region\.bfile";
    $seq{scan} = "$options{o}\/Scan$region\.seq";
    $seq{random} = "$options{o}\/$root$region";
    $genesets->list_downstream($root, $seq{test});
    #$genesets->list_downstream_primers($root, '800', '100', '100', "$outdir");
    $genesets->list_genome_downstream($seq{background});
    $genesets->list_scanace_downstream($seq{scan});
    $genesets->calculate_markov($seq{background}, $options{f});
    $genesets->list_markov($seq{markov}, $options{f});
    if ($options{s}){$genesets->make_downstream_randoms($seq{test}, $options{s})};
    print "Downstream regions done\n";
};

sub err {
    print "Usage\: perl GenomeExtractor.pl \-i <file.seq> <options>\n";
    print"\-i    Text file containing RefSeq accession numbers\n";
    print"\-o    Output directory\n";
    print"\-b    Directory containing background files\n";
    print"\-m    Regions to mask\n";
    print"       r = Repeat\n";
    print"       c = Non\-conserved\n";
    #print"\-s    Number of random sets\n";
    print"\-r    Regions to search\n";
    print"       u = Upstream\n";
    print"       i = Intron\n";
    print"       d = Downstream\n";
    print"\-d    Number of bases to search upstream of the transcriptional start site\n";
    print"\-l    Minimum acceptible sequence length\n";
    die;
};




