Bio::Grep::Backend::BackendI - Superclass for all back-ends

#############################################################################
#   $Author: markus $
#     $Date: 2009-11-12 19:54:03 +0100 (Thu, 12 Nov 2009) $
# $Revision: 1848 $
#############################################################################

package Bio::Grep::Backend::BackendI;

use strict;
use warnings;

use Carp::Assert;
use version; our $VERSION = qv('0.10.6');

use English qw( -no_match_vars );
use autodie qw (open close opendir closedir);

use File::Spec;
use File::Copy;
use File::Temp qw/ tempfile tempdir /;
use File::Basename;

use Scalar::Util qw(reftype);
use Cwd 'abs_path';
use IPC::Open3 'open3';
use Readonly;

use Bio::AlignIO;
use Bio::Factory::EMBOSS;
use Bio::Index::Fasta;
use Bio::Seq;
use Bio::SeqIO;

use Bio::Grep::SearchSettings;
use Bio::Grep::Root;

use UNIVERSAL qw(isa);

use base 'Bio::Grep::Root';

use Class::MethodMaker [
    new    => 'new2',
    scalar => [qw / settings _output_fh _output_fn _current_res_id _tmp_var/],
    array  => [qw / _query_seqs _results/],
    hash   => [qw / features/],
    abstract => [
        qw / search get_sequences get_databases
            generate_database _parse_next_res/
    ],
];

Readonly my $DNA_ALPHABET_SIZE => 5;

sub new {
    my $self = shift->new2;

    # initialize standard settings
    my $settings = Bio::Grep::SearchSettings->new();

    $self->settings($settings);

    # assume back-end binary is in path
    $self->settings->execpath(q{});

    $self->features( %{ $self->_get_all_possible_features() } );
    $self->import;
    return $self;
}

sub _get_all_possible_features {
    my ($self) = @_;
    return {
        MISMATCHES         => 1,
        GUMISMATCHES       => 1,
        EDITDISTANCE       => 1,
        INSERTIONS         => 1,
        DELETIONS          => 1,
        FILTERS            => 1,
        NATIVE_ALIGNMENTS  => 1,
        EVALUE             => 1,
        PERCENT_IDENTITY   => 1,
        PROTEINS           => 1,
        ONLINE             => 1,
        UPSTREAM           => 1,
        DOWNSTREAM         => 1,
        SORT               => 1,
        MAXHITS            => 1,
        COMPLETE           => 1,
        QUERY              => 1,
        QUERY_FILE         => 1,
        QUERY_LENGTH       => 1,
        SHOWDESC           => 1,
        QSPEEDUP           => 1,
        HXDROP             => 1,
        EXDROP             => 1,
        REGEX              => 1,
        REVCOM_DEFAULT     => 1,
        DIRECT_AND_REV_COM => 1,
        NATIVE_D_A_REV_COM => 1,    # backend supports searching for
                                    # direct and revcom matches or
                                    # do we have to emulate it?
    };
}

sub _get_databases {
    my ( $self, $suffix ) = @_;
    my %result = ();
    opendir my $DIR, $self->settings->datapath;
    my @files = grep {/${suffix} \z/xms} readdir $DIR;
    closedir $DIR;
FILE:
    foreach my $file (@files) {
        my $prefix = $file;
        $prefix =~ s/$suffix//xms;
        $result{$prefix} = $prefix;
        $file =~ s/$suffix/\.nfo/xms;
        $file = $self->_cat_path_filename( $self->settings->datapath, $file );
        next FILE if !-e $file;
        open my $FILE, '<', $file;
        my $desc = q{};

        while ( my $line = <$FILE> ) {
            chomp $line;
            $desc .= $line;
        }
        $result{$prefix} = $desc;
        close $FILE;
    }
    return %result;
}

sub _filter_result {
    my ( $self, $res ) = @_;
    if ( $self->settings->filters_isset ) {
        foreach my $filter ( @{ $self->settings->filters } ) {
            $filter->message_reset;

            #$res->_real_query( $self->settings->_real_query );
            $filter->search_result($res);
            if ( !$filter->filter ) {
                if ( $filter->delete ) {
                    return 0;
                }
            }
            if ( $filter->message_isset ) {
                $res->remark( $filter->message );
            }
        }
    }
    return $res;
}

sub results {
    my ($self) = @_;
    $self->deprecated(q{Use 'next_res()' instead.});
    my @results;
    while ( my $res = $self->next_res ) {
        push @results, $res;
    }
    return wantarray ? @results : \@results;
}

