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