← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to scripts/searchio/fastam9_to_table.PLS

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

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