view release on metacpan or  search on metacpan

lib/Bio/MUST/Apps/TwoScalp/Seq2Seq.pm  view on Meta::CPAN

has 'ali' => (
    is       => 'ro',
    isa      => 'Bio::MUST::Core::Ali',
    required => 1,
    coerce   => 1,
);

has 'lookup' => (
    is       => 'ro',
    isa      => 'Bio::MUST::Core::IdList',
    init_arg => undef,
    writer   => '_set_lookup',
);

has 'new_ali' => (
    is       => 'ro',
    isa      => 'Bio::MUST::Core::Ali',
    init_arg => undef,
    writer   => '_set_new_ali',
);

has 'integrator' => (
    is       => 'ro',
    isa      => 'Bio::MUST::Apps::SlaveAligner::Local',
    init_arg => undef,
    writer   => '_set_integrator',
);

has 'blastdb' => (
    is       => 'ro',
    isa      => 'Bio::MUST::Drivers::Blast::Database::Temporary',
    init_arg => undef,
    writer   => '_set_blastdb',
);

has 'query_seqs' => (
    is       => 'ro',
    isa      => 'Bio::MUST::Core::Ali::Temporary',
    init_arg => undef,
    writer   => '_set_query_seqs',
);


sub _align_seqs {
    my $self = shift;

    my $args->{-outfmt} = 5;
       $args->{-max_target_seqs} = 5;

    my $blastdb = $self->blastdb;
    my $query_seqs = $self->query_seqs;
    my $parser = $blastdb->blast($query_seqs, $args);
    #### [S2S] XML BLASTP/N: $parser->filename

    my $bo = $parser->blast_output;
    return unless $bo;

    my $sort_method = $self->single_hsp ? 'score' : 'hit_start';

    QUERY:
    for my $query ( $bo->all_iterations ) {

        my $query_id = $query_seqs->long_id_for( $query->query_def );
        ##### [S2S] Aligning: $query_id

        my @templates;
        my @template_seqs;
        my $best_coverage = 0;

        TEMPLATE:
        for my $template ( $query->all_hits ) {

            # compute query coverage by template HSPs
            my $mask = SeqMask->empty_mask( $query->query_len );
            $mask->mark_block( $_->query_start, $_->query_end )
                for $template->all_hsps;
            my $coverage = $mask->coverage;

            last TEMPLATE if $coverage
                < $best_coverage * $self->coverage_mul;
            $best_coverage = $coverage;

            # fetch template full_id
            my $template_id = SeqId->new(
                full_id => $blastdb->long_id_for( $template->def )
                        // $blastdb->long_id_for( $template->id  )
            );          # workaround to accommodate change in BLAST XML report

            ###### [S2S] template: $template_id->full_id
            ###### [S2S] coverage: $coverage

            # fetch and cache aligned template seq from Ali
            push @templates, $template;
            push @template_seqs, $self->ali->get_seq(
                $self->lookup->index_for( $template_id->full_id )
            );
        }

        unless (@templates) {
            ##### [S2S] skipped alignment due to lack of suitable template
            next QUERY;
        }

        @templates = ( $templates[-1] ) if $self->single_hsp;

        TEMPLATE:
        for my $template (@templates) {

            # use each template in turn for BLAST alignment
            my $template_seq = shift @template_seqs;
            last TEMPLATE unless $template_seq;

            # sort HSPs by descending start coordinate on template or by
            # descending score depending on --single-hsp option
            my @hsps = sort {
                $b->$sort_method <=> $a->$sort_method
            } $template->all_hsps;

            HSP:
            for my $hsp (@hsps) {