sub _prepare_results {
    my ($self) = @_;

    # sorting by dG is a Bio::Grep sortmode,  not a back-end sortmode. so we
    # have to load the results in memory and sort them.
    # next_res() will then shift the next result in _results() (see
    # next_res())
    if ( $self->settings->sort_isset
        && substr( $self->settings->sort, 0, 1 ) eq 'g' )
    {
        my @results;
        while ( my $res = $self->next_res ) {
            push @results, $res;
        }
        @results = $self->_sort_by_dg(@results);
        $self->_results(@results);

        # reset for next_res iterator
        $self->_current_res_id(0);
    }
    else {
        $self->_results_reset;
    }
    return;
}

sub next_res {
    my ($self) = @_;
    my $id = $self->_current_res_id;
    $self->_current_res_id( $id + 1 );
    if ( $self->_results_isset ) {
        if ( $id >= scalar @{ $self->_results } ) {
            return 0;
        }
        else {
            return $self->_results->[$id];
        }
    }
    else {
        return $self->_parse_next_res();
    }
}

sub _sort_by_dg {
    my ( $self, @results ) = @_;
    foreach my $result (@results) {
        if ( !$result->dG_isset ) {
            $self->warn(
                'Not sorting results by dG because some results have no '
                    . 'dG calculated. Use Bio::Grep::RNA.' );
            return @results;
        }
    }
    if ( $self->settings->sort eq 'gd' ) {
        @results = reverse sort { $a->dG <=> $b->dG } @results;
    }
    else {
        @results = sort { $a->dG <=> $b->dG } @results;
    }
    return @results;
}

sub _results_may_have_gaps {
    my ($self) = @_;
    my $s = $self->settings;
    return 1 if ( defined $s->editdistance && $s->editdistance > 0 );
    return 1 if ( defined $s->insertions   && $s->insertions > 0 );
    return 1 if ( defined $s->deletions    && $s->deletions > 0 );
    return 0;
}

# calculates needleman-wunsch global alignment with the EMBOSS
# implementation
sub _get_alignment {
    my ( $self, $seq_a, $seq_b ) = @_;

    # aligments are easy when settings don't allow gaps
    # back-ends must return hit coordinates
    if ( !$self->_results_may_have_gaps
        && length $seq_a->seq == length $seq_b->seq )
    {
        my $alignment = Bio::SimpleAlign->new( -source => 'Bio::Grep' );
        $alignment->add_seq($seq_a);
        $alignment->add_seq($seq_b);
        return $alignment;
    }

    my $factory = Bio::Factory::EMBOSS->new();
    my $prog    = $factory->program('needle');
    my $outfile = $self->_cat_path_filename( $self->settings->tmppath,
        $seq_a->id . '.out' );
    my $gapopen = '5.0';
    if ( !$self->settings->editdistance_isset ) {
        $gapopen = '50.0';
    }

    $prog->run(
        {   -asequence => $seq_a,
            -bsequence => $seq_b,
            -gapopen   => $gapopen,
            -gapextend => '5.0',
            -outfile   => $outfile
        }
    );

    if ( -e $outfile ) {
        my $alignio_fmt = 'emboss';
        my $align_io    = Bio::AlignIO->new(
            -format => $alignio_fmt,
            -file   => $outfile
        );
        unlink $outfile;
        my $alignment = $align_io->next_aln();
        my $s1        = $alignment->get_seq_by_pos(1);
        my $s2        = $alignment->get_seq_by_pos(2);
        $s1->start( $seq_a->start );
        $s1->end( $seq_a->end );
        $alignment = Bio::SimpleAlign->new( -source => 'Bio::Grep' );
        $alignment->add_seq($s1);
        $alignment->add_seq($s2);
        return $alignment;
    }

}

# concatenates a path and a filename platform-independently
sub _cat_path_filename {
    my ( $self, $path, $filename ) = @_;
    if ( $path eq q{} ) {
        return $filename;
    }
    else {
        return File::Spec->catfile( $path, $filename );
    }
}

