forester/java/src/org/forester/application/count_support.java

   1 // $Id:
   2 // FORESTER -- software libraries and applications
   3 // for evolutionary biology research and applications.
   4 //
   5 // Copyright (C) 2008-2009 Christian M. Zmasek
   6 // Copyright (C) 2008-2009 Burnham Institute for Medical Research
   7 // All rights reserved
   8 //
   9 // This library is free software; you can redistribute it and/or
  10 // modify it under the terms of the GNU Lesser General Public
  11 // License as published by the Free Software Foundation; either
  12 // version 2.1 of the License, or (at your option) any later version.
  13 //
  14 // This library is distributed in the hope that it will be useful,
  15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  17 // Lesser General Public License for more details.
  18 //
  19 // You should have received a copy of the GNU Lesser General Public
  20 // License along with this library; if not, write to the Free Software
  21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  22 //
  23 // Contact: phylosoft @ gmail . com
  24 // WWW: https://sites.google.com/site/cmzmasek/home/software/forester
  25
  26 package org.forester.application;
  27
  28 import java.io.File;
  29 import java.io.IOException;
  30 import java.util.ArrayList;
  31 import java.util.Arrays;
  32 import java.util.List;
  33
  34 import org.forester.io.parsers.PhylogenyParser;
  35 import org.forester.io.parsers.util.ParserUtils;
  36 import org.forester.io.writers.PhylogenyWriter;
  37 import org.forester.phylogeny.Phylogeny;
  38 import org.forester.phylogeny.PhylogenyMethods;
  39 import org.forester.phylogeny.PhylogenyNode;
  40 import org.forester.phylogeny.data.Confidence;
  41 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
  42 import org.forester.phylogeny.factories.PhylogenyFactory;
  43 import org.forester.phylogeny.iterators.PhylogenyNodeIterator;
  44 import org.forester.tools.SupportCount;
  45 import org.forester.util.CommandLineArguments;
  46 import org.forester.util.ForesterUtil;
  47
  48 public class count_support {
  49
  50     final static private String  PRG_NAME                = "count_support";
  51     final static private String  PRG_VERSION             = "1.0";
  52     final static private String  PRG_DATE                = "2008.03.04";
  53     final static private String  E_MAIL                  = "phylosoft@gmail.com";
  54     final static private String  WWW                     = "www.phylosoft.org/forester/";
  55     private final static boolean WRITE_EVALUATORS_AS_NHX = false;
  56
  57     public static void main( final String args[] ) {
  58         ForesterUtil.printProgramInformation( PRG_NAME,
  59                                               null,
  60                                               PRG_VERSION,
  61                                               PRG_DATE,
  62                                               E_MAIL,
  63                                               WWW,
  64                                               ForesterUtil.getForesterLibraryInformation() );
  65         if ( ( args.length < 3 ) || ( args.length > 7 ) ) {
  66             System.out.println();
  67             System.out.println( count_support.PRG_NAME + ": wrong number of arguments" );
  68             System.out.println();
  69             System.out
  70             .println( "Usage: \"count_support [options] <file containing phylogeny to be evaluated> <file with phylogenies to be used for evaluation> <outfile> [outfile for evaluator phylogenies, "
  71                     + "always unstripped if -t=<d> option is used, otherwise strippedness is dependent on -s option]\"\n" );
  72             System.out
  73             .println( " Options: -s strip external nodes from evaluator phylogenies not found in phylogeny to be evaluated" );
  74             System.out.println( "        : -t=<d> threshold for similarity (0.0 to 1.0)" );
  75             System.out.println( "        : -n no branch lengths in outfile for evaluator phylogenies" );
  76             System.out.println();
  77             System.exit( -1 );
  78         }
  79         CommandLineArguments cla = null;
  80         try {
  81             cla = new CommandLineArguments( args );
  82         }
  83         catch ( final Exception e ) {
  84             ForesterUtil.fatalError( PRG_NAME, e.getMessage() );
  85         }
  86         final List<String> allowed_options = new ArrayList<String>();
  87         allowed_options.add( "s" );
  88         allowed_options.add( "t" );
  89         allowed_options.add( "n" );
  90         final String dissallowed_options = cla.validateAllowedOptionsAsString( allowed_options );
  91         if ( dissallowed_options.length() > 0 ) {
  92             ForesterUtil.fatalError( count_support.PRG_NAME, "Unknown option(s): " + dissallowed_options );
  93         }
  94         final File phylogeny_infile = cla.getFile( 0 );
  95         final File evaluators_infile = cla.getFile( 1 );
  96         final File phylogeny_outfile = cla.getFile( 2 );
  97         File evaluators_outfile = null;
  98         boolean branch_lengths_in_ev_out = true;
  99         if ( cla.isOptionSet( "n" ) ) {
 100             branch_lengths_in_ev_out = false;
 101         }
 102         if ( cla.getNumberOfNames() == 4 ) {
 103             evaluators_outfile = cla.getFile( 3 );
 104         }
 105         else {
 106             if ( !branch_lengths_in_ev_out ) {
 107                 ForesterUtil.fatalError( count_support.PRG_NAME,
 108                         "Cannot use -n option if no outfile for evaluators specified" );
 109             }
 110         }
 111         Phylogeny p = null;
 112         Phylogeny[] ev = null;
 113         try {
 114             final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 115             final PhylogenyParser pp = ParserUtils.createParserDependingOnFileType( phylogeny_infile, true );
 116             p = factory.create( phylogeny_infile, pp )[ 0 ];
 117         }
 118         catch ( final Exception e ) {
 119             ForesterUtil.fatalError( count_support.PRG_NAME,
 120                                      "Could not read \"" + phylogeny_infile + "\" [" + e.getMessage() + "]" );
 121         }
 122         try {
 123             final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 124             final PhylogenyParser pp = ParserUtils.createParserDependingOnFileType( evaluators_infile, true );
 125             ev = factory.create( evaluators_infile, pp );
 126         }
 127         catch ( final Exception e ) {
 128             ForesterUtil.fatalError( count_support.PRG_NAME,
 129                                      "Could not read \"" + evaluators_infile + "\" [" + e.getMessage() + "]" );
 130         }
 131         boolean strip = false;
 132         if ( cla.isOptionSet( "s" ) ) {
 133             strip = true;
 134         }
 135         double threshhold = -1.0;
 136         if ( cla.isOptionSet( "t" ) ) {
 137             try {
 138                 threshhold = cla.getOptionValueAsDouble( "t" );
 139             }
 140             catch ( final Exception e ) {
 141                 ForesterUtil.fatalError( count_support.PRG_NAME, "error in command line arguments: " + e.getMessage() );
 142             }
 143             if ( ( threshhold < 0 ) || ( threshhold > 1.0 ) ) {
 144                 ForesterUtil.fatalError( count_support.PRG_NAME,
 145                         "support threshold has to be between 0.0 and 1.0 (inclusive)" );
 146             }
 147         }
 148         List<Phylogeny> evaluator_phylogenies_above_threshold = null;
 149         try {
 150             if ( threshhold >= 0 ) {
 151                 evaluator_phylogenies_above_threshold = SupportCount.count( p, ev, strip, threshhold, true );
 152                 if ( evaluator_phylogenies_above_threshold.size() < 1 ) {
 153                     ForesterUtil.fatalError( "count_support", "appears like threshold for similarity is set too high" );
 154                 }
 155             }
 156             else {
 157                 SupportCount.count( p, ev, strip, true );
 158             }
 159         }
 160         catch ( final Exception e ) {
 161             ForesterUtil.fatalError( count_support.PRG_NAME, "Failure during support counting: " + e.getMessage() );
 162         }
 163         if ( threshhold >= 0 ) {
 164             count_support.normalizeSupport( p, 100, evaluator_phylogenies_above_threshold.size() );
 165             System.out.println( evaluator_phylogenies_above_threshold.size() + " out of " + ev.length
 166                                 + " evaluator phylogenies are above threshold of " + threshhold );
 167         }
 168         try {
 169             final PhylogenyWriter w = new PhylogenyWriter();
 170             w.toPhyloXML( phylogeny_outfile, p, 1 );
 171         }
 172         catch ( final IOException e ) {
 173             ForesterUtil.fatalError( count_support.PRG_NAME, "Failure to write output [" + e.getMessage() + "]" );
 174         }
 175         System.out.println();
 176         System.out.println( "Wrote phylogeny with support values to: " + phylogeny_outfile );
 177         if ( evaluators_outfile != null ) {
 178             try {
 179                 final PhylogenyWriter w = new PhylogenyWriter();
 180                 if ( evaluator_phylogenies_above_threshold != null ) {
 181                     System.out.println( "Writing " + evaluator_phylogenies_above_threshold.size()
 182                                         + " evaluator phylogenies above threshold of " + threshhold + " to: " + evaluators_outfile );
 183                     if ( count_support.WRITE_EVALUATORS_AS_NHX ) {
 184                         w.toNewHampshireX( evaluator_phylogenies_above_threshold, evaluators_outfile, ";"
 185                                 + ForesterUtil.getLineSeparator() );
 186                     }
 187                     else {
 188                         w.toNewHampshire( evaluator_phylogenies_above_threshold,
 189                                           branch_lengths_in_ev_out,
 190                                           evaluators_outfile,
 191                                           ";" + ForesterUtil.getLineSeparator() );
 192                     }
 193                 }
 194                 else {
 195                     System.out.println( "Writing " + ev.length + " evaluator phylogenies to :" + evaluators_outfile );
 196                     if ( count_support.WRITE_EVALUATORS_AS_NHX ) {
 197                         w.toNewHampshireX( Arrays.asList( ev ),
 198                                            evaluators_outfile,
 199                                            ";" + ForesterUtil.getLineSeparator() );
 200                     }
 201                     else {
 202                         w.toNewHampshire( Arrays.asList( ev ), branch_lengths_in_ev_out, evaluators_outfile, ";"
 203                                 + ForesterUtil.getLineSeparator() );
 204                     }
 205                 }
 206             }
 207             catch ( final IOException e ) {
 208                 ForesterUtil.fatalError( count_support.PRG_NAME, "Failure to write output [" + e.getMessage() + "]" );
 209             }
 210         }
 211         System.out.println();
 212         System.out.println( "Done." );
 213         System.out.println();
 214     }
 215
 216     private static void normalizeSupport( final Phylogeny p, final double normalized_max, final int number_phylos ) {
 217         double min = Double.MAX_VALUE;
 218         double max = -Double.MAX_VALUE;
 219         double sum = 0.0;
 220         int n = 0;
 221         for( final PhylogenyNodeIterator iter = p.iteratorPostorder(); iter.hasNext(); ) {
 222             final PhylogenyNode node = iter.next();
 223             if ( !node.isRoot() && !node.isExternal() ) {
 224                 final double b = PhylogenyMethods.getConfidenceValue( node );
 225                 if ( b > max ) {
 226                     max = b;
 227                 }
 228                 if ( ( b >= 0 ) && ( b < min ) ) {
 229                     min = b;
 230                 }
 231                 sum += b;
 232                 ++n;
 233             }
 234         }
 235         double av = sum / n;
 236         System.out.println( "Max support before normalization is    : " + max );
 237         System.out.println( "Min support before normalization is    : " + min );
 238         System.out.println( "Average support before normalization is: " + av + " (=" + sum + "/" + n + ")" );
 239         System.out.println( "Normalizing so that theoretical maximum support value is: " + normalized_max );
 240         System.out.println( "Number of phylogenies used in support analysis: " + number_phylos );
 241         final double f = normalized_max / number_phylos;
 242         min = Double.MAX_VALUE;
 243         max = -Double.MAX_VALUE;
 244         sum = 0.0;
 245         n = 0;
 246         for( final PhylogenyNodeIterator iter = p.iteratorPostorder(); iter.hasNext(); ) {
 247             final PhylogenyNode node = iter.next();
 248             if ( node.isRoot() || node.isExternal() ) {
 249                 PhylogenyMethods.setBootstrapConfidence( node, Confidence.CONFIDENCE_DEFAULT_VALUE );
 250             }
 251             else {
 252                 double b = PhylogenyMethods.getConfidenceValue( node );
 253                 b = f * b;
 254                 PhylogenyMethods.setBootstrapConfidence( node, b );
 255                 if ( b > max ) {
 256                     max = b;
 257                 }
 258                 if ( ( b >= 0 ) && ( b < min ) ) {
 259                     min = b;
 260                 }
 261                 sum += b;
 262                 ++n;
 263             }
 264         }
 265         av = sum / n;
 266         System.out.println( "Max support after normalization is    : " + max );
 267         System.out.println( "Min support after normalization is    : " + min );
 268         System.out.println( "Average support after normalization is: " + av + " (=" + sum + "/" + n + ")" );
 269     }
 270 }