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