sub _check_search_settings {
    my ( $self, $arg_ref ) = @_;

    if ( defined $arg_ref ) {
        $self->settings->set_attributes($arg_ref);
    }
    if ( !defined $self->settings->upstream ) {
        $self->settings->upstream(0);
    }
    if ( !defined $self->settings->downstream ) {
        $self->settings->downstream(0);
    }

    $self->_check_int_features();
    $self->_check_sort_mode();
    $self->_check_database();
    $self->_check_for_unsupported_features();

    if (   defined $self->settings->editdistance
        && defined $self->settings->mismatches
        && $self->settings->editdistance > 0
        && $self->settings->mismatches > 0 )
    {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => q{Can't combine editdistance and mismatches.},
            -value => $self->settings->editdistance . ' and '
                . $self->settings->mismatches,
        );
    }
    if ( $ENV{BIOGREPDEBUG} ) {
        carp $self->settings->to_string();
    }

    # reset all filters
    if ( $self->settings->filters_isset ) {
        foreach my $filter ( @{ $self->settings->filters } ) {
            $filter->reset_filter;
        }
    }
    $self->_results_reset;
    $self->_current_res_id(0);
    return;
}

### SOME IMPORTANT SETTINGS TESTS
### called by _check_search_settings()

sub _check_int_features {
    my ($self) = @_;

    my @int_features
        = map { lc $_ } keys %{ $self->_get_all_possible_features() },
        'reverse_complement', 'no_alignments';

    my %skip = (
        FILTERS            => 1,    # checked later
        NATIVE_ALIGNMENTS  => 1,
        NATIVE_D_A_REV_COM => 1,
        PERCENT_IDENTITY   => 1,    # no user setting
        PROTEINS           => 1,    # no user setting
        SORT               => 1,    # checked later
        QUERY_FILE         => 1,    # checked later
        QUERY              => 1,    # checked later
        EVALUE             => 1,    # no user setting
        REVCOM_DEFAULT     => 1,    # no user setting
        REGEX              => 1,    # no user setting
    );

INT_FEATURE:
    for my $option (@int_features) {
        next INT_FEATURE if defined $skip{ uc $option };

        # check if values are integers and untaint them
        $self->settings->$option(
            $self->is_integer( $self->settings->$option, $option ) );
        if ( !defined $self->settings->$option ) {
            my $reset = $option . '_reset';
            $self->settings->$reset;
        }
    }
    return 1;
}

sub _check_for_unsupported_features {
    my ($self) = @_;

    # some warnings if user requests features that are not available
    my %skip = (
        NATIVE_ALIGNMENTS  => 1,
        NATIVE_D_A_REV_COM => 1,
        PERCENT_IDENTITY   => 1,
        REVCOM_DEFAULT     => 1,
        PROTEINS           => 1,
        EVALUE             => 1,
        REGEX              => 1,    # no user setting
    );

FEATURE:
    for my $feature ( keys %{ $self->_get_all_possible_features() } ) {
        next FEATURE if defined $skip{$feature};
        my $value_ok = 0;
        if ( $feature eq 'GUMISMATCHES' ) {
            $value_ok = 1;
        }
        if ( !defined $self->features->{$feature} ) {
            my $lc_f = lc $feature;
            my $lc_i = $lc_f . '_isset';
            if (   $self->settings->$lc_i
                && $self->settings->$lc_f ne $value_ok )
            {
                $self->warn( $lc_f . ' not available in this back-end.' );
            }
        }
    }
    return;
}

sub _check_database {
    my ($self) = @_;

    if ( !$self->settings->datapath_isset ) {
        $self->settings->datapath(q{./});
    }

    # check if database is set and valid
    my $found_database = 0;
    if ( defined $self->settings->database ) {
        $self->settings->database(
            $self->is_word( $self->settings->database ) );
        my %available_dbs = $self->get_databases;
        if ( !defined $available_dbs{ $self->settings->database } ) {
            $self->throw(
                -class => 'Bio::Root::BadParameter',
                -text  => 'Database not found.',
                -value => $self->settings->database,
            );
        }
    }
    else {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => 'Database not defined.',
        );
    }
    return 1;
}

sub _check_sort_mode {
    my ($self) = @_;

    if ( $self->settings->sort_isset ) {
        my %sort_modes = $self->available_sort_modes();
        if ( defined $sort_modes{ $self->settings->sort } ) {
            my ($sort_mode) = $self->settings->sort =~ /(\w+)/xms;
            $self->settings->sort($sort_mode);    #make taint happy
        }
        else {
            $self->throw(
                -class => 'Bio::Root::BadParameter',
                -text  => 'Sort mode not valid.',
                -value => 'sort mode'
            );
        }
    }
    return 1;
}

### END OF SETTINGS TESTS

