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.SDI.ALGORITHM;
  51 import org.forester.sdi.SDI.TaxonomyComparisonBase;
  52 import org.forester.sdi.SDIException;
  53 import org.forester.sdi.SDIR;
  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                    _samples;
  64     private int                    _ext_nodes;
  65     private TaxonomyComparisonBase _gsdir_tax_comp_base;
  66     private StringBuilder          _log;
  67     private boolean                _produce_log;
  68
  69     public RIO( final File gene_trees_file,
  70                 final Phylogeny species_tree,
  71                 final ALGORITHM algorithm,
  72                 final boolean produce_log ) throws IOException, SDIException, RIOException {
  73         init( produce_log );
  74         inferOrthologs( gene_trees_file, species_tree, algorithm );
  75     }
  76
  77     private final void init( final boolean produce_log ) {
  78         _produce_log = produce_log;
  79         _samples = -1;
  80         _ext_nodes = -1;
  81         _log = null;
  82         _gsdir_tax_comp_base = null;
  83         _analyzed_gene_trees = null;
  84         _removed_gene_tree_nodes = null;
  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 int getNumberOfSamples() {
 102         return _samples;
 103     }
 104
 105     public final List<PhylogenyNode> getRemovedGeneTreeNodes() {
 106         return _removed_gene_tree_nodes;
 107     }
 108
 109     private final void inferOrthologs( final File gene_trees_file,
 110                                        final Phylogeny species_tree,
 111                                        final ALGORITHM algorithm ) throws SDIException, RIOException,
 112             FileNotFoundException, IOException {
 113         // Read in first tree to get its sequence names
 114         // and strip species_tree.
 115         final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 116         final PhylogenyParser p = ParserUtils.createParserDependingOnFileType( gene_trees_file, true );
 117         if ( p instanceof NHXParser ) {
 118             final NHXParser nhx = ( NHXParser ) p;
 119             nhx.setReplaceUnderscores( false );
 120             nhx.setIgnoreQuotes( true );
 121             nhx.setTaxonomyExtraction( NHXParser.TAXONOMY_EXTRACTION.YES );
 122         }
 123         final Phylogeny[] gene_trees = factory.create( gene_trees_file, p );
 124         if ( algorithm == ALGORITHM.SDIR ) {
 125             // Removes from species_tree all species not found in gene_tree.
 126             PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( gene_trees[ 0 ], species_tree );
 127             if ( species_tree.isEmpty() ) {
 128                 throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 129             }
 130         }
 131         if ( _produce_log ) {
 132             _log = new StringBuilder();
 133         }
 134         _analyzed_gene_trees = new Phylogeny[ gene_trees.length ];
 135         int i = 0;
 136         int gene_tree_ext_nodes = 0;
 137         if ( _produce_log ) {
 138             _log.append( "#" );
 139             _log.append( "\t" );
 140             _log.append( "with minimal number of duplications" );
 141             _log.append( "/" );
 142             _log.append( "root placements" );
 143             _log.append( "\t[" );
 144             _log.append( "min" );
 145             _log.append( "-" );
 146             _log.append( "max" );
 147             _log.append( "]" );
 148             _log.append( ForesterUtil.LINE_SEPARATOR );
 149         }
 150         for( final Phylogeny gt : gene_trees ) {
 151             if ( algorithm == ALGORITHM.SDIR ) {
 152                 // Removes from gene_tree all species not found in species_tree.
 153                 PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gt );
 154                 if ( gt.isEmpty() ) {
 155                     throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 156                 }
 157                 if ( i == 0 ) {
 158                     gene_tree_ext_nodes = gt.getNumberOfExternalNodes();
 159                 }
 160                 else if ( gene_tree_ext_nodes != gt.getNumberOfExternalNodes() ) {
 161                     throw new RIOException( "(cleaned up) gene tree #" + ( i + 1 )
 162                             + " has a different number of external nodes (" + gt.getNumberOfExternalNodes()
 163                             + ") than those gene trees preceding it (" + gene_tree_ext_nodes + ")" );
 164                 }
 165             }
 166             _analyzed_gene_trees[ i ] = performOrthologInference( gt, species_tree, algorithm, i );
 167             ++i;
 168         }
 169         setNumberOfSamples( gene_trees.length );
 170     }
 171
 172     private final Phylogeny performOrthologInference( final Phylogeny gene_tree,
 173                                                       final Phylogeny species_tree,
 174                                                       final ALGORITHM algorithm,
 175                                                       final int i ) throws SDIException, RIOException {
 176         final Phylogeny assigned_tree;
 177         switch ( algorithm ) {
 178             case SDIR: {
 179                 final SDIR sdir = new SDIR();
 180                 assigned_tree = sdir.infer( gene_tree,
 181                                             species_tree,
 182                                             false,
 183                                             RIO.ROOT_BY_MINIMIZING_SUM_OF_DUPS,
 184                                             RIO.ROOT_BY_MINIMIZING_TREE_HEIGHT,
 185                                             true,
 186                                             1 )[ 0 ];
 187                 break;
 188             }
 189             case GSDIR: {
 190                 //  System.out.println( "gene/species tree size before: " + gene_tree.getNumberOfExternalNodes() + "/"
 191                 //         + species_tree.getNumberOfExternalNodes() );
 192                 final GSDIR gsdir = new GSDIR( gene_tree, species_tree, true, i == 0 );
 193                 // System.out.println( "gene/species tree size before: " + gene_tree.getNumberOfExternalNodes() + "/"
 194                 //         + species_tree.getNumberOfExternalNodes() );
 195                 assigned_tree = gsdir.getMinDuplicationsSumGeneTrees().get( 0 );
 196                 if ( i == 0 ) {
 197                     _removed_gene_tree_nodes = gsdir.getStrippedExternalGeneTreeNodes();
 198                 }
 199                 if ( _produce_log ) {
 200                     final BasicDescriptiveStatistics stats = gsdir.getDuplicationsSumStats();
 201                     _log.append( i );
 202                     _log.append( "\t" );
 203                     _log.append( gsdir.getMinDuplicationsSumGeneTrees().size() );
 204                     _log.append( "/" );
 205                     _log.append( stats.getN() );
 206                     _log.append( "\t[" );
 207                     _log.append( ( int ) stats.getMin() );
 208                     _log.append( "-" );
 209                     _log.append( ( int ) stats.getMax() );
 210                     _log.append( "]" );
 211                     _log.append( ForesterUtil.LINE_SEPARATOR );
 212                 }
 213                 _gsdir_tax_comp_base = gsdir.getTaxCompBase();
 214                 break;
 215             }
 216             default: {
 217                 throw new IllegalArgumentException( "illegal algorithm: " + algorithm );
 218             }
 219         }
 220         setExtNodesOfAnalyzedGeneTrees( assigned_tree.getNumberOfExternalNodes() );
 221         return assigned_tree;
 222     }
 223
 224     private final void setExtNodesOfAnalyzedGeneTrees( final int i ) {
 225         _ext_nodes = i;
 226     }
 227
 228     private final void setNumberOfSamples( final int i ) {
 229         _samples = i;
 230     }
 231
 232     public final static IntMatrix calculateOrthologTable( final Phylogeny[] analyzed_gene_trees, final boolean sort )
 233             throws RIOException {
 234         final List<String> labels = new ArrayList<String>();
 235         final Set<String> labels_set = new HashSet<String>();
 236         String label;
 237         for( final PhylogenyNode n : analyzed_gene_trees[ 0 ].getExternalNodes() ) {
 238             if ( n.getNodeData().isHasSequence() && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getName() ) ) {
 239                 label = n.getNodeData().getSequence().getName();
 240             }
 241             else if ( n.getNodeData().isHasSequence()
 242                     && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getSymbol() ) ) {
 243                 label = n.getNodeData().getSequence().getSymbol();
 244             }
 245             else if ( !ForesterUtil.isEmpty( n.getName() ) ) {
 246                 label = n.getName();
 247             }
 248             else {
 249                 throw new IllegalArgumentException( "node " + n + " has no appropriate label" );
 250             }
 251             if ( labels_set.contains( label ) ) {
 252                 throw new IllegalArgumentException( "label " + label + " is not unique" );
 253             }
 254             labels_set.add( label );
 255             labels.add( label );
 256         }
 257         if ( sort ) {
 258             Collections.sort( labels );
 259         }
 260         final IntMatrix m = new IntMatrix( labels );
 261         int counter = 0;
 262         for( final Phylogeny gt : analyzed_gene_trees ) {
 263             counter++;
 264             PhylogenyMethods.preOrderReId( gt );
 265             final HashMap<String, PhylogenyNode> map = PhylogenyMethods.createNameToExtNodeMap( gt );
 266             for( int x = 0; x < m.size(); ++x ) {
 267                 final String mx = m.getLabel( x );
 268                 final PhylogenyNode nx = map.get( mx );
 269                 if ( nx == null ) {
 270                     throw new RIOException( "node \"" + mx + "\" not present in gene tree #" + counter );
 271                 }
 272                 String my;
 273                 PhylogenyNode ny;
 274                 for( int y = 0; y < m.size(); ++y ) {
 275                     my = m.getLabel( y );
 276                     ny = map.get( my );
 277                     if ( ny == null ) {
 278                         throw new RIOException( "node \"" + my + "\" not present in gene tree #" + counter );
 279                     }
 280                     if ( !PhylogenyMethods.calculateLCAonTreeWithIdsInPreOrder( nx, ny ).isDuplication() ) {
 281                         m.inreaseByOne( x, y );
 282                     }
 283                 }
 284             }
 285         }
 286         return m;
 287     }
 288
 289     public final TaxonomyComparisonBase getGSDIRtaxCompBase() {
 290         return _gsdir_tax_comp_base;
 291     }
 292
 293     public final StringBuilder getLog() {
 294         return _log;
 295     }
 296 }