forester/archive/perl/bootstrapSelector.pl

   1 #!/usr/bin/perl -w
   2 #
   3 # bootstrapSelector.pl
   4 # --------------------
   5 #
   6 # Copyright (C) 2001 Washington University School of Medicine
   7 # and Howard Hughes Medical Institute
   8 # All rights reserved
   9 #
  10 # Author: Christian M. Zmasek
  11 #         zmasek@genetics.wustl.edu
  12 #         http://www.genetics.wustl.edu/eddy/people/zmasek/
  13 #
  14 # Created: 04/06/01
  15 #
  16 # Last modified 09/24/01
  17 #
  18 #
  19 # Objective. Selection of RIO analysis results with top ortholgy
  20 #            bootstrap values greater or less than a threshold.
  21 #
  22 # Usage: "bootstrapSelector.pl <threshold options> <infile = Xrio.pl-output> <outfile>"
  23 # Options: "l" for "less or equal" ("grater or equal" is default)
  24 #          "c" for "all hits must meet threshold in case of
  25 #              multiple copies of the same domain in the query"
  26 #              (default: "at least one")
  27 # Example: "bootstrapSelector.pl 95lc OUTFILE_At_1 At_1_out"
  28 #
  29 # Important. The result of this is meaningful ONLY if the thresholds
  30 #            for output of the RIO analysis are set to zero (L=0 R=0).
  31 #
  32 #
  33 # Format for infile:
  34 #
  35 # ...
  36 #
  37 # # ############################################################################
  38 # # Annotation: B0511.6 CE17345   helicase (ST.LOUIS) TR:O61815 protein_id:AAC17654.1
  39 # # HMM       : ABC_tran
  40 # # score     : -59.6
  41 # # E-value   : 1.1
  42 # # Query has not been aligned (score lower than gathering cutoff).
  43 # # ############################################################################
  44 #
  45 #
  46 # # ############################################################################
  47 # # Annotation: B0511.7 CE17346    (ST.LOUIS) TR:O61817 protein_id:AAC17655.1
  48 # # HMM       : FHA
  49 # # score     : 71.6
  50 # # E-value   : 1.7e-17
  51 # RIO - Resampled Inference of Orthologs
  52 # Version: 1.000
  53 # ------------------------------------------------------------------------------
  54 # Alignment file: /tmp/Xriopl9846081980/Full-FHA
  55 # Alignment     : FHA domain
  56 # HMM           : FHA
  57 # Query file    : /tmp/Xriopl9846081980/__queryfile__
  58 # ==============================================================================
  59 #
  60 # Query         : CE17346.FHA_CAEEL/45-114
  61 #
  62 # Number (in %) of observed orthologies (o) and super orthologies (s) to query
  63 # in bootstrapped trees, evolutionary distance to query:
  64 #
  65 # Sequence              Description                                                                  # o[%] s[%]  distance
  66 # --------              -----------                                                                  ---- ----  --------
  67 # YC67_MYCTU/308-372    -                                                                              20   14  1.577840
  68 # FRAH_ANASP/204-277    FRAH PROTEIN.                                                                  17   16  1.532670
  69 # ABA2_NICPL/557-633    ZEAXANTHIN EPOXIDASE PRECURSOR (EC 1.14.-.-).                                  14   11  1.885700
  70 # ABA2_LYCES/563-639    ZEAXANTHIN EPOXIDASE PRECURSOR (EC 1.14.-.-).                                  14   11  2.140000
  71 #
  72 #
  73 #
  74 # Distance values (based on ML branch length values on consensus tree)
  75 # --------------------------------------------------------------------
  76 # Given the thresholds for distance calculations:
  77 # No sequence is considered orthologous to query.
  78 #
  79 # ...
  80
  81
  82
  83 use strict;
  84
  85 my $VERSION              = 1.000;
  86 my $threshold            = 0;
  87 my $infile               = "";
  88 my $outfile              = "";
  89 my $summary_outfile      = "";
  90 my $return_line          = "";
  91 my $identifier           = "";
  92 my $top1                 = "";
  93 my $analysis_performed   = 0;
  94 my $reading              = 0;
  95 my $i                    = 0;
  96 my @lines                = ();
  97 my $larger               = 1;
  98 my $complete             = 0;
  99 my $total                = 0;
 100
 101 if ( @ARGV != 3 ) {
 102     &errorInCommandLine();
 103     exit ( -1 );
 104 }
 105
 106 $threshold  = $ARGV[ 0 ];
 107 $infile     = $ARGV[ 1 ];
 108 $outfile    = $ARGV[ 2 ];
 109 $summary_outfile = $outfile.".short";
 110
 111 if ( -e $outfile ) {
 112     die "\n$0: <<$outfile>> already exists.\n";
 113 }
 114 unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
 115     die "\n$0: <<$infile>> does not exist, is empty, or is not a plain textfile.\n";
 116 }
 117
 118
 119 if ( $threshold =~ /l/ ) {
 120     $larger = 0;
 121     $threshold =~ s/l//;
 122 }
 123 if ( $threshold =~ /c/ ) {
 124     $complete = 1;
 125     $threshold =~ s/c//;
 126 }
 127
 128 open( IN, "$infile" ) || die "\n$0: Cannot open file <<$infile>>: $!\n";
 129
 130 open( OUT, ">$outfile" ) || die "\n$0: Cannot create file \"$outfile\": $!\n";
 131 open( OUT_SUMMARY, ">$summary_outfile" ) || die "\n$0: Cannot create file \"$summary_outfile\": $!\n";
 132
 133 print OUT "bootstrapSelector.pl version: $VERSION\n\n";
 134 print OUT "Selection of RIO analysis results with top ortholgy\n";
 135 print OUT "bootstrap values greater or less than a threshold.\n\n";
 136 if ( $larger == 1 ) {
 137     print OUT "Threshold  : Grater than or equal to $threshold\n";
 138 }
 139 else {
 140     print OUT "Threshold  : Less than or equal to $threshold\n";
 141 }
 142 print OUT "In case of multiple copies of the same domain in the query:\n";
 143 if ( $complete == 1 ) {
 144     print OUT "All hits must meet threshold.\n";
 145 }
 146 else {
 147     print OUT "At least one hit must meet threshold.\n";
 148 }
 149 print OUT "Input file       : $infile\n";
 150 print OUT "Output file      : $outfile\n";
 151 print OUT "Output file short: $summary_outfile\n";
 152 print OUT "Date             : ".`date`."\n\n\n";
 153
 154 while ( $return_line = <IN> ) {
 155
 156     if ( $return_line =~ /^\s*# Annotation:\s*(.+)/ ) {
 157         $identifier = $1;
 158         $identifier = substr( $identifier, 0, 60);
 159         $analysis_performed = 0;
 160         $reading = 1;
 161         $i = 0;
 162         @lines = ();
 163     }
 164
 165     if ( $reading == 1 && $return_line =~ /^\s*RIO/ ) {
 166         $analysis_performed = 1;
 167     }
 168
 169     if ( $reading == 1
 170     && $return_line =~ /^\s*# ####################################/ ) {
 171         if ( $analysis_performed == 1 ) {
 172             &analyze();
 173         }
 174         $reading = 0;
 175     }
 176
 177     if ( $reading == 1 ) {
 178         $lines[ $i++ ] = $return_line;
 179     }
 180 }
 181
 182 close( IN );
 183
 184 print OUT "\n\nTotal: $total\n";
 185
 186 close( OUT );
 187 close( OUT_SUMMARY );
 188
 189 print "\nTotal: $total\n";
 190 print "Done.\n\n";
 191
 192 exit( 0 );
 193
 194
 195 sub analyze {
 196     my $j            = 0;
 197     my $results      = 0;
 198     my $o_bootstraps = 0;
 199     $top1            = "";
 200
 201     for ( $j = 0; $j < $i; $j++ ) {
 202
 203         if ( $lines[ $j ] =~ /^\s*--------\s+/ ) {
 204             $results = 1;
 205         }
 206         elsif ( $lines[ $j ] =~ /^\s*Distance\s+values\s+/i ) {
 207             $results = 0;
 208         }
 209         elsif ( $results == 1
 210         && ( $lines[ $j ] =~ /\S+\s+\S+\s+\S+\s*$/
 211         || $lines[ $j ] =~ /^\s*!NO\s+ORTHOLOGS/ ) ) {
 212
 213             if ( $lines[ $j ] =~ /^\s*!NO\s+ORTHOLOGS/ ) {
 214                 $o_bootstraps = 0;
 215             }
 216             else {
 217                 $lines[ $j ] =~ /(\S+)\s+\S+\s+\S+\s*$/;
 218                 $o_bootstraps = $1;
 219                 if ( $top1 eq "" ) {
 220                     $top1 = $lines[ $j ];
 221                     $top1 =~ s/\n//;
 222                     $top1 =~ s/\s{2,}/ /g;
 223                 }
 224             }
 225
 226             $results = 0;
 227
 228             if ( $o_bootstraps > 100 || $o_bootstraps < 0 ) {
 229                 print "o bootstraps: $o_bootstraps\n";
 230                 die "\n\n$0: Error: Boostrap value(s) out of range.\n\n";
 231             }
 232
 233             if ( $larger == 1 ) {
 234                 if ( $complete != 1 && $o_bootstraps >= $threshold ) {
 235                     &writeout();
 236                     $total++;
 237                     return;
 238                 }
 239                 elsif ( $complete == 1 && $o_bootstraps < $threshold ) {
 240                     return;
 241                 }
 242             }
 243             else {
 244                 if ( $complete != 1 && $o_bootstraps <= $threshold ) {
 245                     &writeout();
 246                     $total++;
 247                     return;
 248                 }
 249                 elsif ( $complete == 1 && $o_bootstraps > $threshold ) {
 250                     return;
 251                 }
 252             }
 253         }
 254     }
 255     if ( $complete == 1 ) {
 256         &writeout();
 257         $total++;
 258     }
 259     return;
 260 }
 261
 262
 263
 264 sub writeout {
 265     my $j = 0;
 266     print OUT "# ############################################################################\n";
 267     for ( $j = 0; $j < $i; ++$j ) {
 268         print OUT "$lines[ $j ]";
 269     }
 270     print OUT "# ############################################################################\n\n\n";
 271     print OUT_SUMMARY "$identifier [top 1: $top1]\n\n";
 272 }
 273
 274
 275
 276 sub errorInCommandLine {
 277     print "\n";
 278     print " bootstrapCounter.pl version: $VERSION\n";
 279     print " Usage: \"bootstrapSelector.pl <threshold options> <infile = Xrio.pl-output> <outfile>\"\n";
 280     print " Options: \"l\" for \"less or equal\" (\"grater or equal\" is default)\n";
 281     print "          \"c\" for \"all hits must meet threshold in case of\n";
 282     print "          multiple copies of the same domain in the query\"\n";
 283     print "          (default: \"at least one\")\n";
 284     print " Example:\n";
 285     print " \"bootstrapSelector.pl 95lc OUTFILE_At_1 At_1_out\"\n\n";
 286     print " Important: The result of this is meaningful ONLY if the thresholds\n";
 287     print " for output of the RIO analysis are set to zero (L=0 R=0).\n\n";
 288     exit( -1 );
 289 }
 290
 291