# copies the specified fasta file in the data directory and cerate a file
# <databasename>.nfo with the description, specified in the optional 2nd
# argument
# changes directory! so please save the oldpath before calling this method
sub _copy_fasta_file_and_create_nfo {
    my ( $self, $args ) = @_;

    # throw exception if filename looks wrong
    $self->is_word( $args->{basefilename}, 'Fasta filename' );

    my $newfile = $self->_cat_path_filename( $self->settings->datapath,
        $args->{basefilename} );

    $args->{filename} = $newfile;

    my %dbs = $self->get_databases;

    if ( defined $dbs{ $args->{basefilename} } ) {
        $self->warn( "Database with that name already exists.\n"
                . 'Skipping database generation.' );
        $args->{skip} = 1;
        return;
    }

    if ( defined $args->{copy} && $args->{copy} ) {
        copy( $args->{file}, $newfile )
            or $self->throw(
            -class => 'Bio::Root::IOException',
            -text  => q{Can't copy } . $args->{file} . " to $newfile",
            -value => $OS_ERROR,
            );
    }
    else {
        my $abs_path = $self->is_path( abs_path( $args->{file} ) );
        symlink $abs_path, $newfile || $self->throw(
            -class => 'Bio::Root::IOException',
            -text  => q{Can't symlink } . $abs_path . " to $newfile",
            -value => $OS_ERROR,
        );
    }
    if ( defined $args->{description} ) {
        open my $NFOFILE, '>', $newfile . '.nfo';
        print ${NFOFILE} $args->{description}
            or $self->_cannot_print("$newfile.nfo");
        close $NFOFILE;
    }
    return;
}

sub _guess_alphabet_of_file {
    my ( $self, $filename ) = @_;
    my $in = Bio::SeqIO->new( -file => $filename );
    return $in->next_seq->alphabet;
}

sub _bioseq_query {
    my ($self) = @_;
    my $query_obj = $self->settings->query;
    my $query;

    if ( !defined $query_obj ) {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => 'Query not defined.',
        );
    }
    my $db_alphabet
        = $self->get_alphabet_of_database( $self->settings->database );

    if ( eval { $query_obj->isa('Bio::Seq') } ) {
        $query = $query_obj->seq;
    }
    else {
        $query = $query_obj;
        if ( $query =~ m{\A \w+ \z}xms ) {
            if ( $db_alphabet eq 'dna' ) {
                $query =~ tr/uU/tT/;
            }
            $query_obj = Bio::Seq->new(
                -id   => '1',
                -desc => 'Query',
                -seq  => $query
            );
        }
        else {
            $query_obj = Bio::Seq->new( -id => '1', -desc => 'Query' );
        }
    }
    return ( $query, $query_obj, $db_alphabet );
}

# prepares the query, for example calculating the reverse complement
# if necessary
# returns the prepared query. settings->query is unchanged!
sub _prepare_query {
    my $self = shift;
    my ( $query, $seq, $db_alphabet ) = $self->_bioseq_query();

    if ( $seq->alphabet ne $db_alphabet ) {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => 'Alphabet of query and database not equal',
            -value => 'Seq: ' . $seq->alphabet . ", DB: $db_alphabet"
        );
    }
    if (   $self->settings->reverse_complement
        || $self->settings->direct_and_rev_com )
    {
        if ( $db_alphabet eq 'dna' ) {
            if ( defined $self->features->{REVCOM_DEFAULT} ) {
                $seq->desc( $seq->desc . ' (reverse complement)' );
            }
            elsif ( defined $self->features->{NATIVE_D_A_REV_COM}
                && $self->settings->direct_and_rev_com )
            {
            }
            else {
                $query = $seq->revcom->seq;
                $seq->desc( $seq->desc . ' (reverse complement)' );
                $seq->seq($query);
            }
        }
        else {
            $self->throw(
                -class => 'Bio::Root::BadParameter',
                -text =>
                    'Reverse complement only available for DNA databases.',
            );
        }
    }
    elsif ( defined $self->features->{REVCOM_DEFAULT}
        && $db_alphabet eq 'dna' )
    {
        $query = $seq->revcom->seq;
        $seq->seq($query);
    }

    # $self->settings->_real_query( uc $query );
    $self->{_query_obj} = $seq;
    return uc $query;    # $self->settings->_real_query();
}

sub generate_database_out_of_fastafile {
    my ( $self, @args ) = @_;
    $self->deprecated( "generate_database_out_of_fastafile() is deprecated.\n"
            . q{Use generate_database() instead.} );
    return $self->generate_database(@args);
}

