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