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: www.phylosoft.org/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.PhylogenyParser;
  44 import org.forester.io.parsers.nhx.NHXParser;
  45 import org.forester.io.parsers.util.ParserUtils;
  46 import org.forester.phylogeny.Phylogeny;
  47 import org.forester.phylogeny.PhylogenyMethods;
  48 import org.forester.phylogeny.PhylogenyNode;
  49 import org.forester.phylogeny.data.Taxonomy;
  50 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
  51 import org.forester.phylogeny.factories.PhylogenyFactory;
  52 import org.forester.sdi.GSDI;
  53 import org.forester.sdi.GSDIR;
  54 import org.forester.sdi.SDIException;
  55 import org.forester.sdi.SDIR;
  56 import org.forester.sdi.SDIutil.ALGORITHM;
  57 import org.forester.sdi.SDIutil.TaxonomyComparisonBase;
  58 import org.forester.util.BasicDescriptiveStatistics;
  59 import org.forester.util.ForesterUtil;
  60
  61 public final class RIO {
  62
  63     private Phylogeny[]                      _analyzed_gene_trees;
  64     private List<PhylogenyNode>              _removed_gene_tree_nodes;
  65     private int                              _ext_nodes;
  66     private TaxonomyComparisonBase           _gsdir_tax_comp_base;
  67     private final StringBuilder              _log;
  68     private final BasicDescriptiveStatistics _duplications_stats;
  69     private final boolean                    _produce_log;
  70     private final boolean                    _verbose;
  71     private final REROOTING                  _rerooting;
  72
  73     private RIO( final Phylogeny[] gene_trees,
  74                  final Phylogeny species_tree,
  75                  final ALGORITHM algorithm,
  76                  final REROOTING rerooting,
  77                  final String outgroup,
  78                  final int first,
  79                  final int last,
  80                  final boolean produce_log,
  81                  final boolean verbose ) throws IOException, SDIException, RIOException {
  82         if ( !ForesterUtil.isEmpty( outgroup ) && ( rerooting != REROOTING.OUTGROUP ) ) {
  83             throw new IllegalArgumentException( "can only use outgroup when re-rooting by outgroup" );
  84         }
  85         if ( !( ( last == -1 ) && ( first == -1 ) )
  86                 && ( ( last < first ) || ( last >= gene_trees.length ) || ( first >= gene_trees.length ) || ( last < 0 ) || ( first < 0 ) ) ) {
  87             throw new IllegalArgumentException( "gene tree range is out of range: " + first + "-" + last );
  88         }
  89         _produce_log = produce_log;
  90         _verbose = verbose;
  91         _rerooting = rerooting;
  92         _ext_nodes = -1;
  93         _log = new StringBuilder();
  94         _gsdir_tax_comp_base = null;
  95         _analyzed_gene_trees = null;
  96         _removed_gene_tree_nodes = null;
  97         _duplications_stats = new BasicDescriptiveStatistics();
  98         inferOrthologs( gene_trees, species_tree, algorithm, outgroup, first, last );
  99     }
 100
 101     public final Phylogeny[] getAnalyzedGeneTrees() {
 102         return _analyzed_gene_trees;
 103     }
 104
 105     public final BasicDescriptiveStatistics getDuplicationsStatistics() {
 106         return _duplications_stats;
 107     }
 108
 109     /**
 110      * Returns the numbers of number of ext nodes in gene trees analyzed (after
 111      * stripping).
 112      *
 113      * @return number of ext nodes in gene trees analyzed (after stripping)
 114      */
 115     public final int getExtNodesOfAnalyzedGeneTrees() {
 116         return _ext_nodes;
 117     }
 118
 119     public final TaxonomyComparisonBase getGSDIRtaxCompBase() {
 120         return _gsdir_tax_comp_base;
 121     }
 122
 123     public final StringBuilder getLog() {
 124         return _log;
 125     }
 126
 127     public final List<PhylogenyNode> getRemovedGeneTreeNodes() {
 128         return _removed_gene_tree_nodes;
 129     }
 130
 131     private final void inferOrthologs( final Phylogeny[] gene_trees,
 132                                        final Phylogeny species_tree,
 133                                        final ALGORITHM algorithm,
 134                                        final String outgroup,
 135                                        final int first,
 136                                        final int last ) throws SDIException, RIOException, FileNotFoundException,
 137             IOException {
 138         if ( algorithm == ALGORITHM.SDIR ) {
 139             // Removes from species_tree all species not found in gene_tree.
 140             PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( gene_trees[ 0 ], species_tree );
 141             if ( species_tree.isEmpty() ) {
 142                 throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 143             }
 144         }
 145         if ( log() ) {
 146             preLog( gene_trees, species_tree, algorithm, outgroup );
 147         }
 148         final Phylogeny[] my_gene_trees;
 149         if ( ( first >= 0 ) && ( last >= first ) && ( last < gene_trees.length ) ) {
 150             if ( log() ) {
 151                 log( "Gene tree range: " + first + "-" + last );
 152             }
 153             my_gene_trees = new Phylogeny[ 1 + last - first ];
 154             int c = 0;
 155             for( int i = first; i <= last; ++i ) {
 156                 my_gene_trees[ c++ ] = gene_trees[ i ];
 157             }
 158         }
 159         else {
 160             my_gene_trees = gene_trees;
 161         }
 162         if ( _verbose && ( my_gene_trees.length > 10 ) ) {
 163             System.out.println();
 164         }
 165         _analyzed_gene_trees = new Phylogeny[ my_gene_trees.length ];
 166         int gene_tree_ext_nodes = 0;
 167         for( int i = 0; i < my_gene_trees.length; ++i ) {
 168             final Phylogeny gt = my_gene_trees[ i ];
 169             if ( _verbose && ( my_gene_trees.length > 10 ) ) {
 170                 ForesterUtil.updateProgress( ( ( double ) i ) / my_gene_trees.length );
 171             }
 172             if ( i == 0 ) {
 173                 gene_tree_ext_nodes = gt.getNumberOfExternalNodes();
 174             }
 175             else if ( gene_tree_ext_nodes != gt.getNumberOfExternalNodes() ) {
 176                 throw new RIOException( "gene tree #" + ( i + 1 ) + " has a different number of external nodes ("
 177                         + gt.getNumberOfExternalNodes() + ") than the preceding gene trees (" + gene_tree_ext_nodes
 178                         + ")" );
 179             }
 180             if ( algorithm == ALGORITHM.SDIR ) {
 181                 // Removes from gene_tree all species not found in species_tree.
 182                 PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gt );
 183                 if ( gt.isEmpty() ) {
 184                     throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 185                 }
 186             }
 187             _analyzed_gene_trees[ i ] = performOrthologInference( gt, species_tree, algorithm, outgroup, i );
 188         }
 189         if ( log() ) {
 190             postLog( species_tree );
 191         }
 192         if ( _verbose && ( my_gene_trees.length > 10 ) ) {
 193             System.out.println();
 194             System.out.println();
 195         }
 196     }
 197
 198     private final boolean log() {
 199         return _produce_log;
 200     }
 201
 202     private final void log( final String s ) {
 203         _log.append( s );
 204         _log.append( ForesterUtil.LINE_SEPARATOR );
 205     }
 206
 207     private final void logRemovedGeneTreeNodes() {
 208         log( "Species stripped from gene trees:" );
 209         final SortedSet<String> rn = new TreeSet<String>();
 210         for( final PhylogenyNode n : getRemovedGeneTreeNodes() ) {
 211             final Taxonomy t = n.getNodeData().getTaxonomy();
 212             switch ( getGSDIRtaxCompBase() ) {
 213                 case CODE: {
 214                     rn.add( t.getTaxonomyCode() );
 215                     break;
 216                 }
 217                 case ID: {
 218                     rn.add( t.getIdentifier().toString() );
 219                     break;
 220                 }
 221                 case SCIENTIFIC_NAME: {
 222                     rn.add( t.getScientificName() );
 223                     break;
 224                 }
 225             }
 226         }
 227         for( final String s : rn ) {
 228             log( s );
 229         }
 230         log( "" );
 231     }
 232
 233     private final Phylogeny performOrthologInference( final Phylogeny gene_tree,
 234                                                       final Phylogeny species_tree,
 235                                                       final ALGORITHM algorithm,
 236                                                       final String outgroup,
 237                                                       final int i ) throws SDIException, RIOException {
 238         final Phylogeny assigned_tree;
 239         switch ( algorithm ) {
 240             case SDIR: {
 241                 assigned_tree = performOrthologInferenceBySDI( gene_tree, species_tree );
 242                 break;
 243             }
 244             case GSDIR: {
 245                 assigned_tree = performOrthologInferenceByGSDI( gene_tree, species_tree, outgroup, i );
 246                 break;
 247             }
 248             default: {
 249                 throw new IllegalArgumentException( "illegal algorithm: " + algorithm );
 250             }
 251         }
 252         if ( i == 0 ) {
 253             _ext_nodes = assigned_tree.getNumberOfExternalNodes();
 254         }
 255         else if ( _ext_nodes != assigned_tree.getNumberOfExternalNodes() ) {
 256             throw new RIOException( "after stripping gene tree #" + ( i + 1 )
 257                     + " has a different number of external nodes (" + assigned_tree.getNumberOfExternalNodes()
 258                     + ") than the preceding gene trees (" + _ext_nodes + ")" );
 259         }
 260         return assigned_tree;
 261     }
 262
 263     private final Phylogeny performOrthologInferenceByGSDI( final Phylogeny gene_tree,
 264                                                             final Phylogeny species_tree,
 265                                                             final String outgroup,
 266                                                             final int i ) throws SDIException, RIOException {
 267         final Phylogeny assigned_tree;
 268         if ( _rerooting == REROOTING.BY_ALGORITHM ) {
 269             final GSDIR gsdir = new GSDIR( gene_tree, species_tree, true, i == 0 );
 270             final List<Phylogeny> assigned_trees = gsdir.getMinDuplicationsSumGeneTrees();
 271             if ( i == 0 ) {
 272                 _removed_gene_tree_nodes = gsdir.getStrippedExternalGeneTreeNodes();
 273                 for( final PhylogenyNode r : _removed_gene_tree_nodes ) {
 274                     if ( !r.getNodeData().isHasTaxonomy() ) {
 275                         throw new RIOException( "node with no (appropriate) taxonomic information found in gene tree #1: "
 276                                 + r.toString() );
 277                     }
 278                 }
 279             }
 280             final List<Integer> shortests = GSDIR.getIndexesOfShortestTree( assigned_trees );
 281             assigned_tree = assigned_trees.get( shortests.get( 0 ) );
 282             if ( log() ) {
 283                 writeStatsToLog( i, gsdir, shortests );
 284             }
 285             if ( i == 0 ) {
 286                 _gsdir_tax_comp_base = gsdir.getTaxCompBase();
 287             }
 288             _duplications_stats.addValue( gsdir.getMinDuplicationsSum() );
 289         }
 290         else {
 291             if ( _rerooting == REROOTING.MIDPOINT ) {
 292                 PhylogenyMethods.midpointRoot( gene_tree );
 293             }
 294             else if ( _rerooting == REROOTING.OUTGROUP ) {
 295                 PhylogenyNode n;
 296                 try {
 297                     n = gene_tree.getNode( outgroup );
 298                 }
 299                 catch ( final IllegalArgumentException e ) {
 300                     throw new RIOException( "failed to perform re-rooting by outgroup: " + e.getLocalizedMessage() );
 301                 }
 302                 gene_tree.reRoot( n );
 303             }
 304             final GSDI gsdi = new GSDI( gene_tree, species_tree, true, true, true );
 305             _removed_gene_tree_nodes = gsdi.getStrippedExternalGeneTreeNodes();
 306             for( final PhylogenyNode r : _removed_gene_tree_nodes ) {
 307                 if ( !r.getNodeData().isHasTaxonomy() ) {
 308                     throw new RIOException( "node with no (appropriate) taxonomic information found in gene tree #1: "
 309                             + r.toString() );
 310                 }
 311             }
 312             assigned_tree = gene_tree;
 313             if ( i == 0 ) {
 314                 _gsdir_tax_comp_base = gsdi.getTaxCompBase();
 315             }
 316             _duplications_stats.addValue( gsdi.getDuplicationsSum() );
 317         }
 318         return assigned_tree;
 319     }
 320
 321     private final Phylogeny performOrthologInferenceBySDI( final Phylogeny gene_tree, final Phylogeny species_tree )
 322             throws SDIException {
 323         final SDIR sdir = new SDIR();
 324         return sdir.infer( gene_tree, species_tree, false, true, true, true, 1 )[ 0 ];
 325     }
 326
 327     private final void postLog( final Phylogeny species_tree ) {
 328         log( "" );
 329         if ( getRemovedGeneTreeNodes().size() > 0 ) {
 330             logRemovedGeneTreeNodes();
 331         }
 332         log( "Species tree external nodes (after stripping)   : " + species_tree.getNumberOfExternalNodes() );
 333         log( "Species tree polytomies (after stripping)       : "
 334                 + PhylogenyMethods.countNumberOfPolytomies( species_tree ) );
 335         log( "Taxonomy linking based on                       : " + getGSDIRtaxCompBase() );
 336         final java.text.DecimalFormat df = new java.text.DecimalFormat( "0.#" );
 337         log( "Gene trees analyzed                             : " + _duplications_stats.getN() );
 338         log( "Mean number of duplications                     : " + df.format( _duplications_stats.arithmeticMean() )
 339                 + " (sd: " + df.format( _duplications_stats.sampleStandardDeviation() ) + ")" );
 340         if ( _duplications_stats.getN() > 3 ) {
 341             log( "Median number of duplications                   : " + df.format( _duplications_stats.median() ) );
 342         }
 343         log( "Minimum duplications                            : " + ( int ) _duplications_stats.getMin() );
 344         log( "Maximum duplications                            : " + ( int ) _duplications_stats.getMax() );
 345     }
 346
 347     private final void preLog( final Phylogeny[] gene_trees,
 348                                final Phylogeny species_tree,
 349                                final ALGORITHM algorithm,
 350                                final String outgroup ) {
 351         log( "Number of gene tree (total)                     : " + gene_trees.length );
 352         log( "Algorithm                                       : " + algorithm );
 353         log( "Species tree external nodes (prior to stripping): " + species_tree.getNumberOfExternalNodes() );
 354         log( "Species tree polytomies (prior to stripping)    : "
 355                 + PhylogenyMethods.countNumberOfPolytomies( species_tree ) );
 356         String rs = "";
 357         switch ( _rerooting ) {
 358             case BY_ALGORITHM: {
 359                 rs = "minimizing duplications";
 360                 break;
 361             }
 362             case MIDPOINT: {
 363                 rs = "midpoint";
 364                 break;
 365             }
 366             case OUTGROUP: {
 367                 rs = "outgroup: " + outgroup;
 368                 break;
 369             }
 370             case NONE: {
 371                 rs = "none";
 372                 break;
 373             }
 374         }
 375         log( "Re-rooting                                      : " + rs );
 376         if ( _rerooting == REROOTING.BY_ALGORITHM ) {
 377             writeLogSubHeader();
 378         }
 379     }
 380
 381     private final void writeLogSubHeader() {
 382         _log.append( ForesterUtil.LINE_SEPARATOR );
 383         _log.append( "Some information about duplication numbers in gene trees:" );
 384         _log.append( ForesterUtil.LINE_SEPARATOR );
 385         _log.append( "#" );
 386         _log.append( "\t" );
 387         _log.append( "re-rootings with minimal number of duplications" );
 388         _log.append( "/" );
 389         _log.append( "total root placements" );
 390         _log.append( "\t" );
 391         _log.append( "duplications range" );
 392         _log.append( "\t" );
 393         _log.append( "mininal duplication re-rootings with shortest tree heigth" );
 394         _log.append( ForesterUtil.LINE_SEPARATOR );
 395     }
 396
 397     private final void writeStatsToLog( final int i, final GSDIR gsdir, final List<Integer> shortests ) {
 398         final BasicDescriptiveStatistics stats = gsdir.getDuplicationsSumStats();
 399         _log.append( i );
 400         _log.append( "\t" );
 401         _log.append( gsdir.getMinDuplicationsSumGeneTrees().size() );
 402         _log.append( "/" );
 403         _log.append( stats.getN() );
 404         _log.append( "\t" );
 405         _log.append( ( int ) stats.getMin() );
 406         _log.append( "-" );
 407         _log.append( ( int ) stats.getMax() );
 408         _log.append( "\t" );
 409         _log.append( shortests.size() );
 410         _log.append( ForesterUtil.LINE_SEPARATOR );
 411     }
 412
 413     public final static IntMatrix calculateOrthologTable( final Phylogeny[] analyzed_gene_trees, final boolean sort )
 414             throws RIOException {
 415         final List<String> labels = new ArrayList<String>();
 416         final Set<String> labels_set = new HashSet<String>();
 417         String label;
 418         for( final PhylogenyNode n : analyzed_gene_trees[ 0 ].getExternalNodes() ) {
 419             if ( n.getNodeData().isHasSequence() && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getName() ) ) {
 420                 label = n.getNodeData().getSequence().getName();
 421             }
 422             else if ( n.getNodeData().isHasSequence()
 423                     && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getSymbol() ) ) {
 424                 label = n.getNodeData().getSequence().getSymbol();
 425             }
 426             else if ( !ForesterUtil.isEmpty( n.getName() ) ) {
 427                 label = n.getName();
 428             }
 429             else {
 430                 throw new RIOException( "node " + n + " has no appropriate label" );
 431             }
 432             if ( labels_set.contains( label ) ) {
 433                 throw new RIOException( "label " + label + " is not unique" );
 434             }
 435             labels_set.add( label );
 436             labels.add( label );
 437         }
 438         if ( sort ) {
 439             Collections.sort( labels );
 440         }
 441         final IntMatrix m = new IntMatrix( labels );
 442         int counter = 0;
 443         for( final Phylogeny gt : analyzed_gene_trees ) {
 444             counter++;
 445             PhylogenyMethods.preOrderReId( gt );
 446             final HashMap<String, PhylogenyNode> map = PhylogenyMethods.createNameToExtNodeMap( gt );
 447             for( int x = 0; x < m.size(); ++x ) {
 448                 final String mx = m.getLabel( x );
 449                 final PhylogenyNode nx = map.get( mx );
 450                 if ( nx == null ) {
 451                     throw new RIOException( "node \"" + mx + "\" not present in gene tree #" + counter );
 452                 }
 453                 String my;
 454                 PhylogenyNode ny;
 455                 for( int y = 0; y < m.size(); ++y ) {
 456                     my = m.getLabel( y );
 457                     ny = map.get( my );
 458                     if ( ny == null ) {
 459                         throw new RIOException( "node \"" + my + "\" not present in gene tree #" + counter );
 460                     }
 461                     if ( !PhylogenyMethods.calculateLCAonTreeWithIdsInPreOrder( nx, ny ).isDuplication() ) {
 462                         m.inreaseByOne( x, y );
 463                     }
 464                 }
 465             }
 466         }
 467         return m;
 468     }
 469
 470     public final static RIO executeAnalysis( final File gene_trees_file,
 471                                              final Phylogeny species_tree,
 472                                              final ALGORITHM algorithm,
 473                                              final REROOTING rerooting,
 474                                              final String outgroup,
 475                                              final boolean produce_log,
 476                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 477         final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 478         final PhylogenyParser p = ParserUtils.createParserDependingOnFileType( gene_trees_file, true );
 479         if ( p instanceof NHXParser ) {
 480             final NHXParser nhx = ( NHXParser ) p;
 481             nhx.setReplaceUnderscores( false );
 482             nhx.setIgnoreQuotes( true );
 483             nhx.setTaxonomyExtraction( NHXParser.TAXONOMY_EXTRACTION.YES );
 484         }
 485         final Phylogeny[] gene_trees = factory.create( gene_trees_file, p );
 486         return new RIO( gene_trees, species_tree, algorithm, rerooting, outgroup, -1, -1, produce_log, verbose );
 487     }
 488
 489     public final static RIO executeAnalysis( final Phylogeny[] gene_trees, final Phylogeny species_tree )
 490             throws IOException, SDIException, RIOException {
 491         return new RIO( gene_trees, species_tree, ALGORITHM.GSDIR, REROOTING.BY_ALGORITHM, null, -1, -1, false, false );
 492     }
 493
 494     public final static RIO executeAnalysis( final Phylogeny[] gene_trees,
 495                                              final Phylogeny species_tree,
 496                                              final ALGORITHM algorithm,
 497                                              final REROOTING rerooting,
 498                                              final String outgroup,
 499                                              final boolean produce_log,
 500                                              final boolean verbose ) throws IOException, SDIException, RIOException {
 501         return new RIO( gene_trees, species_tree, algorithm, rerooting, outgroup, -1, -1, produce_log, verbose );
 502     }
 503
 504     public enum REROOTING {
 505         NONE, BY_ALGORITHM, MIDPOINT, OUTGROUP;
 506     }
 507 }