sub available_sort_modes {
    my ($self) = @_;
    return (
        ga => 'ascending order of dG',
        gd => 'descending order of dG',
    );
}

sub _get_sequences_from_bio_index {
    my ( $self, $seqid ) = @_;

    $self->is_arrayref_of_size( $seqid, 1 );
    my $indexfile = $self->settings->datapath . q{/}
        . $self->settings->database . '.idx';
    my $idx = Bio::Index::Fasta->new($indexfile);
    my $string;
    my $stringio = IO::String->new($string);
    my $out = Bio::SeqIO->new( -fh => $stringio, -format => 'Fasta' );
    foreach my $seqid ( @{$seqid} ) {
        $out->write_seq( $idx->fetch($seqid) );
    }
    $stringio = IO::String->new($string);
    return Bio::SeqIO->new( -fh => $stringio, -format => 'Fasta' );
}

sub get_alphabet_of_database {
    my ( $self, $db ) = @_;
    my %dbs = $self->get_databases();
    if ( !defined $dbs{$db} ) {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => 'Database not found',
            -value => $db,
        );
    }
    $db = $self->is_word($db);

    open my $ALFILE, '<',
        $self->_cat_path_filename( $self->settings->datapath, $db . '.al1' );
    my $lines = 0;
    while ( my $line = <$ALFILE> ) {
        $lines++;
    }
    close $ALFILE;
    return $lines <= $DNA_ALPHABET_SIZE ? return 'dna' : return 'protein';
}

sub _delete_output {
    my ($self) = @_;
    return 0 if !defined $self->_output_fn;
    return 0 if !-e $self->_output_fn;
    unlink $self->_output_fn
        or $self->throw(
        -class => 'Bio::Root::IOException',
        -text  => q{Can't remove } . $self->_output_fn,
        -value => $OS_ERROR,
        );
    return 1;
}

sub _execute_command {
    my ( $self, $cmd ) = @_;
    $self->_delete_output();
    my ( $tmp_fh, $tmp_fn )
        = tempfile( 'parser_XXXXXXXXXXXX', DIR => $self->settings->tmppath );
    $self->_output_fh($tmp_fh);
    $self->_output_fn($tmp_fn);
    system "$cmd > $tmp_fn";

    return !$CHILD_ERROR;
}

sub _execute_command_and_return_output {
    my ( $self, $cmd ) = @_;
    my ( $writer, $reader, $err );
    my $pid = open3( $writer, $reader, $err, $cmd );
    waitpid $pid, 0;
    my $output = do { local $INPUT_RECORD_SEPARATOR = undef; <$reader> };

    #    my $error = join q{}, <$err>;
    return $output;
}

sub _create_index_and_alphabet_file {
    my ( $self, $filename ) = @_;
    my $in = Bio::SeqIO->new( -file => $filename, -format => 'Fasta' );
    my $alphabet = $in->next_seq()->alphabet();
    if ( !defined $self->features->{PROTEINS} && $alphabet eq 'protein' ) {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => 'Back-end does not support protein data',
            -value => $alphabet,
        );
    }

    # create a vmatch alphabet file
    open my $ALFILE, '>', "$filename.al1";
    if ( $alphabet eq 'dna' ) {
        print ${ALFILE} "aA\ncC\ngG\ntTuU\nnsywrkvbdhmNSYWRKVBDHM\n"
            or $self->_cannot_print("$filename.al1");
    }
    else {
        print ${ALFILE}
            "L\nV\nI\nF\nK\nR\nE\nD\nA\nG\nS\nT\nN\nQ\nY\nW\nP\nH\nM\nC\nXUBZ*\n"
            or $self->_cannot_print("$filename.al1");
    }
    close $ALFILE;

    # create database from directory of fasta files
    my $idx = Bio::Index::Fasta->new(
        -filename   => $filename . '.idx',
        -write_flag => 1
    );
    $idx->make_index( ($filename) );

    return;
}

