forester/java/src/org/forester/sdi/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.sdi;
  29
  30 import java.io.File;
  31 import java.io.IOException;
  32 import java.util.ArrayList;
  33 import java.util.Arrays;
  34 import java.util.HashMap;
  35 import java.util.HashSet;
  36 import java.util.List;
  37 import java.util.Set;
  38
  39 import org.forester.datastructures.IntMatrix;
  40 import org.forester.evoinference.matrix.distance.DistanceMatrix;
  41 import org.forester.io.parsers.PhylogenyParser;
  42 import org.forester.io.parsers.SymmetricalDistanceMatrixParser;
  43 import org.forester.io.parsers.nhx.NHXParser;
  44 import org.forester.io.parsers.util.ParserUtils;
  45 import org.forester.phylogeny.Phylogeny;
  46 import org.forester.phylogeny.PhylogenyMethods;
  47 import org.forester.phylogeny.PhylogenyNode;
  48 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
  49 import org.forester.phylogeny.factories.PhylogenyFactory;
  50 import org.forester.phylogeny.iterators.PhylogenyNodeIterator;
  51 import org.forester.util.ForesterUtil;
  52
  53 /*
  54  * @author Christian M. Zmasek
  55  */
  56 public final class RIO {
  57
  58     private final static boolean                      ROOT_BY_MINIMIZING_MAPPING_COST = false;
  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 final static boolean                      TIME                            = false;
  62     private HashMap<String, HashMap<String, Integer>> _o_hash_maps;
  63     private HashMap<String, HashMap<String, Integer>> _so_hash_maps;
  64     private HashMap<String, HashMap<String, Integer>> _up_hash_maps;
  65     private HashMap<String, HashMap<String, Integer>> _sn_hash_maps;                          // HashMap of HashMaps
  66     private DistanceMatrix                            _m;
  67     private HashMap<String, Double>                   _l;
  68     private List<String>                              _seq_names;
  69     private int                                       _bootstraps;
  70     private int                                       _ext_nodes_;
  71     private long                                      _time;
  72
  73     /**
  74      * Default constructor.
  75      */
  76     public RIO() {
  77         reset();
  78     }
  79
  80     public static IntMatrix calculateOrthologTable( final Phylogeny[] gene_trees ) {
  81         final List<String> labels = new ArrayList<String>();
  82         final Set<String> labels_set = new HashSet<String>();
  83         String label;
  84         for( final PhylogenyNode n : gene_trees[ 0 ].getExternalNodes() ) {
  85             if ( n.getNodeData().isHasSequence() && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getName() ) ) {
  86                 label = n.getNodeData().getSequence().getName();
  87             }
  88             else if ( n.getNodeData().isHasSequence()
  89                     && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getSymbol() ) ) {
  90                 label = n.getNodeData().getSequence().getSymbol();
  91             }
  92             else if ( !ForesterUtil.isEmpty( n.getName() ) ) {
  93                 label = n.getName();
  94             }
  95             else {
  96                 throw new IllegalArgumentException( "node " + n + " has no appropriate label" );
  97             }
  98             if ( labels_set.contains( label ) ) {
  99                 throw new IllegalArgumentException( "label " + label + " is not unique" );
 100             }
 101             labels_set.add( label );
 102             labels.add( label );
 103         }
 104         final IntMatrix m = new IntMatrix( labels );
 105         int counter = 0;
 106         for( final Phylogeny gt : gene_trees ) {
 107             System.out.println( counter );
 108             counter++;
 109             PhylogenyMethods.preOrderReId( gt );
 110             final HashMap<String, PhylogenyNode> map = PhylogenyMethods.createNameToExtNodeMap( gt );
 111             for( int x = 0; x < m.size(); ++x ) {
 112                 final PhylogenyNode nx = map.get( m.getLabel( x ) );
 113                 for( int y = 0; y < m.size(); ++y ) {
 114                     if ( !PhylogenyMethods.calculateLCAonTreeWithIdsInPreOrder( nx, map.get( m.getLabel( y ) ) ).isDuplication() ) {
 115                         m.inreaseByOne( x, y );
 116                     }
 117                 }
 118             }
 119         }
 120         return m;
 121     }
 122
 123     /**
 124      * Returns the numbers of trees analyzed.
 125      *
 126      * @return the numbers of trees analyzed
 127      */
 128     public final int getBootstraps() {
 129         return _bootstraps;
 130     }
 131
 132     // Helper method for inferredOrthologsToString.
 133     // inferredOrthologsToArrayList,
 134     // and inferredUltraParalogsToString.
 135     private final double getBootstrapValueFromHash( final HashMap<String, Integer> h, final String name ) {
 136         if ( !h.containsKey( name ) ) {
 137             return 0.0;
 138         }
 139         final int i = h.get( name );
 140         return ( ( i * 100.0 ) / getBootstraps() );
 141     }
 142
 143     /**
 144      * Returns the distance to a sequences/taxa after a distance list file has
 145      * been read in with readDistanceList(File). Throws an exception if name is
 146      * not found or if no list has been read in.
 147      *
 148      * @param name
 149      *            a sequence name
 150      */
 151     public final double getDistance( String name ) {
 152         double distance = 0.0;
 153         name = name.trim();
 154         if ( _l == null ) {
 155             throw new RuntimeException( "Distance list has probably not been read in (successfully)." );
 156         }
 157         if ( _l.get( name ) == null ) {
 158             throw new IllegalArgumentException( name + " not found." );
 159         }
 160         distance = ( _l.get( name ) ).doubleValue();
 161         return distance;
 162     }
 163
 164     public final double getDistance( final String name1, final String name2 ) {
 165         try {
 166             return _m.getValue( _m.getIndex( name1 ), _m.getIndex( name2 ) );
 167         }
 168         catch ( final Exception e ) {
 169             return 1;
 170         }
 171     }
 172
 173     /**
 174      * Returns the numbers of number of ext nodes in gene trees analyzed (after
 175      * stripping).
 176      *
 177      * @return number of ext nodes in gene trees analyzed (after stripping)
 178      */
 179     public final int getExtNodesOfAnalyzedGeneTrees() {
 180         return _ext_nodes_;
 181     }
 182
 183     /**
 184      * Returns a HashMap containing the inferred orthologs of the external gene
 185      * tree node with the sequence name seq_name. Sequence names are the keys
 186      * (String), numbers of observations are the values (Int). Orthologs are to
 187      * be inferred by method "inferOrthologs". Throws an exception if seq_name
 188      * is not found.
 189      *
 190      * @param seq_name
 191      *            sequence name of a external node of the gene trees
 192      * @return HashMap containing the inferred orthologs
 193      *         (name(String)->value(Int))
 194      */
 195     public final HashMap<String, Integer> getInferredOrthologs( final String seq_name ) {
 196         if ( _o_hash_maps == null ) {
 197             return null;
 198         }
 199         return _o_hash_maps.get( seq_name );
 200     }
 201
 202     private final HashMap<String, Integer> getInferredSubtreeNeighbors( final String seq_name ) {
 203         if ( _sn_hash_maps == null ) {
 204             return null;
 205         }
 206         return _sn_hash_maps.get( seq_name );
 207     }
 208
 209     /**
 210      * Returns a HashMap containing the inferred "super orthologs" of the
 211      * external gene tree node with the sequence name seq_name. Sequence names
 212      * are the keys (String), numbers of observations are the values (Int).
 213      * Super orthologs are to be inferred by method "inferOrthologs". Throws an
 214      * exception if seq_name is not found.
 215      *
 216      * @param seq_name
 217      *            sequence name of a external node of the gene trees
 218      * @return HashMap containing the inferred super orthologs
 219      *         (name(String)->value(Int))
 220      */
 221     public final HashMap<String, Integer> getInferredSuperOrthologs( final String seq_name ) {
 222         if ( _so_hash_maps == null ) {
 223             return null;
 224         }
 225         return _so_hash_maps.get( seq_name );
 226     }
 227
 228     /**
 229      * Returns a HashMap containing the inferred "ultra paralogs" of the
 230      * external gene tree node with the sequence name seq_name. Sequence names
 231      * are the keys (String), numbers of observations are the values (Int).
 232      * "ultra paralogs" are to be inferred by method "inferOrthologs". Throws an
 233      * exception if seq_name is not found.
 234      *
 235      * @param seq_name
 236      *            sequence name of a external node of the gene trees
 237      * @return HashMap containing the inferred ultra paralogs
 238      *         (name(String)->value(Int))
 239      */
 240     public final HashMap<String, Integer> getInferredUltraParalogs( final String seq_name ) {
 241         if ( _up_hash_maps == null ) {
 242             return null;
 243         }
 244         return _up_hash_maps.get( seq_name );
 245     }
 246
 247     /**
 248      * Returns the time (in ms) needed to run "inferOrthologs". Final variable
 249      * TIME needs to be set to true.
 250      *
 251      * @return time (in ms) needed to run method "inferOrthologs"
 252      */
 253     public long getTime() {
 254         return _time;
 255     }
 256
 257     /**
 258      * Infers the orthologs (as well the "super orthologs", the "subtree
 259      * neighbors", and the "ultra paralogs") for each external node of the gene
 260      * Trees in multiple tree File gene_trees_file (=output of PHYLIP NEIGHBOR,
 261      * for example). Tallies how many times each sequence is (super-)
 262      * orthologous towards the query. Tallies how many times each sequence is
 263      * ultra paralogous towards the query. Tallies how many times each sequence
 264      * is a subtree neighbor of the query. Gene duplications are inferred using
 265      * SDI. Modifies its argument species_tree. Is a little faster than
 266      * "inferOrthologs(File,Phylogeny)" since orthologs are only inferred for
 267      * query.
 268      * <p>
 269      * To obtain the results use the methods listed below.
 270      *
 271      * @param gene_trees_file
 272      *            a File containing gene Trees in NH format, which is the result
 273      *            of performing a bootstrap analysis in PHYLIP
 274      * @param species_tree
 275      *            a species Phylogeny, which has species names in its species
 276      *            fields
 277      * @param query
 278      *            the sequence name of the squence whose orthologs are to be
 279      *            inferred
 280      * @throws SDIException
 281      */
 282     public void inferOrthologs( final File gene_trees_file, final Phylogeny species_tree, final String query )
 283             throws IOException, SDIException {
 284         int bs = 0;
 285         if ( RIO.TIME ) {
 286             _time = System.currentTimeMillis();
 287         }
 288         // Read in first tree to get its sequence names
 289         // and strip species_tree.
 290         final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 291         final PhylogenyParser p = ParserUtils.createParserDependingOnFileType( gene_trees_file, true );
 292         if ( p instanceof NHXParser ) {
 293             final NHXParser nhx = ( NHXParser ) p;
 294             nhx.setReplaceUnderscores( false );
 295             nhx.setIgnoreQuotes( true );
 296             nhx.setTaxonomyExtraction( PhylogenyMethods.TAXONOMY_EXTRACTION.YES );
 297         }
 298         final Phylogeny gene_tree = factory.create( gene_trees_file, p )[ 0 ];
 299         System.out.println( "species " + species_tree.toString() );
 300         // Removes from species_tree all species not found in gene_tree.
 301         PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( gene_tree, species_tree );
 302         PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gene_tree );
 303         _seq_names = getAllExternalSequenceNames( gene_tree );
 304         if ( ( _seq_names == null ) || ( _seq_names.size() < 1 ) ) {
 305             throw new IOException( "could not get sequence names" );
 306         }
 307         _o_hash_maps = new HashMap<String, HashMap<String, Integer>>();
 308         _so_hash_maps = new HashMap<String, HashMap<String, Integer>>();
 309         _up_hash_maps = new HashMap<String, HashMap<String, Integer>>();
 310         _sn_hash_maps = new HashMap<String, HashMap<String, Integer>>();
 311         _o_hash_maps.put( query, new HashMap<String, Integer>( _seq_names.size() ) );
 312         _so_hash_maps.put( query, new HashMap<String, Integer>( _seq_names.size() ) );
 313         _up_hash_maps.put( query, new HashMap<String, Integer>( _seq_names.size() ) );
 314         _sn_hash_maps.put( query, new HashMap<String, Integer>( _seq_names.size() ) );
 315         // Go through all gene trees in the file.
 316         final Phylogeny[] gene_trees = factory.create( gene_trees_file, p );
 317         final Phylogeny[] assigned_trees = new Phylogeny[ gene_trees.length ];
 318         int c = 0;
 319         for( final Phylogeny gt : gene_trees ) {
 320             bs++;
 321             // Removes from gene_tree all species not found in species_tree.
 322             PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gt );
 323             assigned_trees[ c++ ] = inferOrthologsHelper( gt, species_tree, query );
 324             // System.out.println( bs );
 325         }
 326         final IntMatrix m = calculateOrthologTable( assigned_trees );
 327         System.out.println( m.toString() );
 328         setBootstraps( bs );
 329         if ( RIO.TIME ) {
 330             _time = ( System.currentTimeMillis() - _time );
 331         }
 332     }
 333
 334     public List<PhylogenyNode> getNodesViaSequenceName( final Phylogeny phy, final String seq_name ) {
 335         final List<PhylogenyNode> nodes = new ArrayList<PhylogenyNode>();
 336         for( final PhylogenyNodeIterator iter = phy.iteratorPreorder(); iter.hasNext(); ) {
 337             final PhylogenyNode n = iter.next();
 338             if ( n.getNodeData().isHasSequence() && n.getNodeData().getSequence().getName().equals( seq_name ) ) {
 339                 nodes.add( n );
 340             }
 341             if ( !n.getNodeData().isHasSequence() && n.getName().equals( seq_name ) ) {
 342                 nodes.add( n );
 343             }
 344         }
 345         return nodes;
 346     }
 347
 348     // Helper method which performs the actual ortholog inference for
 349     // the external node with seqname query.
 350     private Phylogeny inferOrthologsHelper( final Phylogeny gene_tree, final Phylogeny species_tree, final String query )
 351             throws SDIException {
 352         Phylogeny assigned_tree = null;
 353         List<PhylogenyNode> nodes = null;
 354         final SDIR sdiunrooted = new SDIR();
 355         List<PhylogenyNode> orthologs = null;
 356         List<PhylogenyNode> super_orthologs = null;
 357         List<PhylogenyNode> ultra_paralogs = null;
 358         List<PhylogenyNode> subtree_neighbors = null;
 359         assigned_tree = sdiunrooted.infer( gene_tree,
 360                                            species_tree,
 361                                            RIO.ROOT_BY_MINIMIZING_MAPPING_COST,
 362                                            RIO.ROOT_BY_MINIMIZING_SUM_OF_DUPS,
 363                                            RIO.ROOT_BY_MINIMIZING_TREE_HEIGHT,
 364                                            true,
 365                                            1 )[ 0 ];
 366         setExtNodesOfAnalyzedGeneTrees( assigned_tree.getNumberOfExternalNodes() );
 367         nodes = getNodesViaSequenceName( assigned_tree, query );
 368         if ( nodes.size() > 1 ) {
 369             throw new IllegalArgumentException( "node named [" + query + "] not unique" );
 370         }
 371         else if ( nodes.isEmpty() ) {
 372             throw new IllegalArgumentException( "no node containing a sequence named [" + query + "] found" );
 373         }
 374         final PhylogenyNode query_node = nodes.get( 0 );
 375         orthologs = PhylogenyMethods.getOrthologousNodes( assigned_tree, query_node );
 376         updateHash( _o_hash_maps, query, orthologs );
 377         super_orthologs = PhylogenyMethods.getSuperOrthologousNodes( query_node );
 378         updateHash( _so_hash_maps, query, super_orthologs );
 379         subtree_neighbors = getSubtreeNeighbors( query_node, 2 );
 380         updateHash( _sn_hash_maps, query, subtree_neighbors );
 381         ultra_paralogs = PhylogenyMethods.getUltraParalogousNodes( query_node );
 382         updateHash( _up_hash_maps, query, ultra_paralogs );
 383         return assigned_tree;
 384     }
 385
 386     /**
 387      * Returns an ArrayList containg the names of orthologs of the PhylogenyNode
 388      * with seq name seq_name.
 389      *
 390      * @param seq_name
 391      *            sequence name of a external node of the gene trees
 392      * @param threshold_orthologs
 393      *            the minimal number of observations for a a sequence to be
 394      *            reported as orthologous as percentage (0.0-100.0%)
 395      * @return ArrayList containg the names of orthologs of the PhylogenyNode
 396      *         with seq name seq_name
 397      */
 398     public ArrayList<String> inferredOrthologsToArrayList( final String seq_name, double threshold_orthologs ) {
 399         HashMap<String, Integer> o_hashmap = null;
 400         String name = null;
 401         double o = 0.0;
 402         final ArrayList<String> arraylist = new ArrayList<String>();
 403         if ( _o_hash_maps == null ) {
 404             throw new RuntimeException( "Orthologs have not been calculated (successfully)." );
 405         }
 406         if ( threshold_orthologs < 0.0 ) {
 407             threshold_orthologs = 0.0;
 408         }
 409         else if ( threshold_orthologs > 100.0 ) {
 410             threshold_orthologs = 100.0;
 411         }
 412         o_hashmap = getInferredOrthologs( seq_name );
 413         if ( o_hashmap == null ) {
 414             throw new RuntimeException( "Orthologs for " + seq_name + " were not established." );
 415         }
 416         if ( _seq_names.size() > 0 ) {
 417             I: for( int i = 0; i < _seq_names.size(); ++i ) {
 418                 name = _seq_names.get( i );
 419                 if ( name.equals( seq_name ) ) {
 420                     continue I;
 421                 }
 422                 o = getBootstrapValueFromHash( o_hashmap, name );
 423                 if ( o < threshold_orthologs ) {
 424                     continue I;
 425                 }
 426                 arraylist.add( name );
 427             }
 428         }
 429         return arraylist;
 430     }
 431
 432     /**
 433      * Returns a String containg the names of orthologs of the PhylogenyNode
 434      * with seq name query_name. The String also contains how many times a
 435      * particular ortholog has been observed.
 436      * <p>
 437      * <ul>
 438      * The output order is (per line): Name, Ortholog, Subtree neighbor, Super
 439      * ortholog, Distance
 440      * </ul>
 441      * <p>
 442      * The sort priority of this is determined by sort in the following manner:
 443      * <ul>
 444      * <li>0 : Ortholog
 445      * <li>1 : Ortholog, Super ortholog
 446      * <li>2 : Super ortholog, Ortholog
 447      * <li>3 : Ortholog, Distance
 448      * <li>4 : Distance, Ortholog
 449      * <li>5 : Ortholog, Super ortholog, Distance
 450      * <li>6 : Ortholog, Distance, Super ortholog
 451      * <li>7 : Super ortholog, Ortholog, Distance
 452      * <li>8 : Super ortholog, Distance, Ortholog
 453      * <li>9 : Distance, Ortholog, Super ortholog
 454      * <li>10 : Distance, Super ortholog, Ortholog
 455      * <li>11 : Ortholog, Subtree neighbor, Distance
 456      * <li>12 : Ortholog, Subtree neighbor, Super ortholog, Distance (default)
 457      * <li>13 : Ortholog, Super ortholog, Subtree neighbor, Distance
 458      * <li>14 : Subtree neighbor, Ortholog, Super ortholog, Distance
 459      * <li>15 : Subtree neighbor, Distance, Ortholog, Super ortholog
 460      * <li>16 : Ortholog, Distance, Subtree neighbor, Super ortholog
 461      * <li>17 : Ortholog, Subtree neighbor, Distance, Super ortholog
 462      * </ul>
 463      * <p>
 464      * Returns "-" if no putative orthologs have been found (given
 465      * threshold_orthologs).
 466      * <p>
 467      * Orthologs are to be inferred by method "inferOrthologs".
 468      * <p>
 469      * (Last modified: 05/08/01)
 470      *
 471      * @param query_name
 472      *            sequence name of a external node of the gene trees
 473      * @param sort
 474      *            order and sort priority
 475      * @param threshold_orthologs
 476      *            the minimal number of observations for a a sequence to be
 477      *            reported as orthologous, in percents (0.0-100.0%)
 478      * @param threshold_subtreeneighborings
 479      *            the minimal number of observations for a a sequence to be
 480      *            reported as orthologous, in percents (0.0-100.0%)
 481      * @return String containing the inferred orthologs, String containing "-"
 482      *         if no orthologs have been found null in case of error
 483      * @see #inferOrthologs(File,Phylogeny,String)
 484      * @see #inferOrthologs(Phylogeny[],Phylogeny)
 485      * @see #inferOrthologs(File,Phylogeny)
 486      * @see #getOrder(int)
 487      */
 488     public StringBuffer inferredOrthologsToString( final String query_name,
 489                                                    int sort,
 490                                                    double threshold_orthologs,
 491                                                    double threshold_subtreeneighborings ) {
 492         HashMap<String, Integer> o_hashmap = null;
 493         HashMap<String, Integer> s_hashmap = null;
 494         HashMap<String, Integer> n_hashmap = null;
 495         String name = "";
 496         double o = 0.0, // Orthologs.
 497         s = 0.0, // Super orthologs.
 498         sn = 0.0, // Subtree neighbors.
 499         value1 = 0.0, value2 = 0.0, value3 = 0.0, value4 = 0.0, d = 0.0;
 500         final ArrayList<Tuplet> nv = new ArrayList<Tuplet>();
 501         if ( ( _o_hash_maps == null ) || ( _so_hash_maps == null ) || ( _sn_hash_maps == null ) ) {
 502             throw new RuntimeException( "Orthologs have not been calculated (successfully)" );
 503         }
 504         if ( ( sort < 0 ) || ( sort > 17 ) ) {
 505             sort = 12;
 506         }
 507         if ( ( sort > 2 ) && ( _m == null ) && ( _l == null ) ) {
 508             throw new RuntimeException( "Distance list or matrix have not been read in (successfully)" );
 509         }
 510         if ( threshold_orthologs < 0.0 ) {
 511             threshold_orthologs = 0.0;
 512         }
 513         else if ( threshold_orthologs > 100.0 ) {
 514             threshold_orthologs = 100.0;
 515         }
 516         if ( threshold_subtreeneighborings < 0.0 ) {
 517             threshold_subtreeneighborings = 0.0;
 518         }
 519         else if ( threshold_subtreeneighborings > 100.0 ) {
 520             threshold_subtreeneighborings = 100.0;
 521         }
 522         o_hashmap = getInferredOrthologs( query_name );
 523         s_hashmap = getInferredSuperOrthologs( query_name );
 524         n_hashmap = getInferredSubtreeNeighbors( query_name );
 525         if ( ( o_hashmap == null ) || ( s_hashmap == null ) || ( n_hashmap == null ) ) {
 526             throw new RuntimeException( "Orthologs for " + query_name + " were not established" );
 527         }
 528         final StringBuffer orthologs = new StringBuffer();
 529         if ( _seq_names.size() > 0 ) {
 530             I: for( int i = 0; i < _seq_names.size(); ++i ) {
 531                 name = _seq_names.get( i );
 532                 if ( name.equals( query_name ) ) {
 533                     continue I;
 534                 }
 535                 o = getBootstrapValueFromHash( o_hashmap, name );
 536                 if ( o < threshold_orthologs ) {
 537                     continue I;
 538                 }
 539                 sn = getBootstrapValueFromHash( n_hashmap, name );
 540                 if ( sn < threshold_subtreeneighborings ) {
 541                     continue I;
 542                 }
 543                 s = getBootstrapValueFromHash( s_hashmap, name );
 544                 if ( sort >= 3 ) {
 545                     if ( _m != null ) {
 546                         d = getDistance( query_name, name );
 547                     }
 548                     else {
 549                         d = getDistance( name );
 550                     }
 551                 }
 552                 switch ( sort ) {
 553                     case 0:
 554                         nv.add( new Tuplet( name, o, 5 ) );
 555                         break;
 556                     case 1:
 557                         nv.add( new Tuplet( name, o, s, 5 ) );
 558                         break;
 559                     case 2:
 560                         nv.add( new Tuplet( name, s, o, 5 ) );
 561                         break;
 562                     case 3:
 563                         nv.add( new Tuplet( name, o, d, 1 ) );
 564                         break;
 565                     case 4:
 566                         nv.add( new Tuplet( name, d, o, 0 ) );
 567                         break;
 568                     case 5:
 569                         nv.add( new Tuplet( name, o, s, d, 2 ) );
 570                         break;
 571                     case 6:
 572                         nv.add( new Tuplet( name, o, d, s, 1 ) );
 573                         break;
 574                     case 7:
 575                         nv.add( new Tuplet( name, s, o, d, 2 ) );
 576                         break;
 577                     case 8:
 578                         nv.add( new Tuplet( name, s, d, o, 1 ) );
 579                         break;
 580                     case 9:
 581                         nv.add( new Tuplet( name, d, o, s, 0 ) );
 582                         break;
 583                     case 10:
 584                         nv.add( new Tuplet( name, d, s, o, 0 ) );
 585                         break;
 586                     case 11:
 587                         nv.add( new Tuplet( name, o, sn, d, 2 ) );
 588                         break;
 589                     case 12:
 590                         nv.add( new Tuplet( name, o, sn, s, d, 3 ) );
 591                         break;
 592                     case 13:
 593                         nv.add( new Tuplet( name, o, s, sn, d, 3 ) );
 594                         break;
 595                     case 14:
 596                         nv.add( new Tuplet( name, sn, o, s, d, 3 ) );
 597                         break;
 598                     case 15:
 599                         nv.add( new Tuplet( name, sn, d, o, s, 1 ) );
 600                         break;
 601                     case 16:
 602                         nv.add( new Tuplet( name, o, d, sn, s, 1 ) );
 603                         break;
 604                     case 17:
 605                         nv.add( new Tuplet( name, o, sn, d, s, 2 ) );
 606                         break;
 607                     default:
 608                         nv.add( new Tuplet( name, o, 5 ) );
 609                 }
 610             } // End of I for loop.
 611             if ( ( nv != null ) && ( nv.size() > 0 ) ) {
 612                 orthologs.append( "[seq name]\t\t[ortho]\t[st-n]\t[sup-o]\t[dist]" + ForesterUtil.LINE_SEPARATOR );
 613                 final Tuplet[] nv_array = new Tuplet[ nv.size() ];
 614                 for( int j = 0; j < nv.size(); ++j ) {
 615                     nv_array[ j ] = nv.get( j );
 616                 }
 617                 Arrays.sort( nv_array );
 618                 for( final Tuplet element : nv_array ) {
 619                     name = element.getKey();
 620                     value1 = element.getValue1();
 621                     value2 = element.getValue2();
 622                     value3 = element.getValue3();
 623                     value4 = element.getValue4();
 624                     orthologs.append( addNameAndValues( name, value1, value2, value3, value4, sort ) );
 625                 }
 626             }
 627         }
 628         // No orthologs found.
 629         if ( ( orthologs == null ) || ( orthologs.length() < 1 ) ) {
 630             orthologs.append( "-" );
 631         }
 632         return orthologs;
 633     } // inferredOrthologsToString( String, int, double )
 634
 635     /**
 636      * Returns a String containg the names of orthologs of the PhylogenyNode
 637      * with seq name query_name. The String also contains how many times a
 638      * particular ortholog has been observed. Returns "-" if no putative
 639      * orthologs have been found (given threshold_orthologs).
 640      * <p>
 641      * Orthologs are to be inferred by method "inferOrthologs".
 642      *
 643      * @param query_name
 644      *            sequence name of a external node of the gene trees
 645      * @param return_dists
 646      * @param threshold_ultra_paralogs
 647      *            between 1 and 100
 648      * @return String containing the inferred orthologs, String containing "-"
 649      *         if no orthologs have been found null in case of error
 650      */
 651     public String inferredUltraParalogsToString( final String query_name,
 652                                                  final boolean return_dists,
 653                                                  double threshold_ultra_paralogs ) {
 654         HashMap<String, Integer> sp_hashmap = null;
 655         String name = "", ultra_paralogs = "";
 656         int sort = 0;
 657         double sp = 0.0, value1 = 0.0, value2 = 0.0, d = 0.0;
 658         final List<Tuplet> nv = new ArrayList<Tuplet>();
 659         if ( threshold_ultra_paralogs < 1.0 ) {
 660             threshold_ultra_paralogs = 1.0;
 661         }
 662         else if ( threshold_ultra_paralogs > 100.0 ) {
 663             threshold_ultra_paralogs = 100.0;
 664         }
 665         if ( _up_hash_maps == null ) {
 666             throw new RuntimeException( "Ultra paralogs have not been calculated (successfully)." );
 667         }
 668         if ( return_dists && ( _m == null ) && ( _l == null ) ) {
 669             throw new RuntimeException( "Distance list or matrix have not been read in (successfully)." );
 670         }
 671         sp_hashmap = getInferredUltraParalogs( query_name );
 672         if ( sp_hashmap == null ) {
 673             throw new RuntimeException( "Ultra paralogs for " + query_name + " were not established" );
 674         }
 675         if ( _seq_names.size() > 0 ) {
 676             I: for( int i = 0; i < _seq_names.size(); ++i ) {
 677                 name = _seq_names.get( i );
 678                 if ( name.equals( query_name ) ) {
 679                     continue I;
 680                 }
 681                 sp = getBootstrapValueFromHash( sp_hashmap, name );
 682                 if ( sp < threshold_ultra_paralogs ) {
 683                     continue I;
 684                 }
 685                 if ( return_dists ) {
 686                     if ( _m != null ) {
 687                         d = getDistance( query_name, name );
 688                     }
 689                     else {
 690                         d = getDistance( name );
 691                     }
 692                     nv.add( new Tuplet( name, sp, d, 1 ) );
 693                 }
 694                 else {
 695                     nv.add( new Tuplet( name, sp, 5 ) );
 696                 }
 697             } // End of I for loop.
 698             if ( ( nv != null ) && ( nv.size() > 0 ) ) {
 699                 final Tuplet[] nv_array = new Tuplet[ nv.size() ];
 700                 for( int j = 0; j < nv.size(); ++j ) {
 701                     nv_array[ j ] = nv.get( j );
 702                 }
 703                 Arrays.sort( nv_array );
 704                 if ( return_dists ) {
 705                     sort = 91;
 706                 }
 707                 else {
 708                     sort = 90;
 709                 }
 710                 for( final Tuplet element : nv_array ) {
 711                     name = element.getKey();
 712                     value1 = element.getValue1();
 713                     value2 = element.getValue2();
 714                     ultra_paralogs += addNameAndValues( name, value1, value2, 0.0, 0.0, sort );
 715                 }
 716             }
 717         }
 718         // No ultra paralogs found.
 719         if ( ( ultra_paralogs == null ) || ( ultra_paralogs.length() < 1 ) ) {
 720             ultra_paralogs = "-";
 721         }
 722         return ultra_paralogs;
 723     }
 724
 725     public final void readDistanceMatrix( final File matrix_file ) throws IOException {
 726         DistanceMatrix[] matrices = null;
 727         final SymmetricalDistanceMatrixParser parser = SymmetricalDistanceMatrixParser.createInstance();
 728         matrices = parser.parse( matrix_file );
 729         if ( ( matrices == null ) || ( matrices.length == 0 ) ) {
 730             throw new IOException( "failed to parse distance matrix from [" + matrix_file + "]" );
 731         }
 732         if ( matrices.length > 1 ) {
 733             throw new IOException( "[" + matrix_file + "] contains more than once distance matrix" );
 734         }
 735         _m = matrices[ 0 ];
 736     }
 737
 738     /**
 739      * Brings this into the same state as immediately after construction.
 740      */
 741     private final void reset() {
 742         _o_hash_maps = null;
 743         _so_hash_maps = null;
 744         _up_hash_maps = null;
 745         _seq_names = null;
 746         _m = null;
 747         _l = null;
 748         _bootstraps = 1;
 749         _ext_nodes_ = 0;
 750         _time = 0;
 751     }
 752
 753     /**
 754      * Sets the numbers of trees analyzed.
 755      * @param the
 756      *            numbers of trees analyzed
 757      */
 758     private void setBootstraps( int i ) {
 759         if ( i < 1 ) {
 760             i = 1;
 761         }
 762         _bootstraps = i;
 763     }
 764
 765     /**
 766      * Sets number of ext nodes in gene trees analyzed (after stripping).
 767      * @param the
 768      *            number of ext nodes in gene trees analyzed (after stripping)
 769      */
 770     private void setExtNodesOfAnalyzedGeneTrees( int i ) {
 771         if ( i < 1 ) {
 772             i = 0;
 773         }
 774         _ext_nodes_ = i;
 775     }
 776
 777     // Helper for doInferOrthologs( Phylogeny, Phylogeny, String )
 778     // and doInferOrthologs( Phylogeny, Phylogeny ).
 779     private void updateHash( final HashMap<String, HashMap<String, Integer>> counter_map,
 780                              final String query_seq_name,
 781                              final List<PhylogenyNode> nodes ) {
 782         final HashMap<String, Integer> hash_map = counter_map.get( query_seq_name );
 783         if ( hash_map == null ) {
 784             throw new RuntimeException( "Unexpected failure in method updateHash." );
 785         }
 786         for( int j = 0; j < nodes.size(); ++j ) {
 787             String seq_name;
 788             if ( ( nodes.get( j ) ).getNodeData().isHasSequence()
 789                     && !ForesterUtil.isEmpty( ( nodes.get( j ) ).getNodeData().getSequence().getName() ) ) {
 790                 seq_name = ( nodes.get( j ) ).getNodeData().getSequence().getName();
 791             }
 792             else {
 793                 seq_name = ( nodes.get( j ) ).getName();
 794             }
 795             if ( hash_map.containsKey( seq_name ) ) {
 796                 hash_map.put( seq_name, hash_map.get( seq_name ) + 1 );
 797             }
 798             else {
 799                 hash_map.put( seq_name, 1 );
 800             }
 801         }
 802     }
 803
 804     // Helper method for inferredOrthologsToString
 805     // and inferredUltraParalogsToString.
 806     private final static String addNameAndValues( final String name,
 807                                                   final double value1,
 808                                                   final double value2,
 809                                                   final double value3,
 810                                                   final double value4,
 811                                                   final int sort ) {
 812         final java.text.DecimalFormat df = new java.text.DecimalFormat( "0.#####" );
 813         df.setDecimalSeparatorAlwaysShown( false );
 814         String line = "";
 815         if ( name.length() < 8 ) {
 816             line += ( name + "\t\t\t" );
 817         }
 818         else if ( name.length() < 16 ) {
 819             line += ( name + "\t\t" );
 820         }
 821         else {
 822             line += ( name + "\t" );
 823         }
 824         switch ( sort ) {
 825             case 0:
 826                 line += addToLine( value1, df );
 827                 line += "-\t";
 828                 line += "-\t";
 829                 line += "-\t";
 830                 break;
 831             case 1:
 832                 line += addToLine( value1, df );
 833                 line += "-\t";
 834                 line += addToLine( value2, df );
 835                 line += "-\t";
 836                 break;
 837             case 2:
 838                 line += addToLine( value2, df );
 839                 line += "-\t";
 840                 line += addToLine( value1, df );
 841                 line += "-\t";
 842                 break;
 843             case 3:
 844                 line += addToLine( value1, df );
 845                 line += "-\t";
 846                 line += "-\t";
 847                 line += addToLine( value2, df );
 848                 break;
 849             case 4:
 850                 line += addToLine( value2, df );
 851                 line += "-\t";
 852                 line += "-\t";
 853                 line += addToLine( value1, df );
 854                 break;
 855             case 5:
 856                 line += addToLine( value1, df );
 857                 line += "-\t";
 858                 line += addToLine( value2, df );
 859                 line += addToLine( value3, df );
 860                 break;
 861             case 6:
 862                 line += addToLine( value1, df );
 863                 line += "-\t";
 864                 line += addToLine( value3, df );
 865                 line += addToLine( value2, df );
 866                 break;
 867             case 7:
 868                 line += addToLine( value2, df );
 869                 line += "-\t";
 870                 line += addToLine( value1, df );
 871                 line += addToLine( value3, df );
 872                 break;
 873             case 8:
 874                 line += addToLine( value3, df );
 875                 line += "-\t";
 876                 line += addToLine( value1, df );
 877                 line += addToLine( value2, df );
 878                 break;
 879             case 9:
 880                 line += addToLine( value2, df );
 881                 line += "-\t";
 882                 line += addToLine( value3, df );
 883                 line += addToLine( value1, df );
 884                 break;
 885             case 10:
 886                 line += addToLine( value3, df );
 887                 line += "-\t";
 888                 line += addToLine( value2, df );
 889                 line += addToLine( value1, df );
 890                 break;
 891             case 11:
 892                 line += addToLine( value1, df );
 893                 line += addToLine( value2, df );
 894                 line += "-\t";
 895                 line += addToLine( value3, df );
 896                 break;
 897             case 12:
 898                 line += addToLine( value1, df );
 899                 line += addToLine( value2, df );
 900                 line += addToLine( value3, df );
 901                 line += addToLine( value4, df );
 902                 break;
 903             case 13:
 904                 line += addToLine( value1, df );
 905                 line += addToLine( value3, df );
 906                 line += addToLine( value2, df );
 907                 line += addToLine( value4, df );
 908                 break;
 909             case 14:
 910                 line += addToLine( value2, df );
 911                 line += addToLine( value1, df );
 912                 line += addToLine( value3, df );
 913                 line += addToLine( value4, df );
 914                 break;
 915             case 15:
 916                 line += addToLine( value3, df );
 917                 line += addToLine( value1, df );
 918                 line += addToLine( value4, df );
 919                 line += addToLine( value2, df );
 920                 break;
 921             case 16:
 922                 line += addToLine( value1, df );
 923                 line += addToLine( value3, df );
 924                 line += addToLine( value4, df );
 925                 line += addToLine( value2, df );
 926                 break;
 927             case 17:
 928                 line += addToLine( value1, df );
 929                 line += addToLine( value2, df );
 930                 line += addToLine( value4, df );
 931                 line += addToLine( value3, df );
 932                 break;
 933             case 90:
 934                 line += addToLine( value1, df );
 935                 line += "-\t";
 936                 break;
 937             case 91:
 938                 line += addToLine( value1, df );
 939                 line += addToLine( value2, df );
 940                 break;
 941         }
 942         line += ForesterUtil.LINE_SEPARATOR;
 943         return line;
 944     }
 945
 946     // Helper for addNameAndValues.
 947     private final static String addToLine( final double value, final java.text.DecimalFormat df ) {
 948         String s = "";
 949         if ( value != Tuplet.DEFAULT ) {
 950             s = df.format( value ) + "\t";
 951         }
 952         else {
 953             s = "-\t";
 954         }
 955         return s;
 956     }
 957
 958     private static List<String> getAllExternalSequenceNames( final Phylogeny phy ) {
 959         final List<String> names = new ArrayList<String>();
 960         for( final PhylogenyNodeIterator iter = phy.iteratorExternalForward(); iter.hasNext(); ) {
 961             final PhylogenyNode n = iter.next();
 962             if ( n.getNodeData().isHasSequence() && !ForesterUtil.isEmpty( n.getNodeData().getSequence().getName() ) ) {
 963                 names.add( n.getNodeData().getSequence().getName() );
 964             }
 965             else if ( !ForesterUtil.isEmpty( n.getName() ) ) {
 966                 names.add( n.getName() );
 967             }
 968             else {
 969                 throw new IllegalArgumentException( "node has no (sequence) name: " + n );
 970             }
 971         }
 972         return names;
 973     }
 974
 975     /**
 976      * Returns the order in which ortholog (o), "super ortholog" (s) and
 977      * distance (d) are returned and sorted (priority of sort always goes from
 978      * left to right), given sort. For the meaning of sort
 979      *
 980      * @see #inferredOrthologsToString(String,int,double,double)
 981      *
 982      * @param sort
 983      *            determines order and sort priority
 984      * @return String indicating the order
 985      */
 986     public final static String getOrder( final int sort ) {
 987         String order = "";
 988         switch ( sort ) {
 989             case 0:
 990                 order = "orthologies";
 991                 break;
 992             case 1:
 993                 order = "orthologies > super orthologies";
 994                 break;
 995             case 2:
 996                 order = "super orthologies > orthologies";
 997                 break;
 998             case 3:
 999                 order = "orthologies > distance to query";
1000                 break;
1001             case 4:
1002                 order = "distance to query > orthologies";
1003                 break;
1004             case 5:
1005                 order = "orthologies > super orthologies > distance to query";
1006                 break;
1007             case 6:
1008                 order = "orthologies > distance to query > super orthologies";
1009                 break;
1010             case 7:
1011                 order = "super orthologies > orthologies > distance to query";
1012                 break;
1013             case 8:
1014                 order = "super orthologies > distance to query > orthologies";
1015                 break;
1016             case 9:
1017                 order = "distance to query > orthologies > super orthologies";
1018                 break;
1019             case 10:
1020                 order = "distance to query > super orthologies > orthologies";
1021                 break;
1022             case 11:
1023                 order = "orthologies > subtree neighbors > distance to query";
1024                 break;
1025             case 12:
1026                 order = "orthologies > subtree neighbors > super orthologies > distance to query";
1027                 break;
1028             case 13:
1029                 order = "orthologies > super orthologies > subtree neighbors > distance to query";
1030                 break;
1031             case 14:
1032                 order = "subtree neighbors > orthologies > super orthologies > distance to query";
1033                 break;
1034             case 15:
1035                 order = "subtree neighbors > distance to query > orthologies > super orthologies";
1036                 break;
1037             case 16:
1038                 order = "orthologies > distance to query > subtree neighbors > super orthologies";
1039                 break;
1040             case 17:
1041                 order = "orthologies > subtree neighbors > distance to query > super orthologies";
1042                 break;
1043             default:
1044                 order = "orthologies";
1045                 break;
1046         }
1047         return order;
1048     }
1049
1050     public final static StringBuffer getOrderHelp() {
1051         final StringBuffer sb = new StringBuffer();
1052         sb.append( "  0: orthologies" + ForesterUtil.LINE_SEPARATOR );
1053         sb.append( "  1: orthologies > super orthologies" + ForesterUtil.LINE_SEPARATOR );
1054         sb.append( "  2: super orthologies > orthologies" + ForesterUtil.LINE_SEPARATOR );
1055         sb.append( "  3: orthologies > distance to query" + ForesterUtil.LINE_SEPARATOR );
1056         sb.append( "  4: distance to query > orthologies" + ForesterUtil.LINE_SEPARATOR );
1057         sb.append( "  5: orthologies > super orthologies > distance to query" + ForesterUtil.LINE_SEPARATOR );
1058         sb.append( "  6: orthologies > distance to query > super orthologies" + ForesterUtil.LINE_SEPARATOR );
1059         sb.append( "  7: super orthologies > orthologies > distance to query" + ForesterUtil.LINE_SEPARATOR );
1060         sb.append( "  8: super orthologies > distance to query > orthologies" + ForesterUtil.LINE_SEPARATOR );
1061         sb.append( "  9: distance to query > orthologies > super orthologies" + ForesterUtil.LINE_SEPARATOR );
1062         sb.append( " 10: distance to query > super orthologies > orthologies" + ForesterUtil.LINE_SEPARATOR );
1063         sb.append( " 11: orthologies > subtree neighbors > distance to query" + ForesterUtil.LINE_SEPARATOR );
1064         sb.append( " 12: orthologies > subtree neighbors > super orthologies > distance to query"
1065                 + ForesterUtil.LINE_SEPARATOR );
1066         sb.append( " 13: orthologies > super orthologies > subtree neighbors > distance to query"
1067                 + ForesterUtil.LINE_SEPARATOR );
1068         sb.append( " 14: subtree neighbors > orthologies > super orthologies > distance to query"
1069                 + ForesterUtil.LINE_SEPARATOR );
1070         sb.append( " 15: subtree neighbors > distance to query > orthologies > super orthologies"
1071                 + ForesterUtil.LINE_SEPARATOR );
1072         sb.append( " 16: orthologies > distance to query > subtree neighbors > super orthologies"
1073                 + ForesterUtil.LINE_SEPARATOR );
1074         sb.append( " 17: orthologies > subtree neighbors > distance to query > super orthologies"
1075                 + ForesterUtil.LINE_SEPARATOR );
1076         return sb;
1077     }
1078
1079     private final static List<PhylogenyNode> getSubtreeNeighbors( final PhylogenyNode query, final int level ) {
1080         PhylogenyNode node = query;
1081         if ( !node.isExternal() ) {
1082             return null;
1083         }
1084         if ( !node.isRoot() ) {
1085             node = node.getParent();
1086         }
1087         if ( level == 2 ) {
1088             if ( !node.isRoot() ) {
1089                 node = node.getParent();
1090             }
1091         }
1092         else {
1093             throw new IllegalArgumentException( "currently only supporting level 2 subtree neighbors " );
1094         }
1095         final List<PhylogenyNode> sn = node.getAllExternalDescendants();
1096         sn.remove( query );
1097         return sn;
1098     }
1099 }