—#!/usr/bin/perl -w
use
POSIX;
use
Pod::Usage;
use
FAST;
use
FAST::Bio::SeqIO;
use
strict;
$VERSION
=
$FAST::VERSION
;
$DESC
=
"annotate sequence lengths to descriptions"
;
$NAME
= $0;
$NAME
=~ s/^.*\///;
$COMMAND
=
join
" "
,
$NAME
,
@ARGV
;
$DATE
= POSIX::strftime(
"%c"
,
localtime
());
## DEFAULT OPTION VALUES
my
$def_format
=
$FAST::DEF_FORMAT
;
#7/1/13 "fasta";
my
$def_logname
=
$FAST::DEF_LOGNAME
;
#7/1/13 "FAST.log.txt";
my
$def_join_string
=
$FAST::DEF_JOIN_STRING
;
## OPTION VARIABLES
my
$man
=
undef
;
# --man
my
$help
=
undef
;
# -h
my
$moltype
=
undef
;
# -m, in case bioperl can't tell
my
$format
=
$def_format
;
# --format
my
$log
=
undef
;
# -l
my
$logname
=
$def_logname
;
# -L
my
$comment
=
undef
;
# -C
my
$table
=
undef
;
# -t
my
$join
=
$def_join_string
;
# -j
GetOptions(
'help|h'
=> \
$help
,
'man'
=> \
$man
,
'moltype|m=s'
=>
sub
{
my
(
undef
,
$val
) =
@_
;
die
"$NAME: --moltype or -m option argument must be \"dna\", \"rna\" or \"protein\"\n"
unless
$val
=~ /dna|rna|protein/i;
$moltype
=
$val
;
},
'format=s'
=> \
$format
,
'log|l'
=> \
$log
,
'logname|L=s'
=> \
$logname
,
'comment|C=s'
=> \
$comment
,
'table|t'
=> \
$table
,
'join|j=s'
=> \
$join
,
'fastq|q'
=>
sub
{
$format
=
'fastq'
;},
)
or
exit
(1);
$format
=
lc
(
$format
);
pod2usage(
-verbose
=> 1)
if
$help
;
pod2usage(
-verbose
=> 2)
if
$man
;
my
$fromSTDIN
= ((-t STDIN) ? false : true);
pod2usage(
"$NAME: Requires at least one argument FILE [FILE2…FILEN] unless input is from STDIN.\n"
)
if
((-t STDIN) && (
@ARGV
== 0));
pod2usage(
"$NAME: Requires exactly zero arguments if input is from STDIN. Try $NAME -h for help.\n"
)
if
(!(-t STDIN) && (
@ARGV
!= 0));
&FAST::log
(
$logname
,
$DATE
,
$COMMAND
,
$comment
,
$fromSTDIN
)
if
(
$log
);
my
$OUT
= FAST::Bio::SeqIO->newFh(
-fh
=>
*STDOUT
{IO},
'-format'
=>
$format
);
my
$IN
;
unless
(
@ARGV
) {
if
(
$moltype
) {
$IN
= FAST::Bio::SeqIO->new(
-fh
=>
*STDIN
{IO},
'-format'
=>
$format
,
'-alphabet'
=>
$moltype
);
}
else
{
$IN
= FAST::Bio::SeqIO->new(
-fh
=>
*STDIN
{IO},
'-format'
=>
$format
);
}
}
if
(
$join
eq
'\t'
){
$join
=
"\t"
;
}
my
$numseqs
;
my
$numchars
;
my
$file
;
my
$numseqsfile
;
my
$numcharsfile
;
while
(
$IN
or
@ARGV
) {
if
(
@ARGV
) {
$file
=
shift
(
@ARGV
);
$numseqsfile
=
$numcharsfile
= 0;
unless
(-e
$file
) {
warn
"$NAME: Could not find file $file. Skipping.\n"
;
next
;
}
elsif
(
$moltype
) {
$IN
= FAST::Bio::SeqIO->new(
-file
=>
$file
,
'-format'
=>
$format
,
'-alphabet'
=>
$moltype
);
}
else
{
$IN
= FAST::Bio::SeqIO->new(
-file
=>
$file
,
'-format'
=>
$format
);
}
}
else
{
undef
$file
;
}
if
(
$IN
) {
while
(
my
$seq
=
$IN
->next_seq()) {
my
$length
=
$seq
->
length
();
if
(
$table
) {
my
$id
=
$seq
->id();
join
$join
,
$id
,
$length
;
"\n"
;
}
else
{
my
$olddesc
=
$seq
->desc();
$seq
->desc(
join
$join
,
$olddesc
,
"length:$length"
);
$OUT
$seq
;
}
}
undef
$IN
;
}
}
__END__
=head1 NAME
B<faslen> -- annotate sequence lengths to descriptions
=head1 SYNOPSIS
B<faslen> [MULTIFASTA-FILE...]
=head1 DESCRIPTION
B<faslen> takes sequence or alignment data on input, computes sequence
lengths, and annotates sequence record descriptions with their lengths
in the format "length:<value>." All sequence characters get counted,
including gap characters. Identifiers and lengths may optionally be
printed in a table to STDOUT.
Options specific to B<faslen>:
B<-t>, B<--table> print lengths with ids in a table
B<-j>, B<--join>=<string> use <string> to join data
Options general to FAST:
B<-h>, B<--help> print a brief help message
B<--man> print full documentation
B<--version> print version
B<-l>, B<--log> create/append to logfile
B<-L>, B<--logname>=<string> use logfile name <string>
B<-C>, B<--comment>=<string> save comment <string> to log
B<--format>=<format> use alternative format for input
B<--moltype>=<[dna|rna|protein]> specify input sequence type
B<-q>, B<--fastq> use fastq format as input and output
=head1 INPUT AND OUTPUT
B<faslen> is part of FAST, the FAST Analysis of Sequences Toolbox, based
on Bioperl. Most core FAST utilities expect input and return output in
multifasta format. Input can occur in one or more files or on
STDIN. Output occurs to STDOUT.
=head1 OPTIONS
=over 8
=item B<-t>,
B<--table>
Output data in a table to standard output. By default, length data is
annotated into descriptions.
=item B<-j>,
B<--join=[string]>
Use [string] to append length data in descriptions. Use with argument
"\t" to indicate a tab-character.
=item B<-h>,
B<--help>
Print a brief help message and exit.
=item B<--man>
Print the manual page and exit.
=item B<--version>
Print version information and exit.
=item B<-l>,
B<--log>
Creates, or appends to, a generic FAST logfile in the current working
directory. The logfile records date/time of execution, full command
with options and arguments, and an optional comment.
=item B<-L [string]>,
B<--logname=[string]>
Use [string] as the name of the logfile. Default is "FAST.log.txt".
=item B<-C [string]>,
B<--comment=[string]>
Include comment [string] in logfile. No comment is saved by default.
=item B<--format=[format]>
Use alternative format for input. See man page for "B<faslen>" for
allowed formats. This is for convenience; the FAST tools are designed
to exchange data in Fasta format, and "fasta" is the default format
for this tool.
=item B<-m [dna|rna|protein]>,
B<--moltype=[dna|rna|protein]>
Specify the type of sequence on input (should not be needed in most
cases, but sometimes Bioperl cannot guess and complains when
processing data).
=item B<-q>
B<--fastq>
use fastq format as input and output.
=back
=head1 EXAMPLES
Print length of all sequences:
=over 8
B<faslen> sequences.fas > measured.fas
=back
Print tab-delimited table of lengths:
=over 8
B<faslen> -tj "\t" sequences.fas
=back
Filtering sequences to have length >= 300:
=over 8
B<faslen> sequences.fas | B<fasfilter> -t length 300..
=back
Sorting sequences by length:
=over 8
B<faswc> sequences.fas | B<fassort> -t length
=back
=head1 SEE ALSO
=over 8
=item C<man FAST>
=item C<perldoc FAST>
Introduction and cookbook for FAST
=item L<The FAST Home Page|http://compbio.ucmerced.edu/ardell/FAST>"
=back
=head1 CITING
If you use FAST, please cite I<Lawrence et al. (2015). FAST: FAST Analysis of
Sequences Toolbox.> and Bioperl I<Stajich et al.>.
=cut