sub _create_tmp_query_file {
    my ($self)     = @_;
    my $s          = $self->settings;
    my $query_file = 0;
    my ( $tmp_fh, $tmp_query_file );
    if ( defined $s->query_file ) {
        $query_file = $self->is_path( $s->query_file );
    }
    if ( $s->query_isset && $s->query_file_isset ) {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => 'Query and query_file are set. I am confused...',
            -value => $s->query . ' and ' . $s->query_file,
        );
    }
    my @query_seqs = ();

    my $query = q{};
    if ( !$query_file ) {
        $query = $self->_prepare_query();

        ( $tmp_fh, $tmp_query_file ) = tempfile(
            'spatter_XXXXXXXXXXXX',
            DIR    => $s->tmppath,
            UNLINK => !$ENV{BIOGREPDEBUG}
        );    # don't delete when in debug mode, so
              # we can reproduce it.

        # construct a temporary fasta file with the query for vmatch
        my $seqobj = $self->{_query_obj};

        my $outseqio = Bio::SeqIO->new(
            -fh     => $tmp_fh,
            -format => 'fasta'
        );
        $outseqio->write_seq($seqobj);

        push @query_seqs, $seqobj;
        if ( $s->direct_and_rev_com
            && !defined $self->features->{NATIVE_D_A_REV_COM} )
        {
            my $seqobj2 = $seqobj->revcom;
            if ( defined $self->features->{REVCOM_DEFAULT} ) {
                my ($desc)
                    = $seqobj->desc
                    =~ m{\A(.*) \s\(reverse\scomplement\) \z}xms;
                $seqobj2->desc($desc);
            }
            else {
                $seqobj2->desc( $seqobj->desc . ' (reverse complement)' );
            }
            $outseqio->write_seq($seqobj2);
            push @query_seqs, $seqobj2;
        }

    }
    else {
        $tmp_query_file = $query_file;
        my $query_in = Bio::SeqIO->new( -file => $tmp_query_file );
        while ( my $query_seq = $query_in->next_seq ) {
            push @query_seqs, $query_seq;
        }
    }
    $self->_query_seqs(@query_seqs);
    return ( $query, $query_file, $tmp_query_file );
}

sub _prepare_generate_database {
    my ( $self, @args ) = @_;
    my $reftype = reftype $args[0];
    my %args;
    if ( !defined $reftype ) {
        $self->deprecated(
                  'List argument for generate_database() is deprecated. '
                . 'Use hashref instead.' );
        %args = (
            format      => 'Fasta',
            file        => $args[0],
            description => $args[1],
            copy        => 0,
        );
    }
    elsif ( $reftype eq 'HASH' ) {
        %args = %{ $args[0] };
    }
    else {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => 'Invalid arguments. Hashref assumed.',
        );
    }

    if ( !-e $args{file} ) {
        $self->throw(
            -class => 'Bio::Root::BadParameter',
            -text  => 'No such file.',
            -value => $args{file},
        );
    }
    my ($filename) = fileparse( $args{file} );
    $args{basefilename} = $filename;

    if ( !defined $args{format} ) {
        $args{format} = 'Fasta';
    }
    if ( defined $args{prefix_length} ) {

        # untaint pl
        $args{prefix_length} = $self->is_integer( $args{prefix_length} );
    }
    if ( defined $args{datapath} ) {
        $self->settings->datapath( $args{datapath} );
    }
    if ( !defined $self->settings->datapath ) {
        $self->settings->datapath(q{./});
    }
    $self->_copy_fasta_file_and_create_nfo( \%args );

    return %args;
}

sub _parse_regions {
    my ( $self, $args ) = @_;
    my $upstream_seq = q{};
    my $subject_seq = substr $args->{complete_seq}, $args->{subject_begin},
        $args->{subject_end} - $args->{subject_begin};
    my $downstream_seq = q{};

    # initialized for assert below
    my $upstream_begin = $args->{subject_begin};

    if ( $self->settings->upstream > 0 || $self->settings->downstream > 0 ) {

        # coordinates of upstream region, check if available region is
        # as large as requested
        $upstream_begin = $args->{subject_begin} - $self->settings->upstream;
        if ( $upstream_begin < 0 ) {
            $upstream_begin = 0;
        }

        $upstream_seq = substr $args->{complete_seq}, $upstream_begin,
            $args->{subject_begin} - $upstream_begin;

        # and same for the downstream region
        my $downstream_end
            = $args->{subject_end} + $self->settings->downstream;

        if ( $downstream_end > length $args->{complete_seq} ) {
            $downstream_end = length $args->{complete_seq};
        }

        $downstream_seq = substr $args->{complete_seq}, $args->{subject_end},
            $downstream_end - $args->{subject_end};

    }
    my @ret = ( $upstream_seq, $subject_seq, $downstream_seq );
    assert(
        index( $args->{complete_seq}, join( q{}, @ret ), $upstream_begin )
            == $upstream_begin )
        if DEBUG;
    return @ret;
}

1;    # Magic true value required at end of module
__END__

# vim: ft=perl sw=4 ts=4 expandtab