~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to scripts/searchio/parse_hmmsearch.PLS

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
#!/usr/bin/perl
2
 
 
3
 
use strict;
4
 
 
5
 
=head1 NAME
6
 
 
7
 
parse_hmmsearch - parse single/multiple HMMSEARCH results file(s) with
8
 
                  different output options
9
 
 
10
 
=head1 SYNOPSIS
11
 
 
12
 
parse_hmmsearch [--po] [--ps] -s hmmsearch_file
13
 
 
14
 
parse_hmmsearch [--po] [--ps] -m index_file
15
 
 
16
 
=head1 DESCRIPTION
17
 
 
18
 
=head2 Mandatory Options:
19
 
 
20
 
  -s  HMMSEARCH file to parse.
21
 
  -m  INDEX file that contains a list of HMMSEARCH files for multiple
22
 
      parsing.
23
 
 
24
 
=head2 Special Options:
25
 
 
26
 
  --po    Print only the hits that have positive scores.
27
 
  --ps    Print the total of positive scores found.
28
 
  --help  Show this documentation.
29
 
 
30
 
=head1 FEEDBACK
31
 
 
32
 
=head2 Mailing Lists
33
 
 
34
 
User feedback is an integral part of the evolution of this and other
35
 
Bioperl modules. Send your comments and suggestions preferably to the
36
 
Bioperl mailing list. Your participation is much appreciated.
37
 
 
38
 
  bioperl-l@bioperl.org                  - General discussion
39
 
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
40
 
 
41
 
=head2 Reporting Bugs
42
 
 
43
 
Report bugs to the Bioperl bug tracking system to help us keep track
44
 
of the bugs and their resolution. Bug reports can be submitted via the
45
 
web:
46
 
 
47
 
  https://redmine.open-bio.org/projects/bioperl/
48
 
 
49
 
=head1 AUTHOR
50
 
 
51
 
  Mauricio Herrera Cuadra <mauricio at open-bio.org>
52
 
 
53
 
=cut
54
 
 
55
 
# Modules, pragmas and variables to use
56
 
use Bio::SearchIO;
57
 
use Getopt::Long;
58
 
use vars qw($opt_s $opt_m $opt_po $opt_ps $opt_help);
59
 
 
60
 
# Gets options from the command line
61
 
GetOptions qw(-s:s -m:s --po --ps --help);
62
 
 
63
 
# Print documentation if help switch was given
64
 
exec('perldoc', $0) and exit() if $opt_help;
65
 
 
66
 
# If no mandatory options are given prints an error and exits
67
 
if (!$opt_s && !$opt_m) {
68
 
    print "ERROR: No HMMSEARCH or INDEX file has been specified.\n       Use
69
 
'--help' switch for documentation.\n" and exit();
70
 
} elsif ($opt_s && $opt_m) {
71
 
    print "ERROR: You must select only one option (-s or -m) for input.\n      
72
 
Use '--help' switch for documentation.\n" and exit();
73
 
}
74
 
 
75
 
# Initializes a counter for the domain positive scores if the option
76
 
# was given
77
 
my $pos_scores = 0 if $opt_ps;
78
 
 
79
 
# If single file mode was selected
80
 
if ($opt_s) {
81
 
    parse_hmmsearch($opt_s);
82
 
 
83
 
    # Prints the total domain positive scores if the option was given
84
 
    if ($opt_ps) {
85
 
        print "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86
 
- - - -\n";
87
 
        print "Total domain positive scores: $pos_scores\n";
88
 
    }
89
 
 
90
 
# If multiple files mode was selected
91
 
} elsif ($opt_m) {
92
 
 
93
 
    # Opens the INDEX file sent as input
94
 
    open(FH, "<", $opt_m) or die("Unable to open INDEX file: $opt_m ($!)");
95
 
 
96
 
    # Cycle that extracts one line for every loop until finding the
97
 
    # end of file
98
 
    while (my $line = <FH>) {
99
 
 
100
 
        # Deletes the new line characters from the line
101
 
        chomp $line;
102
 
 
103
 
        # Parses the result file in turn
104
 
        parse_hmmsearch($line);
105
 
        print "= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
106
 
= = = =\n";
107
 
    }
108
 
 
109
 
    # Prints the total domain positive scores if the option was given
110
 
    print "Total domain positive scores: $pos_scores\n" if $opt_ps;
111
 
 
112
 
    # Closes INDEX files
113
 
    close(FH);
114
 
}
115
 
 
116
 
