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