/usr/bin/bp_fastam9_to_table is in bioperl 1.6.924-3.
This file is owned by root:root, with mode 0o755.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 | #!/usr/bin/perl
eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
if 0; # not running under some shell
=head1 NAME
fastm9_to_table - turn FASTA -m 9 output into NCBI -m 9 tabular output
=head1 SYNOPSIS
fastm9_to_table [-e evaluefilter] [-b bitscorefilter] [--header] [-o outfile] inputfile1 inputfile2 ...
=head1 DESCRIPTION
Command line options:
-e/--evalue evalue -- filter by evalue
-b/--bitscore bitscore -- filter by bitscore
--header -- boolean flag to print column header
-o/--out -- optional outputfile to write data,
otherwise will write to STDOUT
-h/--help -- show this documentation
Not technically a SearchIO script as this doesn't use any Bioperl
components but is a useful and fast. The output is tabular output
with the standard NCBI -m9 columns.
queryname
hit name
percent identity
alignment length
number mismatches
number gaps
query start (if on rev-strand start > end)
query end
hit start (if on rev-strand start > end)
hit end
evalue
bit score
Additionally 3 more columns are provided
fasta score
sw-score
percent similar
query length
hit length
query gaps
hit gaps
=head1 AUTHOR - Jason Stajich
Jason Stajich jason_at_bioperl-dot-org
=cut
use strict;
use warnings;
use Getopt::Long;
my $hitsection = 0;
my %data;
my ($evalue,$bitscore,$header,$outfile) = ( 10, 0);
GetOptions(
'b|bitscore|bits:f' => \$bitscore,
'e|evalue:f' => \$evalue,
'header' => \$header,
'o|out|outfile:s' => \$outfile,
'h|help' => sub { exec('perldoc',$0); exit; }
);
my $outfh;
if( $outfile ) {
open $outfh, '>', $outfile or die "Could not write file '$outfile': $!\n";
} else {
$outfh = \*STDOUT;
}
# query start -- an0
# query en -- ax0
# hit start -- an1
# hit end -- ax1
my @fields = qw(qname hname percid alen mmcount gapcount
qstart qend hstart hend evalue score bits fs sw-score
percsim qlen hlen qgap hgap);
print $outfh "#",uc(join("", map{ sprintf("%-10s",$_) } @fields)), "\n" if $header;
while(<>) {
my $linestr = $_;
if( /^\s*\d+>>>(\S+).+/ ) {
$data{'qname'} = $1;
if( /\-?\s+(\d+)\s+(aa|nt)\s+$/ ){
$data{'qlen'} = $1;
}
} elsif( $hitsection && (/^>>>\Q$data{'qname'}/ || /^>>>/) ) {
$hitsection = 0;
} elsif( /^The best scores are:/ ) {
$hitsection = 1;
} elsif( /^\s+$/ ) {
} elsif( $hitsection ) {
if( s/^(\S+)\s+(.+)\(\s*(\d+)\)\s+// ) {
my ($hit, $desc,$hitlen) = ($1,$2,$3);
my ($dir) = ( s/^\[(r|f)\]\s+// );
my @line = split(/\s+/,$_);
$data{'hname'} = $hit;
$data{'hlen'} = $hitlen;
$data{'score'} = shift @line;
$data{'bits'} = shift @line;
$data{'evalue'} = shift @line;
$data{'percid'} = shift @line;
$data{'percsim'} = shift @line;
$data{'sw-score'} = shift @line;
$data{'alen'} = shift @line;
$data{'qstart'} = shift @line;
$data{'qend'} = shift @line;
$data{'pn0'} = shift @line; # pn0
$data{'px0'} = shift @line; # px0
$data{'hstart'} = shift @line; # an1
$data{'hend'} = shift @line; # ax1
$data{'pn1'} = shift @line; # pn1
$data{'px1'} = shift @line; # px1
# query + hit gaps
$data{'qgap'} = shift @line;
$data{'hgap'} = shift @line;
$data{'gapcount'} = $data{'qgap'} + $data{'hgap'};
$data{'fs'} = shift @line;
$data{'mmcount'} = $data{'alen'} - ( int($data{'percid'} * $data{'alen'}) + $data{'gapcount'});
#next if( $data{'evalue'} > $evalue ||
# $data{'bits'} < $bitscore );
for ( $data{'percid'}, $data{'percsim'} ) {
$_ = sprintf("%.2f",$_*100);
}
print $outfh join( "\t",map { $data{$_} } @fields),"\n";
} else {
# print STDERR "unrecognized line \n$linestr";
}
} else {
# warn("skipping a line like this: $_");
}
}
|