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
  40 import org.forester.datastructures.IntMatrix;
  41 import org.forester.io.parsers.PhylogenyParser;
  42 import org.forester.io.parsers.nhx.NHXParser;
  43 import org.forester.io.parsers.util.ParserUtils;
  44 import org.forester.phylogeny.Phylogeny;
  45 import org.forester.phylogeny.PhylogenyMethods;
  46 import org.forester.phylogeny.PhylogenyNode;
  47 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
  48 import org.forester.phylogeny.factories.PhylogenyFactory;
  49 import org.forester.sdi.GSDIR;
  50 import org.forester.sdi.SDIException;
  51 import org.forester.sdi.SDIR;
  52 import org.forester.sdi.SDIutil.ALGORITHM;
  53 import org.forester.sdi.SDIutil.TaxonomyComparisonBase;
  54 import org.forester.util.BasicDescriptiveStatistics;
  55 import org.forester.util.ForesterUtil;
  56
  57 public final class RIO {
  58
  59     private final static boolean   ROOT_BY_MINIMIZING_SUM_OF_DUPS = true;
  60     private final static boolean   ROOT_BY_MINIMIZING_TREE_HEIGHT = true;
  61     private Phylogeny[]            _analyzed_gene_trees;
  62     private List<PhylogenyNode>    _removed_gene_tree_nodes;
  63     private int                    _ext_nodes;
  64     private TaxonomyComparisonBase _gsdir_tax_comp_base;
  65     private StringBuilder          _log;
  66     private boolean                _produce_log;
  67     private boolean                _verbose;
  68
  69     public RIO( final File gene_trees_file,
  70                 final Phylogeny species_tree,
  71                 final ALGORITHM algorithm,
  72                 final boolean produce_log,
  73                 final boolean verbose ) throws IOException, SDIException, RIOException {
  74         init( produce_log, verbose );
  75         inferOrthologs( gene_trees_file, species_tree, algorithm );
  76     }
  77
  78     public RIO( final Phylogeny[] gene_trees,
  79                 final Phylogeny species_tree,
  80                 final ALGORITHM algorithm,
  81                 final boolean produce_log,
  82                 final boolean verbose ) throws IOException, SDIException, RIOException {
  83         init( produce_log, verbose );
  84         inferOrthologs( gene_trees, species_tree, algorithm );
  85     }
  86
  87     public final Phylogeny[] getAnalyzedGeneTrees() {
  88         return _analyzed_gene_trees;
  89     }
  90
  91     /**
  92      * Returns the numbers of number of ext nodes in gene trees analyzed (after
  93      * stripping).
  94      *
  95      * @return number of ext nodes in gene trees analyzed (after stripping)
  96      */
  97     public final int getExtNodesOfAnalyzedGeneTrees() {
  98         return _ext_nodes;
  99     }
 100
 101     public final TaxonomyComparisonBase getGSDIRtaxCompBase() {
 102         return _gsdir_tax_comp_base;
 103     }
 104
 105     public final StringBuilder getLog() {
 106         return _log;
 107     }
 108
 109     public final List<PhylogenyNode> getRemovedGeneTreeNodes() {
 110         return _removed_gene_tree_nodes;
 111     }
 112
 113     private final void inferOrthologs( final File gene_trees_file,
 114                                        final Phylogeny species_tree,
 115                                        final ALGORITHM algorithm ) throws SDIException, RIOException,
 116             FileNotFoundException, IOException {
 117         final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 118         final PhylogenyParser p = ParserUtils.createParserDependingOnFileType( gene_trees_file, true );
 119         if ( p instanceof NHXParser ) {
 120             final NHXParser nhx = ( NHXParser ) p;
 121             nhx.setReplaceUnderscores( false );
 122             nhx.setIgnoreQuotes( true );
 123             nhx.setTaxonomyExtraction( NHXParser.TAXONOMY_EXTRACTION.YES );
 124         }
 125         final Phylogeny[] gene_trees = factory.create( gene_trees_file, p );
 126         inferOrthologs( gene_trees, species_tree, algorithm );
 127     }
 128
 129     private final void inferOrthologs( final Phylogeny[] gene_trees,
 130                                        final Phylogeny species_tree,
 131                                        final ALGORITHM algorithm ) throws SDIException, RIOException,
 132             FileNotFoundException, IOException {
 133         if ( algorithm == ALGORITHM.SDIR ) {
 134             // Removes from species_tree all species not found in gene_tree.
 135             PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( gene_trees[ 0 ], species_tree );
 136             if ( species_tree.isEmpty() ) {
 137                 throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 138             }
 139         }
 140         if ( _produce_log ) {
 141             _log = new StringBuilder();
 142             writeLogSubHeader();
 143         }
 144         _analyzed_gene_trees = new Phylogeny[ gene_trees.length ];
 145         int gene_tree_ext_nodes = 0;
 146         if ( _verbose ) {
 147             System.out.println();
 148         }
 149         for( int i = 0; i < gene_trees.length; ++i ) {
 150             final Phylogeny gt = gene_trees[ i ];
 151             if ( _verbose ) {
 152                 ForesterUtil.updateProgress( ( double ) i / gene_trees.length );
 153             }
 154             if ( i == 0 ) {
 155                 gene_tree_ext_nodes = gt.getNumberOfExternalNodes();
 156             }
 157             else if ( gene_tree_ext_nodes != gt.getNumberOfExternalNodes() ) {
 158                 throw new RIOException( "gene tree #" + ( i + 1 ) + " has a different number of external nodes ("
 159                         + gt.getNumberOfExternalNodes() + ") than the preceding gene trees (" + gene_tree_ext_nodes
 160                         + ")" );
 161             }
 162             if ( algorithm == ALGORITHM.SDIR ) {
 163                 // Removes from gene_tree all species not found in species_tree.
 164                 PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gt );
 165                 if ( gt.isEmpty() ) {
 166                     throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 167                 }
 168             }
 169             _analyzed_gene_trees[ i ] = performOrthologInference( gt, species_tree, algorithm, i );
 170         }
 171         if ( _verbose ) {
 172             System.out.println();
 173             System.out.println();
 174         }
 175     }
 176
 177     private final void init( final boolean produce_log, final boolean verbose ) {
 178         _produce_log = produce_log;
 179         _verbose = verbose;
 180         _ext_nodes = -1;
 181         _log = null;
 182         _gsdir_tax_comp_base = null;
 183         _analyzed_gene_trees = null;
 184         _removed_gene_tree_nodes = null;
 185     }
 186
 187     private final Phylogeny performOrthologInference( final Phylogeny gene_tree,
 188                                                       final Phylogeny species_tree,
 189                                                       final ALGORITHM algorithm,
 190                                                       final int i ) throws SDIException, RIOException {
 191         final Phylogeny assigned_tree;
 192         switch ( algorithm ) {
 193             case SDIR: {
 194                 final SDIR sdir = new SDIR();
 195                 assigned_tree = sdir.infer( gene_tree,
 196                                             species_tree,
 197                                             false,
 198                                             RIO.ROOT_BY_MINIMIZING_SUM_OF_DUPS,
 199                                             RIO.ROOT_BY_MINIMIZING_TREE_HEIGHT,
 200                                             true,
 201                                             1 )[ 0 ];
 202                 break;
 203             }
 204             case GSDIR: {
 205                 //  System.out.println( "gene/species tree size before: " + gene_tree.getNumberOfExternalNodes() + "/"
 206                 //         + species_tree.getNumberOfExternalNodes() );
 207                 final GSDIR gsdir = new GSDIR( gene_tree, species_tree, true, i == 0 );
 208                 // System.out.println( "gene/species tree size before: " + gene_tree.getNumberOfExternalNodes() + "/"
 209                 //         + species_tree.getNumberOfExternalNodes() );
 210                 assigned_tree = gsdir.getMinDuplicationsSumGeneTrees().get( 0 );
 211                 if ( i == 0 ) {
 212                     _removed_gene_tree_nodes = gsdir.getStrippedExternalGeneTreeNodes();
 213                     for( final PhylogenyNode r : _removed_gene_tree_nodes ) {
 214                         if ( !r.getNodeData().isHasTaxonomy() ) {
 215                             throw new RIOException( "node with no (appropriate) taxonomic information found in gene tree #1: "
 216                                     + r.toString() );
 217                         }
 218                     }
 219                 }
 220                 if ( _produce_log ) {
 221                     writeStatsToLog( i, gsdir );
 222                 }
 223                 _gsdir_tax_comp_base = gsdir.getTaxCompBase();
 224                 break;
 225             }
 226             default: {
 227                 throw new IllegalArgumentException( "illegal algorithm: " + algorithm );
 228             }
 229         }
 230         if ( i == 0 ) {
 231             _ext_nodes = assigned_tree.getNumberOfExternalNodes();
 232         }
 233         else if ( _ext_nodes != assigned_tree.getNumberOfExternalNodes() ) {
 234             throw new RIOException( "after stripping gene tree #" + ( i + 1 )
 235                     + " has a different number of external nodes (" + assigned_tree.getNumberOfExternalNodes()
 236                     + ") than the preceding gene trees (" + _ext_nodes + ")" );
 237         }
 238         return assigned_tree;
 239     }
 240
 241     private void writeLogSubHeader() {
 242         _log.append( "#" );
 243         _log.append( "\t" );
 244         _log.append( "with minimal number of duplications" );
 245         _log.append( "/" );
 246         _log.append( "root placements" );
 247         _log.append( "\t[" );
 248         _log.append( "min" );
 249         _log.append( "-" );
 250         _log.append( "max" );
 251         _log.append( "]" );
 252         _log.append( ForesterUtil.LINE_SEPARATOR );
 253     }
 254
 255     private final void writeStatsToLog( final int i, final GSDIR gsdir ) {
 256         final BasicDescriptiveStatistics stats = gsdir.getDuplicationsSumStats();
 257         _log.append( i );
 258         _log.append( "\t" );
 259         _log.append( gsdir.getMinDuplicationsSumGeneTrees().size() );
 260         _log.append( "/" );
 261         _log.append( stats.getN() );
 262         _log.append( "\t[" );
 263         _log.append( ( int ) stats.getMin() );
 264         _log.append( "-" );
 265         _log.append( ( int ) stats.getMax() );
 266         _log.append( "]" );
 267         _log.append( ForesterUtil.LINE_SEPARATOR );
 268     }
 269
 270     public final static IntMatrix calculateOrthologTable( final Phylogeny[] analyzed_gene_trees, final boolean sort )
 271             throws RIOException {
 272         final List<String> labels = new ArrayList<String>();
 273         final Set<String> labels_set = new HashSet<String>();
 274         String label;
 275         for( final PhylogenyNode n : analyzed_gene_trees[ 0 ].getExternalNodes() ) {
 276             if ( n.getNodeData().isHasSequence() && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getName() ) ) {
 277                 label = n.getNodeData().getSequence().getName();
 278             }
 279             else if ( n.getNodeData().isHasSequence()
 280                     && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getSymbol() ) ) {
 281                 label = n.getNodeData().getSequence().getSymbol();
 282             }
 283             else if ( !ForesterUtil.isEmpty( n.getName() ) ) {
 284                 label = n.getName();
 285             }
 286             else {
 287                 throw new RIOException( "node " + n + " has no appropriate label" );
 288             }
 289             if ( labels_set.contains( label ) ) {
 290                 throw new RIOException( "label " + label + " is not unique" );
 291             }
 292             labels_set.add( label );
 293             labels.add( label );
 294         }
 295         if ( sort ) {
 296             Collections.sort( labels );
 297         }
 298         final IntMatrix m = new IntMatrix( labels );
 299         int counter = 0;
 300         for( final Phylogeny gt : analyzed_gene_trees ) {
 301             counter++;
 302             PhylogenyMethods.preOrderReId( gt );
 303             final HashMap<String, PhylogenyNode> map = PhylogenyMethods.createNameToExtNodeMap( gt );
 304             for( int x = 0; x < m.size(); ++x ) {
 305                 final String mx = m.getLabel( x );
 306                 final PhylogenyNode nx = map.get( mx );
 307                 if ( nx == null ) {
 308                     throw new RIOException( "node \"" + mx + "\" not present in gene tree #" + counter );
 309                 }
 310                 String my;
 311                 PhylogenyNode ny;
 312                 for( int y = 0; y < m.size(); ++y ) {
 313                     my = m.getLabel( y );
 314                     ny = map.get( my );
 315                     if ( ny == null ) {
 316                         throw new RIOException( "node \"" + my + "\" not present in gene tree #" + counter );
 317                     }
 318                     if ( !PhylogenyMethods.calculateLCAonTreeWithIdsInPreOrder( nx, ny ).isDuplication() ) {
 319                         m.inreaseByOne( x, y );
 320                     }
 321                 }
 322             }
 323         }
 324         return m;
 325     }
 326 }