forester/java/src/org/forester/application/gsdi.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: www.phylosoft.org/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.Date;
  32 import java.util.List;
  33
  34
  35 import org.forester.io.parsers.PhylogenyParser;
  36 import org.forester.io.parsers.nhx.NHXParser;
  37 import org.forester.io.parsers.phyloxml.PhyloXmlParser;
  38 import org.forester.io.parsers.util.ParserUtils;
  39 import org.forester.io.writers.PhylogenyWriter;
  40 import org.forester.phylogeny.Phylogeny;
  41 import org.forester.phylogeny.PhylogenyMethods;
  42 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
  43 import org.forester.phylogeny.factories.PhylogenyFactory;
  44 import org.forester.sdi.GSDI;
  45 import org.forester.sdi.SDI;
  46 import org.forester.sdi.SDIse;
  47 import org.forester.util.CommandLineArguments;
  48 import org.forester.util.ForesterUtil;
  49
  50 public final class gsdi {
  51
  52     final static public boolean REPLACE_UNDERSCORES_IN_NH_SPECIES_TREE = true;
  53     final static private String STRIP_OPTION                 = "s";
  54     final static private String SDI_OPTION                   = "b";
  55     final static private String MOST_PARSIMONIOUS_OPTION     = "m";
  56     final static private String GUESS_FORMAT_OF_SPECIES_TREE = "q";
  57     final static private String HELP_OPTION_1                = "help";
  58     final static private String HELP_OPTION_2                = "h";
  59     final static private String DEFAULT_OUTFILE              = "gsdi_out.phylo.xml";
  60     final static private String PRG_NAME                     = "gsdi";
  61     final static private String PRG_VERSION                  = "0.5";
  62     final static private String PRG_DATE                     = "120608";
  63     final static private String PRG_DESC                     = "general speciation duplication inference";
  64     final static private String E_MAIL                       = "phylosoft@gmail.com";
  65     final static private String WWW                          = "www.phylosoft.org/forester";
  66
  67     public static void main( final String args[] ) {
  68         ForesterUtil.printProgramInformation( PRG_NAME,
  69                                               PRG_DESC,
  70                                               PRG_VERSION,
  71                                               PRG_DATE,
  72                                               E_MAIL,
  73                                               WWW,
  74                                               ForesterUtil.getForesterLibraryInformation() );
  75         CommandLineArguments cla = null;
  76         try {
  77             cla = new CommandLineArguments( args );
  78         }
  79         catch ( final Exception e ) {
  80             ForesterUtil.fatalError( PRG_NAME, e.getMessage() );
  81         }
  82         if ( cla.isOptionSet( gsdi.HELP_OPTION_1 ) || cla.isOptionSet( gsdi.HELP_OPTION_2 ) ) {
  83             System.out.println();
  84             gsdi.print_help();
  85             System.exit( 0 );
  86         }
  87         else if ( ( args.length < 2 ) || ( cla.getNumberOfNames() < 2 ) || ( cla.getNumberOfNames() > 3 ) ) {
  88             System.out.println();
  89             System.out.println( "Wrong number of arguments." );
  90             System.out.println();
  91             gsdi.print_help();
  92             System.exit( -1 );
  93         }
  94         final List<String> allowed_options = new ArrayList<String>();
  95         allowed_options.add( gsdi.STRIP_OPTION );
  96         allowed_options.add( gsdi.SDI_OPTION );
  97         allowed_options.add( gsdi.GUESS_FORMAT_OF_SPECIES_TREE );
  98         allowed_options.add( gsdi.MOST_PARSIMONIOUS_OPTION );
  99         final String dissallowed_options = cla.validateAllowedOptionsAsString( allowed_options );
 100         if ( dissallowed_options.length() > 0 ) {
 101             ForesterUtil.fatalError( gsdi.PRG_NAME, "unknown option(s): " + dissallowed_options );
 102         }
 103         boolean use_sdise = false;
 104         boolean strip = false;
 105         boolean most_parsimonous_duplication_model = false;
 106         boolean species_tree_in_phyloxml = true;
 107         if ( cla.isOptionSet( gsdi.STRIP_OPTION ) ) {
 108             strip = true;
 109         }
 110         if ( cla.isOptionSet( gsdi.SDI_OPTION ) ) {
 111             use_sdise = true;
 112         }
 113         if ( cla.isOptionSet( gsdi.MOST_PARSIMONIOUS_OPTION ) ) {
 114             if ( use_sdise ) {
 115                 ForesterUtil.fatalError( gsdi.PRG_NAME, "Can only use most parsimonious duplication mode with GSDI" );
 116             }
 117             most_parsimonous_duplication_model = true;
 118         }
 119         if ( cla.isOptionSet( gsdi.GUESS_FORMAT_OF_SPECIES_TREE ) ) {
 120             species_tree_in_phyloxml = false;
 121         }
 122         Phylogeny species_tree = null;
 123         Phylogeny gene_tree = null;
 124         File gene_tree_file = null;
 125         File species_tree_file = null;
 126         File out_file = null;
 127         try {
 128             gene_tree_file = cla.getFile( 0 );
 129             species_tree_file = cla.getFile( 1 );
 130             if ( cla.getNumberOfNames() == 3 ) {
 131                 out_file = cla.getFile( 2 );
 132             }
 133             else {
 134                 out_file = new File( gsdi.DEFAULT_OUTFILE );
 135             }
 136         }
 137         catch ( final IllegalArgumentException e ) {
 138             ForesterUtil.fatalError( gsdi.PRG_NAME, "error in command line: " + e.getMessage() );
 139         }
 140         if ( ForesterUtil.isReadableFile( gene_tree_file ) != null ) {
 141             ForesterUtil.fatalError( gsdi.PRG_NAME, ForesterUtil.isReadableFile( gene_tree_file ) );
 142         }
 143         if ( ForesterUtil.isReadableFile( species_tree_file ) != null ) {
 144             ForesterUtil.fatalError( gsdi.PRG_NAME, ForesterUtil.isReadableFile( species_tree_file ) );
 145         }
 146         if ( ForesterUtil.isWritableFile( out_file ) != null ) {
 147             ForesterUtil.fatalError( gsdi.PRG_NAME, ForesterUtil.isWritableFile( out_file ) );
 148         }
 149         try {
 150             final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 151             if ( species_tree_in_phyloxml ) {
 152                 species_tree = factory.create( species_tree_file, new PhyloXmlParser() )[ 0 ];
 153             }
 154             else {
 155                 final PhylogenyParser p = ParserUtils.createParserDependingOnFileType( species_tree_file, true );
 156                 if ( REPLACE_UNDERSCORES_IN_NH_SPECIES_TREE && p instanceof NHXParser ) {
 157                     (( NHXParser ) p ).setReplaceUnderscores( true );
 158                 }
 159                 species_tree = factory.create( species_tree_file, p )[ 0 ];
 160                 PhylogenyMethods.transferNodeNameToField( species_tree,
 161                                                           PhylogenyMethods.PhylogenyNodeField.TAXONOMY_SCIENTIFIC_NAME );
 162
 163             }
 164         }
 165         catch ( final IOException e ) {
 166             ForesterUtil.fatalError( gsdi.PRG_NAME,
 167                                      "Failed to read species tree from [" + gene_tree_file + "]: " + e.getMessage() );
 168         }
 169
 170         try {
 171             final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 172             gene_tree = factory.create( gene_tree_file, new PhyloXmlParser() )[ 0 ];
 173         }
 174         catch ( final IOException e ) {
 175             ForesterUtil.fatalError( gsdi.PRG_NAME,
 176                                      "Failed to read gene tree from [" + gene_tree_file + "]: " + e.getMessage() );
 177         }
 178         gene_tree.setRooted( true );
 179         species_tree.setRooted( true );
 180         if ( !gene_tree.isCompletelyBinary() ) {
 181             ForesterUtil.fatalError( gsdi.PRG_NAME, "gene tree (\"" + gene_tree_file + "\") is not completely binary" );
 182         }
 183         if ( use_sdise ) {
 184             if ( !species_tree.isCompletelyBinary() ) {
 185                 ForesterUtil.fatalError( gsdi.PRG_NAME, "species tree (\"" + species_tree_file
 186                         + "\") is not completely binary" );
 187             }
 188         }
 189         // For timing.
 190         // gene_tree = Helper.createBalancedTree( 10 );
 191         // species_tree = Helper.createBalancedTree( 13 );
 192         // species_tree = Helper.createUnbalancedTree( 1024 );
 193         // gene_tree = Helper.createUnbalancedTree( 8192 );
 194         // species_tree = gene_tree.copyTree();
 195         // gene_tree = species_tree.copyTree();
 196         // Helper.numberSpeciesInOrder( species_tree );
 197         // Helper.numberSpeciesInOrder( gene_tree );
 198         // Helper.randomizeSpecies( 1, 8192, gene_tree );
 199         // Helper.intervalNumberSpecies( gene_tree, 4096 );
 200         // Helper.numberSpeciesInDescOrder( gene_tree );
 201         System.out.println();
 202         System.out.println( "Strip species tree: " + strip );
 203         SDI sdi = null;
 204         final long start_time = new Date().getTime();
 205         try {
 206             if ( use_sdise ) {
 207                 System.out.println();
 208                 System.out.println( "Using SDIse algorithm" );
 209                 sdi = new SDIse( gene_tree, species_tree );
 210             }
 211             else {
 212                 System.out.println();
 213                 System.out.println( "Using GSDI algorithm" );
 214                 System.out.println();
 215                 System.out.println( "Use most parsimonous duplication model: " + most_parsimonous_duplication_model );
 216                 sdi = new GSDI( gene_tree, species_tree, most_parsimonous_duplication_model, true );
 217             }
 218         }
 219         catch ( final Exception e ) {
 220             ForesterUtil.fatalError( PRG_NAME, e.getLocalizedMessage() );
 221         }
 222         System.out.println();
 223         System.out.println( "Running time (excluding I/O): " + ( new Date().getTime() - start_time ) + "ms" );
 224         try {
 225             final PhylogenyWriter writer = new PhylogenyWriter();
 226             writer.toPhyloXML( out_file, gene_tree, 0 );
 227         }
 228         catch ( final IOException e ) {
 229             ForesterUtil.fatalError( PRG_NAME, "Failed to write to [" + out_file + "]: " + e.getMessage() );
 230         }
 231         System.out.println();
 232         System.out.println( "Successfully wrote resulting gene tree to: " + out_file );
 233         System.out.println();
 234         if ( use_sdise ) {
 235             sdi.computeMappingCostL();
 236             System.out.println( "Mapping cost                    : " + sdi.computeMappingCostL() );
 237         }
 238         System.out.println( "Number of duplications          : " + sdi.getDuplicationsSum() );
 239         if ( !use_sdise && !most_parsimonous_duplication_model ) {
 240             System.out.println( "Number of potential duplications: "
 241                     + ( ( GSDI ) sdi ).getSpeciationOrDuplicationEventsSum() );
 242         }
 243         if ( !use_sdise ) {
 244             System.out.println( "Number speciations              : " + ( ( GSDI ) sdi ).getSpeciationsSum() );
 245         }
 246         System.out.println();
 247     }
 248
 249     private static void print_help() {
 250         System.out.println( "Usage: " + gsdi.PRG_NAME
 251                 + " [-options] <gene tree in phyloXML format> <species tree in phyloXML format> [outfile]" );
 252         System.out.println();
 253         System.out.println( "Options:" );
 254         System.out.println( " -" + gsdi.STRIP_OPTION + ": to strip the species tree prior to duplication inference" );
 255         System.out.println( " -" + gsdi.SDI_OPTION + ": to use SDI algorithm instead of GSDI algorithm" );
 256         System.out.println( " -" + gsdi.MOST_PARSIMONIOUS_OPTION
 257                 + ": use most parimonious duplication model for GSDI: " );
 258         System.out.println( "     assign nodes as speciations which would otherwise be assiged" );
 259         System.out.println( "     as unknown because of polytomies in the species tree" );
 260         System.out.println( " -" + gsdi.GUESS_FORMAT_OF_SPECIES_TREE + ": to allow species tree in other formats than" );
 261         System.out.println( "     phyloXML (Newick, NHX, Nexus)" );
 262         System.out.println();
 263         System.out.println( "Species tree:" );
 264         System.out.println( " In phyloXML format (unless option -" + gsdi.GUESS_FORMAT_OF_SPECIES_TREE
 265                 + " is used), with taxonomy data in appropriate fields" );
 266         System.out.println();
 267         System.out.println( "Gene tree:" );
 268         System.out.println( " In phyloXM format, with taxonomy and sequence data in appropriate fields" );
 269         System.out.println();
 270         System.out.println( "Note -- GSDI algorithm is under development" );
 271         System.out.println();
 272     }
 273 }