forester/java/src/org/forester/rio/RIO.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 // Copyright (C) 2000-2001 Washington University School of Medicine
   8 // and Howard Hughes Medical Institute
   9 // All rights reserved
  10 //
  11 // This library is free software; you can redistribute it and/or
  12 // modify it under the terms of the GNU Lesser General Public
  13 // License as published by the Free Software Foundation; either
  14 // version 2.1 of the License, or (at your option) any later version.
  15 //
  16 // This library is distributed in the hope that it will be useful,
  17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  19 // Lesser General Public License for more details.
  20 //
  21 // You should have received a copy of the GNU Lesser General Public
  22 // License along with this library; if not, write to the Free Software
  23 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  24 //
  25 // Contact: phylosoft @ gmail . com
  26 // WWW: https://sites.google.com/site/cmzmasek/home/software/forester
  27
  28 package org.forester.rio;
  29
  30 import java.io.File;
  31 import java.io.FileNotFoundException;
  32 import java.io.IOException;
  33 import java.text.DecimalFormat;
  34 import java.util.ArrayList;
  35 import java.util.Collections;
  36 import java.util.HashMap;
  37 import java.util.HashSet;
  38 import java.util.List;
  39 import java.util.Set;
  40 import java.util.SortedSet;
  41 import java.util.TreeSet;
  42
  43 import org.forester.datastructures.IntMatrix;
  44 import org.forester.io.parsers.IteratingPhylogenyParser;
  45 import org.forester.io.parsers.PhylogenyParser;
  46 import org.forester.io.parsers.nexus.NexusPhylogeniesParser;
  47 import org.forester.io.parsers.nhx.NHXParser;
  48 import org.forester.io.parsers.nhx.NHXParser.TAXONOMY_EXTRACTION;
  49 import org.forester.io.parsers.util.ParserUtils;
  50 import org.forester.phylogeny.Phylogeny;
  51 import org.forester.phylogeny.PhylogenyMethods;
  52 import org.forester.phylogeny.PhylogenyNode;
  53 import org.forester.phylogeny.data.Taxonomy;
  54 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
  55 import org.forester.phylogeny.factories.PhylogenyFactory;
  56 import org.forester.sdi.GSDI;
  57 import org.forester.sdi.GSDIR;
  58 import org.forester.sdi.SDIException;
  59 import org.forester.sdi.SDIR;
  60 import org.forester.sdi.SDIutil;
  61 import org.forester.sdi.SDIutil.ALGORITHM;
  62 import org.forester.sdi.SDIutil.TaxonomyComparisonBase;
  63 import org.forester.util.BasicDescriptiveStatistics;
  64 import org.forester.util.ForesterUtil;
  65
  66 public final class RIO {
  67
  68     public static final int                  DEFAULT_RANGE = -1;
  69     private static final int                 END_OF_GT     = Integer.MAX_VALUE;
  70     private static IntMatrix                 _m;
  71     private Phylogeny[]                      _analyzed_gene_trees;
  72     private List<PhylogenyNode>              _removed_gene_tree_nodes;
  73     private int                              _ext_nodes;
  74     private int                              _int_nodes;
  75     private TaxonomyComparisonBase           _gsdir_tax_comp_base;
  76     private final StringBuilder              _log;
  77     private final BasicDescriptiveStatistics _duplications_stats;
  78     private final boolean                    _produce_log;
  79     private final boolean                    _verbose;
  80     private final REROOTING                  _rerooting;
  81     private final Phylogeny                  _species_tree;
  82     private Phylogeny                        _min_dub_gene_tree;
  83
  84     private RIO( final IteratingPhylogenyParser p,
  85                  final Phylogeny species_tree,
  86                  final ALGORITHM algorithm,
  87                  final REROOTING rerooting,
  88                  final String outgroup,
  89                  int first,
  90                  int last,
  91                  final boolean produce_log,
  92                  final boolean verbose ) throws IOException, SDIException, RIOException {
  93         if ( ( last == DEFAULT_RANGE ) && ( first >= 0 ) ) {
  94             last = END_OF_GT;
  95         }
  96         else if ( ( first == DEFAULT_RANGE ) && ( last >= 0 ) ) {
  97             first = 0;
  98         }
  99         removeSingleDescendentsNodes( species_tree, verbose );
 100         p.reset();
 101         checkPreconditions( p, species_tree, rerooting, outgroup, first, last );
 102         _produce_log = produce_log;
 103         _verbose = verbose;
 104         _rerooting = rerooting;
 105         _ext_nodes = -1;
 106         _int_nodes = -1;
 107         _log = new StringBuilder();
 108         _gsdir_tax_comp_base = null;
 109         _analyzed_gene_trees = null;
 110         _removed_gene_tree_nodes = null;
 111         _duplications_stats = new BasicDescriptiveStatistics();
 112         p.reset();
 113         inferOrthologs( p, species_tree, algorithm, outgroup, first, last );
 114         _species_tree = species_tree;
 115     }
 116
 117     private RIO( final Phylogeny[] gene_trees,
 118                  final Phylogeny species_tree,
 119                  final ALGORITHM algorithm,
 120                  final REROOTING rerooting,
 121                  final String outgroup,
 122                  int first,
 123                  int last,
 124                  final boolean produce_log,
 125                  final boolean verbose ) throws IOException, SDIException, RIOException {
 126         if ( ( last == DEFAULT_RANGE ) && ( first >= 0 ) ) {
 127             last = gene_trees.length - 1;
 128         }
 129         else if ( ( first == DEFAULT_RANGE ) && ( last >= 0 ) ) {
 130             first = 0;
 131         }
 132         removeSingleDescendentsNodes( species_tree, verbose );
 133         checkPreconditions( gene_trees, species_tree, rerooting, outgroup, first, last );
 134         _produce_log = produce_log;
 135         _verbose = verbose;
 136         _rerooting = rerooting;
 137         _ext_nodes = -1;
 138         _int_nodes = -1;
 139         _log = new StringBuilder();
 140         _gsdir_tax_comp_base = null;
 141         _analyzed_gene_trees = null;
 142         _removed_gene_tree_nodes = null;
 143         _duplications_stats = new BasicDescriptiveStatistics();
 144         inferOrthologs( gene_trees, species_tree, algorithm, outgroup, first, last );
 145         _species_tree = species_tree;
 146     }
 147
 148     public final Phylogeny[] getAnalyzedGeneTrees() {
 149         return _analyzed_gene_trees;
 150     }
 151
 152     public final BasicDescriptiveStatistics getDuplicationsStatistics() {
 153         return _duplications_stats;
 154     }
 155
 156     /**
 157      * Returns the numbers of number of ext nodes in gene trees analyzed (after
 158      * stripping).
 159      *
 160      * @return number of ext nodes in gene trees analyzed (after stripping)
 161      */
 162     public final int getExtNodesOfAnalyzedGeneTrees() {
 163         return _ext_nodes;
 164     }
 165
 166     public final TaxonomyComparisonBase getGSDIRtaxCompBase() {
 167         return _gsdir_tax_comp_base;
 168     }
 169
 170     /**
 171      * Returns the numbers of number of int nodes in gene trees analyzed (after
 172      * stripping).
 173      *
 174      * @return number of int nodes in gene trees analyzed (after stripping)
 175      */
 176     public final int getIntNodesOfAnalyzedGeneTrees() {
 177         return _int_nodes;
 178     }
 179
 180     public final StringBuilder getLog() {
 181         return _log;
 182     }
 183
 184     final public Phylogeny getMinDuplicationsGeneTree() {
 185         return _min_dub_gene_tree;
 186     }
 187
 188     public final IntMatrix getOrthologTable() {
 189         return _m;
 190     }
 191
 192     public final List<PhylogenyNode> getRemovedGeneTreeNodes() {
 193         return _removed_gene_tree_nodes;
 194     }
 195
 196     public final Phylogeny getSpeciesTree() {
 197         return _species_tree;
 198     }
 199
 200     private final void inferOrthologs( final IteratingPhylogenyParser parser,
 201                                        final Phylogeny species_tree,
 202                                        final ALGORITHM algorithm,
 203                                        final String outgroup,
 204                                        int first,
 205                                        final int last ) throws SDIException, RIOException, FileNotFoundException,
 206             IOException {
 207         if ( !parser.hasNext() ) {
 208             throw new RIOException( "no gene trees to analyze" );
 209         }
 210         if ( log() ) {
 211             preLog( -1, species_tree, algorithm, outgroup );
 212         }
 213         if ( _verbose ) {
 214             System.out.println();
 215         }
 216         final DecimalFormat pf = new java.text.DecimalFormat( "000" );
 217         int gene_tree_ext_nodes = 0;
 218         int i = 0;
 219         int counter = 0;
 220         final boolean no_range = ( first < 0 ) || ( last < first );
 221         while ( parser.hasNext() ) {
 222             final Phylogeny gt = parser.next();
 223             if ( no_range || ( ( i >= first ) && ( i <= last ) ) ) {
 224                 if ( _verbose ) {
 225                     ForesterUtil.updateProgress( i, pf );
 226                 }
 227                 if ( counter == 0 ) {
 228                     if ( algorithm == ALGORITHM.SDIR ) {
 229                         // Removes from species_tree all species not found in gene_tree.
 230                         PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( gt, species_tree );
 231                         if ( species_tree.isEmpty() ) {
 232                             throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 233                         }
 234                     }
 235                     gene_tree_ext_nodes = gt.getNumberOfExternalNodes();
 236                 }
 237                 else if ( gene_tree_ext_nodes != gt.getNumberOfExternalNodes() ) {
 238                     throw new RIOException( "gene tree #" + i + " has a different number of external nodes ("
 239                             + gt.getNumberOfExternalNodes() + ") than the preceding gene tree(s) ("
 240                             + gene_tree_ext_nodes + ")" );
 241                 }
 242                 if ( algorithm == ALGORITHM.SDIR ) {
 243                     // Removes from gene_tree all species not found in species_tree.
 244                     PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gt );
 245                     if ( gt.isEmpty() ) {
 246                         throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 247                     }
 248                 }
 249                 final Phylogeny analyzed_gt = performOrthologInference( gt, species_tree, algorithm, outgroup, counter );
 250                 RIO.calculateOrthologTable( analyzed_gt, true, counter );
 251                 ++counter;
 252             }
 253             ++i;
 254         }
 255         if ( ( first >= 0 ) && ( counter == 0 ) && ( i > 0 ) ) {
 256             throw new RIOException( "attempt to analyze first gene tree #" + first + " in a set of " + i );
 257         }
 258         if ( no_range ) {
 259             first = 0;
 260         }
 261         if ( log() ) {
 262             postLog( species_tree, first, first + counter - 1 );
 263         }
 264         if ( _verbose ) {
 265             System.out.println();
 266             System.out.println();
 267         }
 268     }
 269
 270     private final void inferOrthologs( final Phylogeny[] gene_trees,
 271                                        final Phylogeny species_tree,
 272                                        final ALGORITHM algorithm,
 273                                        final String outgroup,
 274                                        final int first,
 275                                        final int last ) throws SDIException, RIOException, FileNotFoundException,
 276             IOException {
 277         if ( algorithm == ALGORITHM.SDIR ) {
 278             // Removes from species_tree all species not found in gene_tree.
 279             PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( gene_trees[ 0 ], species_tree );
 280             if ( species_tree.isEmpty() ) {
 281                 throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 282             }
 283         }
 284         final Phylogeny[] my_gene_trees;
 285         if ( ( first >= 0 ) && ( last >= first ) && ( last < gene_trees.length ) ) {
 286             my_gene_trees = new Phylogeny[ ( 1 + last ) - first ];
 287             int c = 0;
 288             for( int i = first; i <= last; ++i ) {
 289                 my_gene_trees[ c++ ] = gene_trees[ i ];
 290             }
 291         }
 292         else {
 293             my_gene_trees = gene_trees;
 294         }
 295         if ( log() ) {
 296             preLog( gene_trees.length, species_tree, algorithm, outgroup );
 297         }
 298         if ( _verbose && ( my_gene_trees.length >= 4 ) ) {
 299             System.out.println();
 300         }
 301         _analyzed_gene_trees = new Phylogeny[ my_gene_trees.length ];
 302         int gene_tree_ext_nodes = 0;
 303         for( int i = 0; i < my_gene_trees.length; ++i ) {
 304             final Phylogeny gt = my_gene_trees[ i ];
 305             if ( _verbose && ( my_gene_trees.length > 4 ) ) {
 306                 ForesterUtil.updateProgress( ( ( double ) i ) / my_gene_trees.length );
 307             }
 308             if ( i == 0 ) {
 309                 gene_tree_ext_nodes = gt.getNumberOfExternalNodes();
 310             }
 311             else if ( gene_tree_ext_nodes != gt.getNumberOfExternalNodes() ) {
 312                 throw new RIOException( "gene tree #" + i + " has a different number of external nodes ("
 313                         + gt.getNumberOfExternalNodes() + ") than the preceding gene tree(s) (" + gene_tree_ext_nodes
 314                         + ")" );
 315             }
 316             if ( algorithm == ALGORITHM.SDIR ) {
 317                 // Removes from gene_tree all species not found in species_tree.
 318                 PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gt );
 319                 if ( gt.isEmpty() ) {
 320                     throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 321                 }
 322             }
 323             _analyzed_gene_trees[ i ] = performOrthologInference( gt, species_tree, algorithm, outgroup, i );
 324         }
 325         if ( log() ) {
 326             postLog( species_tree, first, last );
 327         }
 328         if ( _verbose && ( my_gene_trees.length > 4 ) ) {
 329             System.out.println();
 330             System.out.println();
 331         }
 332     }
 333
 334     private final boolean log() {
 335         return _produce_log;
 336     }
 337
 338     private final void log( final String s ) {
 339         _log.append( s );
 340         _log.append( ForesterUtil.LINE_SEPARATOR );
 341     }
 342
 343     private final void logRemovedGeneTreeNodes() {
 344         log( "Species stripped from gene trees:" );
 345         final SortedSet<String> rn = new TreeSet<String>();
 346         for( final PhylogenyNode n : getRemovedGeneTreeNodes() ) {
 347             final Taxonomy t = n.getNodeData().getTaxonomy();
 348             switch ( getGSDIRtaxCompBase() ) {
 349                 case CODE: {
 350                     rn.add( t.getTaxonomyCode() );
 351                     break;
 352                 }
 353                 case ID: {
 354                     rn.add( t.getIdentifier().toString() );
 355                     break;
 356                 }
 357                 case SCIENTIFIC_NAME: {
 358                     rn.add( t.getScientificName() );
 359                     break;
 360                 }
 361             }
 362         }
 363         for( final String s : rn ) {
 364             log( s );
 365         }
 366         log( "" );
 367     }
 368
 369     private final Phylogeny performOrthologInference( final Phylogeny gene_tree,
 370                                                       final Phylogeny species_tree,
 371                                                       final ALGORITHM algorithm,
 372                                                       final String outgroup,
 373                                                       final int i ) throws SDIException, RIOException {
 374         final Phylogeny assigned_tree;
 375         switch ( algorithm ) {
 376             case SDIR: {
 377                 assigned_tree = performOrthologInferenceBySDI( gene_tree, species_tree );
 378                 break;
 379             }
 380             case GSDIR: {
 381                 assigned_tree = performOrthologInferenceByGSDI( gene_tree, species_tree, outgroup, i );
 382                 break;
 383             }
 384             default: {
 385                 throw new IllegalArgumentException( "illegal algorithm: " + algorithm );
 386             }
 387         }
 388         if ( i == 0 ) {
 389             _ext_nodes = assigned_tree.getNumberOfExternalNodes();
 390             _int_nodes = assigned_tree.getNumberOfInternalNodes();
 391         }
 392         else if ( _ext_nodes != assigned_tree.getNumberOfExternalNodes() ) {
 393             throw new RIOException( "after stripping gene tree #" + i + " has a different number of external nodes ("
 394                     + assigned_tree.getNumberOfExternalNodes() + ") than the preceding gene tree(s) (" + _ext_nodes
 395                     + ")" );
 396         }
 397         return assigned_tree;
 398     }
 399
 400     private final Phylogeny performOrthologInferenceByGSDI( final Phylogeny gene_tree,
 401                                                             final Phylogeny species_tree,
 402                                                             final String outgroup,
 403                                                             final int i ) throws SDIException, RIOException {
 404         final Phylogeny assigned_tree;
 405         final int dups;
 406         if ( _rerooting == REROOTING.BY_ALGORITHM ) {
 407             final GSDIR gsdir = new GSDIR( gene_tree, species_tree, true, i == 0 );
 408             assigned_tree = gsdir.getMinDuplicationsSumGeneTree();
 409             if ( i == 0 ) {
 410                 _removed_gene_tree_nodes = gsdir.getStrippedExternalGeneTreeNodes();
 411                 for( final PhylogenyNode r : _removed_gene_tree_nodes ) {
 412                     if ( !r.getNodeData().isHasTaxonomy() ) {
 413                         throw new RIOException( "node with no (appropriate) taxonomic information found in gene tree #"
 414                                 + i + ": " + r.toString() );
 415                     }
 416                 }
 417             }
 418             if ( i == 0 ) {
 419                 _gsdir_tax_comp_base = gsdir.getTaxCompBase();
 420             }
 421             dups = gsdir.getMinDuplicationsSum();
 422         }
 423         else {
 424             if ( _rerooting == REROOTING.MIDPOINT ) {
 425                 PhylogenyMethods.midpointRoot( gene_tree );
 426             }
 427             else if ( _rerooting == REROOTING.OUTGROUP ) {
 428                 final PhylogenyNode n = gene_tree.getNode( outgroup );
 429                 gene_tree.reRoot( n );
 430             }
 431             final GSDI gsdi = new GSDI( gene_tree, species_tree, true, true, true );
 432             _removed_gene_tree_nodes = gsdi.getStrippedExternalGeneTreeNodes();
 433             for( final PhylogenyNode r : _removed_gene_tree_nodes ) {
 434                 if ( !r.getNodeData().isHasTaxonomy() ) {
 435                     throw new RIOException( "node with no (appropriate) taxonomic information found in gene tree #" + i
 436                             + ": " + r.toString() );
 437                 }
 438             }
 439             assigned_tree = gene_tree;
 440             if ( i == 0 ) {
 441                 _gsdir_tax_comp_base = gsdi.getTaxCompBase();
 442             }
 443             dups = gsdi.getDuplicationsSum();
 444         }
 445         if ( ( i == 0 ) || ( dups < _duplications_stats.getMin() ) ) {
 446             _min_dub_gene_tree = assigned_tree;
 447         }
 448         _duplications_stats.addValue( dups );
 449         return assigned_tree;
 450     }
 451
 452     private final Phylogeny performOrthologInferenceBySDI( final Phylogeny gene_tree, final Phylogeny species_tree )
 453             throws SDIException {
 454         final SDIR sdir = new SDIR();
 455         return sdir.infer( gene_tree, species_tree, false, true, true, true, 1 )[ 0 ];
 456     }
 457
 458     private final void postLog( final Phylogeny species_tree, final int first, final int last ) {
 459         log( "" );
 460         if ( ( getRemovedGeneTreeNodes() != null ) && ( getRemovedGeneTreeNodes().size() > 0 ) ) {
 461             logRemovedGeneTreeNodes();
 462         }
 463         log( "Species tree external nodes (after stripping)   : " + species_tree.getNumberOfExternalNodes() );
 464         log( "Species tree polytomies (after stripping)       : "
 465                 + PhylogenyMethods.countNumberOfPolytomies( species_tree ) );
 466         log( "Taxonomy linking based on                       : " + getGSDIRtaxCompBase() );
 467         final java.text.DecimalFormat df = new java.text.DecimalFormat( "0.#" );
 468         if ( ( first >= 0 ) && ( last >= 0 ) ) {
 469             log( "Gene trees analyzed range                       : " + first + "-" + last );
 470         }
 471         log( "Gene trees analyzed                             : " + _duplications_stats.getN() );
 472         log( "Mean number of duplications                     : " + df.format( _duplications_stats.arithmeticMean() )
 473                 + " (sd: " + df.format( _duplications_stats.sampleStandardDeviation() ) + ")" + " ("
 474                 + df.format( ( 100.0 * _duplications_stats.arithmeticMean() ) / getIntNodesOfAnalyzedGeneTrees() )
 475                 + "%)" );
 476         if ( _duplications_stats.getN() > 3 ) {
 477             log( "Median number of duplications                   : " + df.format( _duplications_stats.median() )
 478                     + " (" + df.format( ( 100.0 * _duplications_stats.median() ) / getIntNodesOfAnalyzedGeneTrees() )
 479                     + "%)" );
 480         }
 481         log( "Minimum duplications                            : " + ( int ) _duplications_stats.getMin() + " ("
 482                 + df.format( ( 100.0 * _duplications_stats.getMin() ) / getIntNodesOfAnalyzedGeneTrees() ) + "%)" );
 483         log( "Maximum duplications                            : " + ( int ) _duplications_stats.getMax() + " ("
 484                 + df.format( ( 100.0 * _duplications_stats.getMax() ) / getIntNodesOfAnalyzedGeneTrees() ) + "%)" );
 485         log( "Gene tree internal nodes                        : " + getIntNodesOfAnalyzedGeneTrees() );
 486         log( "Gene tree external nodes                        : " + getExtNodesOfAnalyzedGeneTrees() );
 487     }
 488
 489     private final void preLog( final int gene_trees,
 490                                final Phylogeny species_tree,
 491                                final ALGORITHM algorithm,
 492                                final String outgroup ) {
 493         if ( gene_trees > 0 ) {
 494             log( "Number of gene trees (total)                    : " + gene_trees );
 495         }
 496         log( "Algorithm                                       : " + algorithm );
 497         log( "Species tree external nodes (prior to stripping): " + species_tree.getNumberOfExternalNodes() );
 498         log( "Species tree polytomies (prior to stripping)    : "
 499                 + PhylogenyMethods.countNumberOfPolytomies( species_tree ) );
 500         String rs = "";
 501         switch ( _rerooting ) {
 502             case BY_ALGORITHM: {
 503                 rs = "minimizing duplications";
 504                 break;
 505             }
 506             case MIDPOINT: {
 507                 rs = "midpoint";
 508                 break;
 509             }
 510             case OUTGROUP: {
 511                 rs = "outgroup: " + outgroup;
 512                 break;
 513             }
 514             case NONE: {
 515                 rs = "none";
 516                 break;
 517             }
 518         }
 519         log( "Re-rooting                                      : " + rs );
 520     }
 521
 522     public final static IntMatrix calculateOrthologTable( final Phylogeny[] analyzed_gene_trees, final boolean sort )
 523             throws RIOException {
 524         final List<String> labels = new ArrayList<String>();
 525         final Set<String> labels_set = new HashSet<String>();
 526         for( final PhylogenyNode n : analyzed_gene_trees[ 0 ].getExternalNodes() ) {
 527             final String label = obtainLabel( labels_set, n );
 528             labels_set.add( label );
 529             labels.add( label );
 530         }
 531         if ( sort ) {
 532             Collections.sort( labels );
 533         }
 534         final IntMatrix m = new IntMatrix( labels );
 535         int counter = 0;
 536         for( final Phylogeny gt : analyzed_gene_trees ) {
 537             counter++;
 538             updateCounts( m, counter, gt );
 539         }
 540         return m;
 541     }
 542
 543     public final static RIO executeAnalysis( final File gene_trees_file,
 544                                              final File species_tree_file,
 545                                              final ALGORITHM algorithm,
 546                                              final REROOTING rerooting,
 547                                              final String outgroup,
 548                                              final int first,
 549                                              final int last,
 550                                              final boolean produce_log,
 551                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 552         final Phylogeny[] gene_trees = parseGeneTrees( gene_trees_file );
 553         if ( gene_trees.length < 1 ) {
 554             throw new RIOException( "\"" + gene_trees_file + "\" is devoid of appropriate gene trees" );
 555         }
 556         final Phylogeny species_tree = SDIutil.parseSpeciesTree( gene_trees[ 0 ],
 557                                                                  species_tree_file,
 558                                                                  false,
 559                                                                  true,
 560                                                                  TAXONOMY_EXTRACTION.NO );
 561         return new RIO( gene_trees, species_tree, algorithm, rerooting, outgroup, first, last, produce_log, verbose );
 562     }
 563
 564     public final static RIO executeAnalysis( final File gene_trees_file,
 565                                              final Phylogeny species_tree,
 566                                              final ALGORITHM algorithm,
 567                                              final REROOTING rerooting,
 568                                              final String outgroup,
 569                                              final boolean produce_log,
 570                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 571         return new RIO( parseGeneTrees( gene_trees_file ),
 572                         species_tree,
 573                         algorithm,
 574                         rerooting,
 575                         outgroup,
 576                         DEFAULT_RANGE,
 577                         DEFAULT_RANGE,
 578                         produce_log,
 579                         verbose );
 580     }
 581
 582     public final static RIO executeAnalysis( final File gene_trees_file,
 583                                              final Phylogeny species_tree,
 584                                              final ALGORITHM algorithm,
 585                                              final REROOTING rerooting,
 586                                              final String outgroup,
 587                                              final int first,
 588                                              final int last,
 589                                              final boolean produce_log,
 590                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 591         return new RIO( parseGeneTrees( gene_trees_file ),
 592                         species_tree,
 593                         algorithm,
 594                         rerooting,
 595                         outgroup,
 596                         first,
 597                         last,
 598                         produce_log,
 599                         verbose );
 600     }
 601
 602     public final static RIO executeAnalysis( final IteratingPhylogenyParser p,
 603                                              final File species_tree_file,
 604                                              final ALGORITHM algorithm,
 605                                              final REROOTING rerooting,
 606                                              final String outgroup,
 607                                              final int first,
 608                                              final int last,
 609                                              final boolean produce_log,
 610                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 611         final Phylogeny g0 = p.next();
 612         if ( ( g0 == null ) || g0.isEmpty() || ( g0.getNumberOfExternalNodes() < 2 ) ) {
 613             throw new RIOException( "input file does not seem to contain any gene trees" );
 614         }
 615         final Phylogeny species_tree = SDIutil.parseSpeciesTree( g0,
 616                                                                  species_tree_file,
 617                                                                  false,
 618                                                                  true,
 619                                                                  TAXONOMY_EXTRACTION.NO );
 620         p.reset();
 621         return new RIO( p, species_tree, algorithm, rerooting, outgroup, first, last, produce_log, verbose );
 622     }
 623
 624     public final static RIO executeAnalysis( final IteratingPhylogenyParser p,
 625                                              final Phylogeny species_tree,
 626                                              final ALGORITHM algorithm,
 627                                              final REROOTING rerooting,
 628                                              final String outgroup,
 629                                              final boolean produce_log,
 630                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 631         return new RIO( p,
 632                         species_tree,
 633                         algorithm,
 634                         rerooting,
 635                         outgroup,
 636                         DEFAULT_RANGE,
 637                         DEFAULT_RANGE,
 638                         produce_log,
 639                         verbose );
 640     }
 641
 642     public final static RIO executeAnalysis( final IteratingPhylogenyParser p,
 643                                              final Phylogeny species_tree,
 644                                              final ALGORITHM algorithm,
 645                                              final REROOTING rerooting,
 646                                              final String outgroup,
 647                                              final int first,
 648                                              final int last,
 649                                              final boolean produce_log,
 650                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 651         return new RIO( p, species_tree, algorithm, rerooting, outgroup, first, last, produce_log, verbose );
 652     }
 653
 654     public final static RIO executeAnalysis( final Phylogeny[] gene_trees, final Phylogeny species_tree )
 655             throws IOException, SDIException, RIOException {
 656         return new RIO( gene_trees,
 657                         species_tree,
 658                         ALGORITHM.GSDIR,
 659                         REROOTING.BY_ALGORITHM,
 660                         null,
 661                         DEFAULT_RANGE,
 662                         DEFAULT_RANGE,
 663                         false,
 664                         false );
 665     }
 666
 667     public final static RIO executeAnalysis( final Phylogeny[] gene_trees,
 668                                              final Phylogeny species_tree,
 669                                              final ALGORITHM algorithm,
 670                                              final REROOTING rerooting,
 671                                              final String outgroup,
 672                                              final boolean produce_log,
 673                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 674         return new RIO( gene_trees,
 675                         species_tree,
 676                         algorithm,
 677                         rerooting,
 678                         outgroup,
 679                         DEFAULT_RANGE,
 680                         DEFAULT_RANGE,
 681                         produce_log,
 682                         verbose );
 683     }
 684
 685     public final static RIO executeAnalysis( final Phylogeny[] gene_trees,
 686                                              final Phylogeny species_tree,
 687                                              final ALGORITHM algorithm,
 688                                              final REROOTING rerooting,
 689                                              final String outgroup,
 690                                              final int first,
 691                                              final int last,
 692                                              final boolean produce_log,
 693                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 694         return new RIO( gene_trees, species_tree, algorithm, rerooting, outgroup, first, last, produce_log, verbose );
 695     }
 696
 697     private final static void calculateOrthologTable( final Phylogeny g, final boolean sort, final int counter )
 698             throws RIOException {
 699         if ( counter == 0 ) {
 700             final List<String> labels = new ArrayList<String>();
 701             final Set<String> labels_set = new HashSet<String>();
 702             for( final PhylogenyNode n : g.getExternalNodes() ) {
 703                 final String label = obtainLabel( labels_set, n );
 704                 labels_set.add( label );
 705                 labels.add( label );
 706             }
 707             if ( sort ) {
 708                 Collections.sort( labels );
 709             }
 710             _m = new IntMatrix( labels );
 711         }
 712         updateCounts( _m, counter, g );
 713     }
 714
 715     private final static void checkPreconditions( final IteratingPhylogenyParser p,
 716                                                   final Phylogeny species_tree,
 717                                                   final REROOTING rerooting,
 718                                                   final String outgroup,
 719                                                   final int first,
 720                                                   final int last ) throws RIOException, IOException {
 721         final Phylogeny g0 = p.next();
 722         if ( ( g0 == null ) || g0.isEmpty() ) {
 723             throw new RIOException( "input file does not seem to contain any gene trees" );
 724         }
 725         if ( g0.getNumberOfExternalNodes() < 2 ) {
 726             throw new RIOException( "input file does not seem to contain any useable gene trees" );
 727         }
 728         if ( !species_tree.isRooted() ) {
 729             throw new RIOException( "species tree is not rooted" );
 730         }
 731         if ( !( ( last == DEFAULT_RANGE ) && ( first == DEFAULT_RANGE ) )
 732                 && ( ( last < first ) || ( last < 0 ) || ( first < 0 ) ) ) {
 733             throw new RIOException( "attempt to set range (0-based) of gene to analyze to: from " + first + " to "
 734                     + last );
 735         }
 736         if ( ( rerooting == REROOTING.OUTGROUP ) && ForesterUtil.isEmpty( outgroup ) ) {
 737             throw new RIOException( "outgroup not set for midpoint rooting" );
 738         }
 739         if ( ( rerooting != REROOTING.OUTGROUP ) && !ForesterUtil.isEmpty( outgroup ) ) {
 740             throw new RIOException( "outgroup only used for midpoint rooting" );
 741         }
 742         if ( ( rerooting == REROOTING.MIDPOINT ) && ( PhylogenyMethods.calculateMaxDistanceToRoot( g0 ) <= 0 ) ) {
 743             throw new RIOException( "attempt to use midpoint rooting on gene trees which seem to have no (positive) branch lengths (cladograms)" );
 744         }
 745         if ( rerooting == REROOTING.OUTGROUP ) {
 746             try {
 747                 g0.getNode( outgroup );
 748             }
 749             catch ( final IllegalArgumentException e ) {
 750                 throw new RIOException( "cannot perform re-rooting by outgroup: " + e.getLocalizedMessage() );
 751             }
 752         }
 753     }
 754
 755     private final static void checkPreconditions( final Phylogeny[] gene_trees,
 756                                                   final Phylogeny species_tree,
 757                                                   final REROOTING rerooting,
 758                                                   final String outgroup,
 759                                                   final int first,
 760                                                   final int last ) throws RIOException {
 761         if ( !species_tree.isRooted() ) {
 762             throw new RIOException( "species tree is not rooted" );
 763         }
 764         if ( !( ( last == DEFAULT_RANGE ) && ( first == DEFAULT_RANGE ) )
 765                 && ( ( last < first ) || ( last >= gene_trees.length ) || ( last < 0 ) || ( first < 0 ) ) ) {
 766             throw new RIOException( "attempt to set range (0-based) of gene to analyze to: from " + first + " to "
 767                     + last + " (out of " + gene_trees.length + ")" );
 768         }
 769         if ( ( rerooting == REROOTING.OUTGROUP ) && ForesterUtil.isEmpty( outgroup ) ) {
 770             throw new RIOException( "outgroup not set for midpoint rooting" );
 771         }
 772         if ( ( rerooting != REROOTING.OUTGROUP ) && !ForesterUtil.isEmpty( outgroup ) ) {
 773             throw new RIOException( "outgroup only used for midpoint rooting" );
 774         }
 775         if ( ( rerooting == REROOTING.MIDPOINT )
 776                 && ( PhylogenyMethods.calculateMaxDistanceToRoot( gene_trees[ 0 ] ) <= 0 ) ) {
 777             throw new RIOException( "attempt to use midpoint rooting on gene trees which seem to have no (positive) branch lengths (cladograms)" );
 778         }
 779         if ( rerooting == REROOTING.OUTGROUP ) {
 780             try {
 781                 gene_trees[ 0 ].getNode( outgroup );
 782             }
 783             catch ( final IllegalArgumentException e ) {
 784                 throw new RIOException( "cannot perform re-rooting by outgroup: " + e.getLocalizedMessage() );
 785             }
 786         }
 787     }
 788
 789     private final static String obtainLabel( final Set<String> labels_set, final PhylogenyNode n ) throws RIOException {
 790         String label;
 791         if ( n.getNodeData().isHasSequence() && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getName() ) ) {
 792             label = n.getNodeData().getSequence().getName();
 793         }
 794         else if ( n.getNodeData().isHasSequence() && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getSymbol() ) ) {
 795             label = n.getNodeData().getSequence().getSymbol();
 796         }
 797         else if ( !ForesterUtil.isEmpty( n.getName() ) ) {
 798             label = n.getName();
 799         }
 800         else {
 801             throw new RIOException( "node " + n + " has no appropriate label" );
 802         }
 803         if ( labels_set.contains( label ) ) {
 804             throw new RIOException( "label " + label + " is not unique" );
 805         }
 806         return label;
 807     }
 808
 809     private final static Phylogeny[] parseGeneTrees( final File gene_trees_file ) throws FileNotFoundException,
 810             IOException {
 811         final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 812         final PhylogenyParser p = ParserUtils.createParserDependingOnFileType( gene_trees_file, true );
 813         if ( p instanceof NHXParser ) {
 814             final NHXParser nhx = ( NHXParser ) p;
 815             nhx.setReplaceUnderscores( false );
 816             nhx.setIgnoreQuotes( true );
 817             nhx.setTaxonomyExtraction( NHXParser.TAXONOMY_EXTRACTION.YES );
 818         }
 819         else if ( p instanceof NexusPhylogeniesParser ) {
 820             final NexusPhylogeniesParser nex = ( NexusPhylogeniesParser ) p;
 821             nex.setReplaceUnderscores( false );
 822             nex.setIgnoreQuotes( true );
 823             nex.setTaxonomyExtraction( TAXONOMY_EXTRACTION.YES );
 824         }
 825         return factory.create( gene_trees_file, p );
 826     }
 827
 828     private final static void removeSingleDescendentsNodes( final Phylogeny species_tree, final boolean verbose ) {
 829         final int o = PhylogenyMethods.countNumberOfOneDescendantNodes( species_tree );
 830         if ( o > 0 ) {
 831             if ( verbose ) {
 832                 System.out.println( "warning: species tree has " + o
 833                         + " internal nodes with only one descendent which are therefore going to be removed" );
 834             }
 835             PhylogenyMethods.deleteInternalNodesWithOnlyOneDescendent( species_tree );
 836         }
 837     }
 838
 839     private final static void updateCounts( final IntMatrix m, final int counter, final Phylogeny g )
 840             throws RIOException {
 841         PhylogenyMethods.preOrderReId( g );
 842         final HashMap<String, PhylogenyNode> map = PhylogenyMethods.createNameToExtNodeMap( g );
 843         for( int x = 0; x < m.size(); ++x ) {
 844             final String mx = m.getLabel( x );
 845             final PhylogenyNode nx = map.get( mx );
 846             if ( nx == null ) {
 847                 throw new RIOException( "node \"" + mx + "\" not present in gene tree #" + counter );
 848             }
 849             String my;
 850             PhylogenyNode ny;
 851             for( int y = 0; y < m.size(); ++y ) {
 852                 my = m.getLabel( y );
 853                 ny = map.get( my );
 854                 if ( ny == null ) {
 855                     throw new RIOException( "node \"" + my + "\" not present in gene tree #" + counter );
 856                 }
 857                 if ( !PhylogenyMethods.calculateLCAonTreeWithIdsInPreOrder( nx, ny ).isDuplication() ) {
 858                     m.inreaseByOne( x, y );
 859                 }
 860             }
 861         }
 862     }
 863
 864     public enum REROOTING {
 865         NONE, BY_ALGORITHM, MIDPOINT, OUTGROUP;
 866     }
 867 }