view release on metacpan or  search on metacpan

bin/annotate-ali.pl  view on Meta::CPAN

        my $ali = Bio::MUST::Core::Ali->load($infile);
        $ali->dont_guess if $ARGV_noguessing;
        my $outfile = secure_outfile($infile, $ARGV_out_suffix);
        #### Writing annotated file: $outfile->stringify
        prefix_ids($ali, \%ann_for)->store_fasta($outfile);
    }
}

# TODO: replace by or add some --store-id-mapper option?
# TODO: move into BMC::Ali
sub prefix_ids {
    my $ali     = shift;
    my $ann_for = shift;

    for my $seq ($ali->all_seqs) {
        my $prefix = $ann_for->{$seq->full_id};
        $seq->set_seq_id( $prefix . q{-} . $seq->full_id ) if $prefix;
    }

    return $ali;
}

# TODO: check coherence of option layout with cdhit-tax-filter.pl
# e.g., replace --ann-file by --store-id-mapper

__END__

=pod

=head1 NAME

annotate-ali.pl - Annotate sequences by homology search using BLAST

=head1 VERSION

version 0.252830

=head1 USAGE

    annotate-ali.pl <infiles> --database=<file> --regex=<regex> [options]

=head1 REQUIRED ARGUMENTS

=over

=item <infiles>

Path to input ALI files [repeatable argument].

=for Euclid: infiles.type: readable
    repeatable

=item --ref-file [=] <file> | --database [=] <file>

Path to the FASTA file containing the sequence database.

=for Euclid: file.type: readable

=item --[ref-]regex [=] <regex>

Regular expression for capturing the reference seq id fragment that has to be
used for annotating infile seq ids.

=for Euclid: regex.type: string

=back

=head1 OPTIONS

=over

=item --evalue [=] <float>

E-value threshold for annotating a sequence [default: 1e-10].

=for Euclid: float.type: number
    float.default: 1e-10

=item --identity [=] <number>

Identity threshold for annotating a sequence [default: 0]. When specified as a
fraction between 0 and 1 (included), it is first multiplied by 100 to be
interpreted in percentage.

=for Euclid: number.type: number
    number.default: 0

=item --max-hits [=] <number>

Number of hits to return for each query (BLAST -max_target_seqs option)
[default: 10]. Mostly useful in conjunction with the C<--hit-list> option.

=for Euclid: number.type: number
    number.default: 10

=item --tblastx

Enforce using TBLASTX (instead of the auto-selected BLASTN) when both the
infile and the database contain nucleotide sequences [default: no].

=item --ann-file

Write an annotated version (with prefixed ids) of the infile [default: no].

=item --hit-list

Print a list of id/hit pairs (at the specified E-value and identity thresholds)
instead of the standard annotation report [default: no].

=item --fields [=] <str>...

List of whitespace-separated BLAST fields to be displayed in final report
[default: no].

Valid fields are: percent_identity, hsp_length, mismatches, gaps, query_from,
query_to, hit_from, hit_to, evalue, bit_score, query_start, query_end,
hit_start, hit_end.

=for Euclid: str.type: string, str eq "percent_identity" || str eq "hsp_length" || str eq "mismatches" || str eq "gaps" || str eq "query_from" || str eq "query_to" || str eq "hit_from" || str eq "hit_to" || str eq "evalue" || str eq "bit_score" || st...
    str.default: []