# Exits the program
117
 
exit();
118
 
 
119
 
# Subroutine that parses a HMMSEARCH results file
120
 
sub parse_hmmsearch {
121
 
 
122
 
    # Gets the parameters sent to the function
123
 
    my ($file) = @_;
124
 
 
125
 
    # Creates a new Bio::SearchIO object
126
 
    my $in = new Bio::SearchIO(
127
 
        -format => 'hmmer',
128
 
        -file   => $file,
129
 
    );
130
 
 
131
 
    # Loops through the results file
132
 
    while (my $result = $in->next_result()) {
133
 
 
134
 
        # Prints program name and version (these are values from
135
 
        # Bio::Search::Result::GenericResult methods)
136
 
        print $result->algorithm(), " ", $result->algorithm_version(), "\n";
137
 
        print "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138
 
- - - -\n";
139
 
 
140
 
        # Prints HMM file and sequence database (these are values from
141
 
        # Bio::Search::Result::HMMERResult methods)
142
 
        print "HMM file:\t\t\t", $result->hmm_name(), "\n";
143
 
        print "Sequence database:\t\t", $result->sequence_file(), "\n";
144
 
        print "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145
 
-\n";
146
 
 
147
 
        # Prints some values from Bio::Search::Result::GenericResult
148
 
        # methods
149
 
        print "Query HMM:\t\t\t", $result->query_name(), "\n";
150
 
        print "Accession:\t\t\t", $result->query_accession(), "\n";
151
 
        print "Description:\t\t\t", $result->query_description(), "\n";
152
 
        print "Total hits:\t\t\t", $result->num_hits(), "\n";
153
 
 
154
 
        # Loops through the sequence in turn
155
 
        while (my $hit = $result->next_hit()) {
156
 
 
157
 
            # If only positive scores option was given and the score
158
 
            # in turn is greater than zero
159
 
            if ($opt_po) {
160
 
                printHits($hit) if ($hit->score() >= 0);
161
 
 
162
 
            # Prints all hits otherwise
163
 
            } else {
164
 
                printHits($hit);
165
 
            }
166
 
        }
167
 
    }
168
 
}
169
 
 
170
 
# Subroutine that prints the values from a Bio::Search::Hit::HitI
171
 
# object
172
 
sub printHits {
173
 
 
174
 
    # Gets the parameters sent to the function
175
 
    my ($hit) = @_;
176
 
 
177
 
    # Prints some values from Bio::Search::Hit::HitI methods
178
 
    print "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n";
179
 
    print "Hit ", $hit->rank(), "\n";
180
 
    print "Sequence:\t\t\t", $hit->name(), "\n";
181
 
    print "Description:\t\t\t", $hit->description(), "\n";
182
 
    print "Score:\t\t\t\t", $hit->score(), "\n";
183
 
    print "E-value:\t\t\t", $hit->significance(), "\n";
184
 
    print "Number of domains:\t\t", $hit->num_hsps(), "\n";
185
 
 
186
 
    # Loops through the domain in turn
187
 
    while (my $hsp = $hit->next_hsp()) {
188
 
 
189
 
        # Prints some values from Bio::Search::HSP::HSPI methods
190
 
        print "   - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n";
191
 
        print "   Domain:\t\t\t", $hsp->rank(), " of ", $hit->num_hsps(), "\n";
192
 
        print "   seq-f:\t\t\t", $hsp->start('hit'), "\n";
193
 
        print "   seq-t:\t\t\t", $hsp->end('hit'), "\n";
194
 
        print "   hmm-f:\t\t\t", $hsp->start(), "\n";
195
 
        print "   hmm-t:\t\t\t", $hsp->end(), "\n";
196
 
        print "   score:\t\t\t", $hsp->score(), "\n";
197
 
        $pos_scores++ if ($hsp->score() >= 0) && $opt_ps;
198
 
        print "   E-value:\t\t\t", $hsp->evalue(), "\n";
199
 
        my $hmm_string = $hsp->query_string();
200
 
        $hmm_string =~ s/<-\*$//;
201
 
        print "   hmm string:\t\t\t", $hmm_string, "\n";
202
 
        print "   homology string:\t\t", $hsp->homology_string(), "\n";
203
 
        print "   hit string:\t\t\t", $hsp->hit_string(), "\n";
204
 
    }
205
 
}