scripts/multi_hmmsearch.PLS

   1 #!/usr/bin/perl -w
   2 # $Id: multi_hmmsearch.PLS,v 1.3 2006-07-04 22:23:36 mauricio Exp $
   3
   4 use strict;
   5
   6 =head1 NAME
   7
   8 multi_hmmsearch - perform a hmmsearch into multiple FASTA files using
   9                   an INDEX file
  10
  11 =head1 SYNOPSIS
  12
  13 multi_hmmsearch -p hmm_file [-i] -f index_file
  14
  15 =head1 DESCRIPTION
  16
  17 Not technically a Bio::Tools::Run script as this doesn't use any
  18 Bioperl or Bioperl-run components but it's useful.
  19
  20 =head2 Mandatory Options:
  21
  22   -p  HMM profile to use for the search.
  23   -f  INDEX file that contains a list of FASTA files for the multiple
  24       search.
  25
  26 =head2 Special Options:
  27
  28   -i  Create a new index file with the resulting hmms files. This is
  29       useful if you want to pass this list as input arguments into
  30       another programs.
  31   -h  Show this documentation.
  32
  33 =head1 FEEDBACK
  34
  35 =head2 Mailing Lists
  36
  37 User feedback is an integral part of the evolution of this and other
  38 Bioperl modules. Send your comments and suggestions preferably to the
  39 Bioperl mailing list. Your participation is much appreciated.
  40
  41   bioperl-l@bioperl.org                  - General discussion
  42   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  43
  44 =head2 Reporting Bugs
  45
  46 Report bugs to the Bioperl bug tracking system to help us keep track
  47 of the bugs and their resolution. Bug reports can be submitted via the
  48 web:
  49
  50   http://redmine.open-bio.org/projects/bioperl/
  51
  52 =head1 AUTHOR
  53
  54   Mauricio Herrera Cuadra <mauricio-at-open-bio.org>
  55
  56 =cut
  57
  58 # Modules, pragmas and variables to use
  59 use Getopt::Long;
  60 use vars qw($opt_p $opt_i $opt_f $opt_h $index_file);
  61
  62 # Gets options from the command line
  63 GetOptions qw(-p=s -i -f=s -h);
  64
  65 # Print documentation if help switch was given
  66 exec('perldoc', $0) and exit() if $opt_h;
  67
  68 # If no mandatory options are given prints an error and exits
  69 if (!$opt_p) {
  70     print "ERROR: No HMM profile has been specified.\n       Use '-h' switch
  71 for documentation.\n" and exit();
  72 } elsif (!$opt_f) {
  73     print "ERROR: No INDEX file has been specified.\n       Use '-h' switch for
  74 documentation.\n" and exit();
  75 }
  76
  77 # Locates hmmsearch in the filesystem
  78 my $hmmsearch = `which hmmsearch`;
  79 chomp $hmmsearch;
  80
  81 # Creates a directory for writing the resulting files
  82 mkdir("multi", 0755) unless -e "multi" and -d "multi";
  83
  84 # Creates a new INDEX file if the option was given
  85 if ($opt_i) {
  86     my $prefix = $opt_f;
  87     $prefix =~ s/\.INDEX$//;
  88     $index_file = "$prefix.hmms.INDEX";
  89     open(HMMSINDEX, ">", $index_file) or die("Unable to create file:
  90 $index_file ($!)");
  91 }
  92
  93 # Opens the INDEX file sent as input
  94 open(FH, "<", $opt_f) or die("Unable to open INDEX file: $opt_f ($!)");
  95 print "==> Opening INDEX file:\t\t\t\t$opt_f\n";
  96 print "==> HMM profile file is:\t\t\t$opt_p\n";
  97
  98 # Cycle that extracts one line for every loop until finding the end of
  99 # file
 100 while (my $line = <FH>) {
 101
 102     # Deletes the new line characters from the line
 103     chomp $line;
 104
 105     # Gets the name for the result file
 106     my $out = $line;
 107     $out =~ s/^split\///;
 108     $out =~ s/\.faa$//;
 109
 110     # Performs the hmmsearch for the FASTA file in turn
 111     print "--> Performing hmmsearch in file:\t\t$line\n";
 112     system("$hmmsearch $opt_p $line > multi/$out.hmms");
 113     print "==> hmmsearch results stored in file:\t\tmulti/$out.hmms\n";
 114
 115     # Prints the result file name into the new INDEX file if the
 116     # option was given
 117     print HMMSINDEX "multi/$out.hmms\n" if $opt_i;
 118 }
 119
 120 # Closes INDEX files
 121 close(FH);
 122 if ($opt_i) {
 123     print "==> New INDEX stored in file:\t\t\t$index_file\n";
 124     close(HMMSINDEX);
 125 }
 126
 127 # Exits the program
 128 exit();