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.ForesterUtil;
  55
  56 public final class RIO {
  57
  58     private final static boolean   ROOT_BY_MINIMIZING_SUM_OF_DUPS = true;
  59     private final static boolean   ROOT_BY_MINIMIZING_TREE_HEIGHT = true;
  60     private Phylogeny[]            _analyzed_gene_trees;
  61     private List<PhylogenyNode>    _removed_gene_tree_nodes;
  62     private int                    _samples;
  63     private int                    _ext_nodes;
  64     private TaxonomyComparisonBase _gsdir_tax_comp_base;
  65
  66     public RIO( final File gene_trees_file, final Phylogeny species_tree, final ALGORITHM algorithm )
  67             throws IOException, SDIException, RIOException {
  68         init();
  69         inferOrthologs( gene_trees_file, species_tree, algorithm );
  70     }
  71
  72     public final Phylogeny[] getAnalyzedGeneTrees() {
  73         return _analyzed_gene_trees;
  74     }
  75
  76     /**
  77      * Returns the numbers of number of ext nodes in gene trees analyzed (after
  78      * stripping).
  79      *
  80      * @return number of ext nodes in gene trees analyzed (after stripping)
  81      */
  82     public final int getExtNodesOfAnalyzedGeneTrees() {
  83         return _ext_nodes;
  84     }
  85
  86     public final int getNumberOfSamples() {
  87         return _samples;
  88     }
  89
  90     public final List<PhylogenyNode> getRemovedGeneTreeNodes() {
  91         return _removed_gene_tree_nodes;
  92     }
  93
  94     private final void inferOrthologs( final File gene_trees_file,
  95                                        final Phylogeny species_tree,
  96                                        final ALGORITHM algorithm ) throws SDIException, RIOException,
  97             FileNotFoundException, IOException {
  98         // Read in first tree to get its sequence names
  99         // and strip species_tree.
 100         final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 101         final PhylogenyParser p = ParserUtils.createParserDependingOnFileType( gene_trees_file, true );
 102         if ( p instanceof NHXParser ) {
 103             final NHXParser nhx = ( NHXParser ) p;
 104             nhx.setReplaceUnderscores( false );
 105             nhx.setIgnoreQuotes( true );
 106             nhx.setTaxonomyExtraction( NHXParser.TAXONOMY_EXTRACTION.YES );
 107         }
 108         final Phylogeny[] gene_trees = factory.create( gene_trees_file, p );
 109         // Removes from species_tree all species not found in gene_tree.
 110         final List<PhylogenyNode> _removed_species_tree_ext_nodes = PhylogenyMethods
 111                 .taxonomyBasedDeletionOfExternalNodes( gene_trees[ 0 ], species_tree );
 112         if ( species_tree.isEmpty() ) {
 113             throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 114         }
 115         _analyzed_gene_trees = new Phylogeny[ gene_trees.length ];
 116         int i = 0;
 117         int gene_tree_ext_nodes = 0;
 118         for( final Phylogeny gt : gene_trees ) {
 119             if ( algorithm == ALGORITHM.SDIR ) {
 120                 // Removes from gene_tree all species not found in species_tree.
 121                 PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gt );
 122                 if ( gt.isEmpty() ) {
 123                     throw new RIOException( "failed to establish species based mapping between gene and species trees" );
 124                 }
 125                 if ( i == 0 ) {
 126                     gene_tree_ext_nodes = gt.getNumberOfExternalNodes();
 127                 }
 128                 else if ( gene_tree_ext_nodes != gt.getNumberOfExternalNodes() ) {
 129                     throw new RIOException( "(cleaned up) gene tree #" + ( i + 1 )
 130                             + " has a different number of external nodes (" + gt.getNumberOfExternalNodes()
 131                             + ") than those gene trees preceding it (" + gene_tree_ext_nodes + ")" );
 132                 }
 133             }
 134             _analyzed_gene_trees[ i ] = performOrthologInference( gt, species_tree, algorithm, i );
 135             ++i;
 136         }
 137         setNumberOfSamples( gene_trees.length );
 138     }
 139
 140     private final void init() {
 141         _samples = 1;
 142         _ext_nodes = 0;
 143     }
 144
 145     private final Phylogeny performOrthologInference( final Phylogeny gene_tree,
 146                                                       final Phylogeny species_tree,
 147                                                       final ALGORITHM algorithm,
 148                                                       final int i ) throws SDIException, RIOException {
 149         final Phylogeny assigned_tree;
 150         switch ( algorithm ) {
 151             case SDIR: {
 152                 final SDIR sdir = new SDIR();
 153                 assigned_tree = sdir.infer( gene_tree,
 154                                             species_tree,
 155                                             false,
 156                                             RIO.ROOT_BY_MINIMIZING_SUM_OF_DUPS,
 157                                             RIO.ROOT_BY_MINIMIZING_TREE_HEIGHT,
 158                                             true,
 159                                             1 )[ 0 ];
 160                 break;
 161             }
 162             case GSDIR: {
 163                 final GSDIR gsdir = new GSDIR( gene_tree, species_tree, true, i == 0 );
 164                 assigned_tree = gsdir.getMinDuplicationsSumGeneTrees().get( 0 );
 165                 _gsdir_tax_comp_base = gsdir.getTaxCompBase();
 166                 break;
 167             }
 168             default: {
 169                 throw new IllegalArgumentException( "illegal algorithm: " + algorithm );
 170             }
 171         }
 172         setExtNodesOfAnalyzedGeneTrees( assigned_tree.getNumberOfExternalNodes() );
 173         return assigned_tree;
 174     }
 175
 176     private final void setExtNodesOfAnalyzedGeneTrees( final int i ) {
 177         _ext_nodes = i;
 178     }
 179
 180     private final void setNumberOfSamples( int i ) {
 181         if ( i < 1 ) {
 182             i = 1;
 183         }
 184         _samples = i;
 185     }
 186
 187     public final static IntMatrix calculateOrthologTable( final Phylogeny[] analyzed_gene_trees, final boolean sort )
 188             throws RIOException {
 189         final List<String> labels = new ArrayList<String>();
 190         final Set<String> labels_set = new HashSet<String>();
 191         String label;
 192         for( final PhylogenyNode n : analyzed_gene_trees[ 0 ].getExternalNodes() ) {
 193             if ( n.getNodeData().isHasSequence() && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getName() ) ) {
 194                 label = n.getNodeData().getSequence().getName();
 195             }
 196             else if ( n.getNodeData().isHasSequence()
 197                     && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getSymbol() ) ) {
 198                 label = n.getNodeData().getSequence().getSymbol();
 199             }
 200             else if ( !ForesterUtil.isEmpty( n.getName() ) ) {
 201                 label = n.getName();
 202             }
 203             else {
 204                 throw new IllegalArgumentException( "node " + n + " has no appropriate label" );
 205             }
 206             if ( labels_set.contains( label ) ) {
 207                 throw new IllegalArgumentException( "label " + label + " is not unique" );
 208             }
 209             labels_set.add( label );
 210             labels.add( label );
 211         }
 212         if ( sort ) {
 213             Collections.sort( labels );
 214         }
 215         final IntMatrix m = new IntMatrix( labels );
 216         int counter = 0;
 217         for( final Phylogeny gt : analyzed_gene_trees ) {
 218             counter++;
 219             PhylogenyMethods.preOrderReId( gt );
 220             final HashMap<String, PhylogenyNode> map = PhylogenyMethods.createNameToExtNodeMap( gt );
 221             for( int x = 0; x < m.size(); ++x ) {
 222                 final String mx = m.getLabel( x );
 223                 final PhylogenyNode nx = map.get( mx );
 224                 if ( nx == null ) {
 225                     throw new RIOException( "node \"" + mx + "\" not present in gene tree #" + counter );
 226                 }
 227                 String my;
 228                 PhylogenyNode ny;
 229                 for( int y = 0; y < m.size(); ++y ) {
 230                     my = m.getLabel( y );
 231                     ny = map.get( my );
 232                     if ( ny == null ) {
 233                         throw new RIOException( "node \"" + my + "\" not present in gene tree #" + counter );
 234                     }
 235                     if ( !PhylogenyMethods.calculateLCAonTreeWithIdsInPreOrder( nx, ny ).isDuplication() ) {
 236                         m.inreaseByOne( x, y );
 237                     }
 238                 }
 239             }
 240         }
 241         return m;
 242     }
 243 }