view release on metacpan or  search on metacpan

lib/Bio/Gonzales/Domain/Identification/HMMER.pm  view on Meta::CPAN

    ( my $basename_no_suffix = file($file)->basename ) =~ s/\.\w+?$//;
    return $basename_no_suffix;
}

sub _gen_intermediate_result_name {
    my ( $self, $suffix ) = @_;

    return File::Spec->catfile( $self->result_dir, $self->{_intermediate_result_id}++ . "_" . $suffix );
}

sub _catmyfile {
    my ( $self, $filename ) = @_;

    return File::Spec->catfile( $self->result_dir, $filename );
}

=head2 $i->unlink_destination_files

Deletes all destination files (the ones required in the constructor, except the cache file

=cut

=head2 _transform_hmm_hits

Transforms the best hits result from
Bio::Gonzales::Util::SearchIO::HMMResult->get_best_hits to a more accessible structure
and determines the maximum possible spanning region. This function also builds
up a cache of sequence-id_size, found-domain and position cache for faster access
later on

=cut

sub _transform_hmm_hits {
    my ( $self, $best_hits, $sequence_file ) = @_;

    open my $cache, '>>', $self->discovered_cache_file
        or croak "Can't open filehandle: $!";

    my %result_for_sequence;
    for my $k ( @{$best_hits} ) {

        #write to cache
        say $cache $self->_create_cache_string( $sequence_file, $k )
            unless ( $self->from_cache );

        #find maximum possible spanning region of a group in one sequence
        $result_for_sequence{ $k->{seq_id} } //= Bio::Gonzales::Util::FunCon::Domains::Identification::HMMER::SeqMarks->new(
            num_marks => scalar @{ $self->domain_groups } );
        $self->_update_sequence_mark( $result_for_sequence{ $k->{seq_id} }, $k );
    }
    $cache->close;

    return \%result_for_sequence;
}

sub _create_cache_string {
    my ( $self, $sequence_file, $k ) = @_;

    return join "\t",
        (
        _basename_no_suffix($sequence_file) . "_" . stat($sequence_file)->size,
        $k->{seq_id}, $k->{hmm_acc}, $k->{hmm_score}, $k->{from}, $k->{to}
        );
}

sub _update_sequence_mark {
    my ( $self, $seq_result, $best_hit ) = @_;

    my @domain_groups = @{ $self->domain_groups };

    for ( my $i = 0; $i < @domain_groups; $i++ ) {
        $seq_result->update_mark( $i, $best_hit->{from}, $best_hit->{to} )
            if (
            any { $best_hit->{hmm_acc} eq $_ }
            keys %{ $domain_groups[$i] }
            );
    }
}

sub _get_cached_hits {
    my ( $self, $sequence_file ) = @_;

    my @best_hits;

    my $seq_file_id = _basename_no_suffix($sequence_file) . "_" . stat($sequence_file)->size;

    open my $cache, '<', $self->discovered_cache_file
        or croak "Can't open filehandle: $!";
    while ( my $l = <$cache> ) {
        my @rows = split /\t/, $l;
        push @best_hits,
            {
            seq_id    => $rows[1],
            hmm_acc   => $rows[2],
            hmm_score => $rows[3],
            from      => $rows[4],
            to        => $rows[5]
            }
            if ( $rows[0] eq $seq_file_id );
    }
    return \@best_hits;
}

sub identify {
    my ( $self, $sequence_file, $tag ) = @_;

    #run hmmsearch and get the best hits for each domain and seq_id

    my $best_hits;
    unless ( $self->from_cache ) {
        $best_hits
            = Bio::Gonzales::Util::SearchIO::HMMResult->new( file => $self->_run_hmmsearch($sequence_file) )->get_best_hits();

    } elsif ( -f $self->discovered_cache_file ) {
        $best_hits = $self->_get_cached_hits($sequence_file);
    } else {
        croak "Cache file " . $self->discovered_cache_file . " not found";
    }

    #use standard groups if attribute not set, one group with all domains
    unless ( @{ $self->domain_groups } > 0 ) {
        $self->_create_standard_domain_groups($best_hits);
    }

    my $result_for_sequence = $self->_transform_hmm_hits( $best_hits, $sequence_file );

    #open input sequence file
    my $snf2 = Bio::SeqIO->new(
        -format => 'fasta',
        -file   => $sequence_file,
    );

    #open file for spanning domains, if set
    my $domain_spanning_region;
    $domain_spanning_region = Bio::SeqIO->new(
        -format => 'fasta',
        -file   => ">>" . $self->domain_spanning_region_file,
    ) if ( $self->domain_spanning_region_file );

    my $domain_spanning_region_masked;
    $domain_spanning_region_masked = Bio::SeqIO->new(
        -format => 'fasta',
        -file   => ">>" . $self->domain_spanning_region_masked_file,
    ) if ( $self->domain_spanning_region_masked_file );