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.FileWriter;
  32 import java.io.IOException;
  33 import java.io.PrintWriter;
  34 import java.util.ArrayList;
  35 import java.util.Arrays;
  36 import java.util.HashMap;
  37 import java.util.List;
  38
  39 import org.forester.evoinference.matrix.distance.DistanceMatrix;
  40 import org.forester.io.parsers.SymmetricalDistanceMatrixParser;
  41 import org.forester.io.parsers.phyloxml.PhyloXmlParser;
  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 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      */
 235     public void inferOrthologs( final File gene_trees_file, final Phylogeny species_tree, final String query )
 236             throws IOException {
 237         int bs = 0;
 238         if ( RIO.TIME ) {
 239             _time = System.currentTimeMillis();
 240         }
 241         if ( !gene_trees_file.exists() ) {
 242             throw new IllegalArgumentException( gene_trees_file.getAbsolutePath() + " does not exist." );
 243         }
 244         else if ( !gene_trees_file.isFile() ) {
 245             throw new IllegalArgumentException( gene_trees_file.getAbsolutePath() + " is not a file." );
 246         }
 247         // Read in first tree to get its sequence names
 248         // and strip species_tree.
 249         final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 250         final Phylogeny gene_tree = factory.create( gene_trees_file, new PhyloXmlParser() )[ 0 ];
 251         // Removes from species_tree all species not found in gene_tree.
 252         PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( gene_tree, species_tree );
 253         // Removes from gene_tree all species not found in species_tree.
 254         PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gene_tree );
 255         _seq_names = getAllExternalSequenceNames( gene_tree );
 256         if ( ( _seq_names == null ) || ( _seq_names.length < 1 ) ) {
 257             return;
 258         }
 259         _o_hash_maps = new HashMap<String, HashMap<String, Integer>>();
 260         _so_hash_maps = new HashMap<String, HashMap<String, Integer>>();
 261         _up_hash_maps = new HashMap<String, HashMap<String, Integer>>();
 262         _sn_hash_maps = new HashMap<String, HashMap<String, Integer>>();
 263         _o_hash_maps.put( query, new HashMap<String, Integer>( _seq_names.length ) );
 264         _so_hash_maps.put( query, new HashMap<String, Integer>( _seq_names.length ) );
 265         _up_hash_maps.put( query, new HashMap<String, Integer>( _seq_names.length ) );
 266         _sn_hash_maps.put( query, new HashMap<String, Integer>( _seq_names.length ) );
 267         // Go through all gene trees in the file.
 268         final Phylogeny[] gene_trees = factory.create( gene_trees_file, new PhyloXmlParser() );
 269         for( final Phylogeny gt : gene_trees ) {
 270             bs++;
 271             // Removes from gene_tree all species not found in species_tree.
 272             PhylogenyMethods.taxonomyBasedDeletionOfExternalNodes( species_tree, gt );
 273             inferOrthologsHelper( gt, species_tree, query );
 274             // System.out.println( bs );
 275         }
 276         setBootstraps( bs );
 277         if ( RIO.TIME ) {
 278             _time = ( System.currentTimeMillis() - _time );
 279         }
 280     }
 281
 282     // Helper method which performs the actual ortholog inference for
 283     // the external node with seqname query.
 284     private void inferOrthologsHelper( final Phylogeny gene_tree, final Phylogeny species_tree, final String query ) {
 285         Phylogeny assigned_tree = null;
 286         List<PhylogenyNode> nodes = null;
 287         final SDIR sdiunrooted = new SDIR();
 288         List<PhylogenyNode> orthologs = null;
 289         List<PhylogenyNode> super_orthologs = null;
 290         List<PhylogenyNode> ultra_paralogs = null;
 291         List<PhylogenyNode> subtree_neighbors = null;
 292         assigned_tree = sdiunrooted.infer( gene_tree,
 293                                            species_tree,
 294                                            RIO.ROOT_BY_MINIMIZING_MAPPING_COST,
 295                                            RIO.ROOT_BY_MINIMIZING_SUM_OF_DUPS,
 296                                            RIO.ROOT_BY_MINIMIZING_TREE_HEIGHT,
 297                                            true,
 298                                            1 )[ 0 ];
 299         setExtNodesOfAnalyzedGeneTrees( assigned_tree.getNumberOfExternalNodes() );
 300         nodes = assigned_tree.getNodesViaSequenceName( query );
 301         if ( nodes.size() > 1 ) {
 302             throw new IllegalArgumentException( "node named [" + query + "] not unique" );
 303         }
 304         else if ( nodes.isEmpty() ) {
 305             throw new IllegalArgumentException( "no node containing a sequence named [" + query + "] found" );
 306         }
 307         final PhylogenyNode query_node = nodes.get( 0 );
 308         final PhylogenyMethods methods = PhylogenyMethods.getInstance();
 309         orthologs = methods.getOrthologousNodes( assigned_tree, query_node );
 310         updateHash( _o_hash_maps, query, orthologs );
 311         super_orthologs = PhylogenyMethods.getSuperOrthologousNodes( query_node );
 312         updateHash( _so_hash_maps, query, super_orthologs );
 313         subtree_neighbors = getSubtreeNeighbors( query_node, 2 );
 314         updateHash( _sn_hash_maps, query, subtree_neighbors );
 315         ultra_paralogs = PhylogenyMethods.getUltraParalogousNodes( query_node );
 316         updateHash( _up_hash_maps, query, ultra_paralogs );
 317     }
 318
 319     /**
 320      * Returns an ArrayList containg the names of orthologs of the PhylogenyNode
 321      * with seq name seq_name.
 322      *
 323      * @param seq_name
 324      *            sequence name of a external node of the gene trees
 325      * @param threshold_orthologs
 326      *            the minimal number of observations for a a sequence to be
 327      *            reported as orthologous as percentage (0.0-100.0%)
 328      * @return ArrayList containg the names of orthologs of the PhylogenyNode
 329      *         with seq name seq_name
 330      */
 331     public ArrayList<String> inferredOrthologsToArrayList( final String seq_name, double threshold_orthologs ) {
 332         HashMap<String, Integer> o_hashmap = null;
 333         String name = null;
 334         double o = 0.0;
 335         final ArrayList<String> arraylist = new ArrayList<String>();
 336         if ( _o_hash_maps == null ) {
 337             throw new RuntimeException( "Orthologs have not been calculated (successfully)." );
 338         }
 339         if ( threshold_orthologs < 0.0 ) {
 340             threshold_orthologs = 0.0;
 341         }
 342         else if ( threshold_orthologs > 100.0 ) {
 343             threshold_orthologs = 100.0;
 344         }
 345         o_hashmap = getInferredOrthologs( seq_name );
 346         if ( o_hashmap == null ) {
 347             throw new RuntimeException( "Orthologs for " + seq_name + " were not established." );
 348         }
 349         if ( _seq_names.length > 0 ) {
 350             I: for( int i = 0; i < _seq_names.length; ++i ) {
 351                 name = _seq_names[ i ];
 352                 if ( name.equals( seq_name ) ) {
 353                     continue I;
 354                 }
 355                 o = getBootstrapValueFromHash( o_hashmap, name );
 356                 if ( o < threshold_orthologs ) {
 357                     continue I;
 358                 }
 359                 arraylist.add( name );
 360             }
 361         }
 362         return arraylist;
 363     }
 364
 365     /**
 366      * Returns a String containg the names of orthologs of the PhylogenyNode
 367      * with seq name query_name. The String also contains how many times a
 368      * particular ortholog has been observed.
 369      * <p>
 370      * <ul>
 371      * The output order is (per line): Name, Ortholog, Subtree neighbor, Super
 372      * ortholog, Distance
 373      * </ul>
 374      * <p>
 375      * The sort priority of this is determined by sort in the following manner:
 376      * <ul>
 377      * <li>0 : Ortholog
 378      * <li>1 : Ortholog, Super ortholog
 379      * <li>2 : Super ortholog, Ortholog
 380      * <li>3 : Ortholog, Distance
 381      * <li>4 : Distance, Ortholog
 382      * <li>5 : Ortholog, Super ortholog, Distance
 383      * <li>6 : Ortholog, Distance, Super ortholog
 384      * <li>7 : Super ortholog, Ortholog, Distance
 385      * <li>8 : Super ortholog, Distance, Ortholog
 386      * <li>9 : Distance, Ortholog, Super ortholog
 387      * <li>10 : Distance, Super ortholog, Ortholog
 388      * <li>11 : Ortholog, Subtree neighbor, Distance
 389      * <li>12 : Ortholog, Subtree neighbor, Super ortholog, Distance (default)
 390      * <li>13 : Ortholog, Super ortholog, Subtree neighbor, Distance
 391      * <li>14 : Subtree neighbor, Ortholog, Super ortholog, Distance
 392      * <li>15 : Subtree neighbor, Distance, Ortholog, Super ortholog
 393      * <li>16 : Ortholog, Distance, Subtree neighbor, Super ortholog
 394      * <li>17 : Ortholog, Subtree neighbor, Distance, Super ortholog
 395      * </ul>
 396      * <p>
 397      * Returns "-" if no putative orthologs have been found (given
 398      * threshold_orthologs).
 399      * <p>
 400      * Orthologs are to be inferred by method "inferOrthologs".
 401      * <p>
 402      * (Last modified: 05/08/01)
 403      *
 404      * @param query_name
 405      *            sequence name of a external node of the gene trees
 406      * @param sort
 407      *            order and sort priority
 408      * @param threshold_orthologs
 409      *            the minimal number of observations for a a sequence to be
 410      *            reported as orthologous, in percents (0.0-100.0%)
 411      * @param threshold_subtreeneighborings
 412      *            the minimal number of observations for a a sequence to be
 413      *            reported as orthologous, in percents (0.0-100.0%)
 414      * @return String containing the inferred orthologs, String containing "-"
 415      *         if no orthologs have been found null in case of error
 416      * @see #inferOrthologs(File,Phylogeny,String)
 417      * @see #inferOrthologs(Phylogeny[],Phylogeny)
 418      * @see #inferOrthologs(File,Phylogeny)
 419      * @see #getOrder(int)
 420      */
 421     public StringBuffer inferredOrthologsToString( final String query_name,
 422                                                    int sort,
 423                                                    double threshold_orthologs,
 424                                                    double threshold_subtreeneighborings ) {
 425         HashMap<String, Integer> o_hashmap = null;
 426         HashMap<String, Integer> s_hashmap = null;
 427         HashMap<String, Integer> n_hashmap = null;
 428         String name = "";
 429         double o = 0.0, // Orthologs.
 430         s = 0.0, // Super orthologs.
 431         sn = 0.0, // Subtree neighbors.
 432         value1 = 0.0, value2 = 0.0, value3 = 0.0, value4 = 0.0, d = 0.0;
 433         final ArrayList<Tuplet> nv = new ArrayList<Tuplet>();
 434         if ( ( _o_hash_maps == null ) || ( _so_hash_maps == null ) || ( _sn_hash_maps == null ) ) {
 435             throw new RuntimeException( "Orthologs have not been calculated (successfully)" );
 436         }
 437         if ( ( sort < 0 ) || ( sort > 17 ) ) {
 438             sort = 12;
 439         }
 440         if ( ( sort > 2 ) && ( _m == null ) && ( _l == null ) ) {
 441             throw new RuntimeException( "Distance list or matrix have not been read in (successfully)" );
 442         }
 443         if ( threshold_orthologs < 0.0 ) {
 444             threshold_orthologs = 0.0;
 445         }
 446         else if ( threshold_orthologs > 100.0 ) {
 447             threshold_orthologs = 100.0;
 448         }
 449         if ( threshold_subtreeneighborings < 0.0 ) {
 450             threshold_subtreeneighborings = 0.0;
 451         }
 452         else if ( threshold_subtreeneighborings > 100.0 ) {
 453             threshold_subtreeneighborings = 100.0;
 454         }
 455         o_hashmap = getInferredOrthologs( query_name );
 456         s_hashmap = getInferredSuperOrthologs( query_name );
 457         n_hashmap = getInferredSubtreeNeighbors( query_name );
 458         if ( ( o_hashmap == null ) || ( s_hashmap == null ) || ( n_hashmap == null ) ) {
 459             throw new RuntimeException( "Orthologs for " + query_name + " were not established" );
 460         }
 461         final StringBuffer orthologs = new StringBuffer();
 462         if ( _seq_names.length > 0 ) {
 463             I: for( int i = 0; i < _seq_names.length; ++i ) {
 464                 name = _seq_names[ i ];
 465                 if ( name.equals( query_name ) ) {
 466                     continue I;
 467                 }
 468                 o = getBootstrapValueFromHash( o_hashmap, name );
 469                 if ( o < threshold_orthologs ) {
 470                     continue I;
 471                 }
 472                 sn = getBootstrapValueFromHash( n_hashmap, name );
 473                 if ( sn < threshold_subtreeneighborings ) {
 474                     continue I;
 475                 }
 476                 s = getBootstrapValueFromHash( s_hashmap, name );
 477                 if ( sort >= 3 ) {
 478                     if ( _m != null ) {
 479                         d = getDistance( query_name, name );
 480                     }
 481                     else {
 482                         d = getDistance( name );
 483                     }
 484                 }
 485                 switch ( sort ) {
 486                     case 0:
 487                         nv.add( new Tuplet( name, o, 5 ) );
 488                         break;
 489                     case 1:
 490                         nv.add( new Tuplet( name, o, s, 5 ) );
 491                         break;
 492                     case 2:
 493                         nv.add( new Tuplet( name, s, o, 5 ) );
 494                         break;
 495                     case 3:
 496                         nv.add( new Tuplet( name, o, d, 1 ) );
 497                         break;
 498                     case 4:
 499                         nv.add( new Tuplet( name, d, o, 0 ) );
 500                         break;
 501                     case 5:
 502                         nv.add( new Tuplet( name, o, s, d, 2 ) );
 503                         break;
 504                     case 6:
 505                         nv.add( new Tuplet( name, o, d, s, 1 ) );
 506                         break;
 507                     case 7:
 508                         nv.add( new Tuplet( name, s, o, d, 2 ) );
 509                         break;
 510                     case 8:
 511                         nv.add( new Tuplet( name, s, d, o, 1 ) );
 512                         break;
 513                     case 9:
 514                         nv.add( new Tuplet( name, d, o, s, 0 ) );
 515                         break;
 516                     case 10:
 517                         nv.add( new Tuplet( name, d, s, o, 0 ) );
 518                         break;
 519                     case 11:
 520                         nv.add( new Tuplet( name, o, sn, d, 2 ) );
 521                         break;
 522                     case 12:
 523                         nv.add( new Tuplet( name, o, sn, s, d, 3 ) );
 524                         break;
 525                     case 13:
 526                         nv.add( new Tuplet( name, o, s, sn, d, 3 ) );
 527                         break;
 528                     case 14:
 529                         nv.add( new Tuplet( name, sn, o, s, d, 3 ) );
 530                         break;
 531                     case 15:
 532                         nv.add( new Tuplet( name, sn, d, o, s, 1 ) );
 533                         break;
 534                     case 16:
 535                         nv.add( new Tuplet( name, o, d, sn, s, 1 ) );
 536                         break;
 537                     case 17:
 538                         nv.add( new Tuplet( name, o, sn, d, s, 2 ) );
 539                         break;
 540                     default:
 541                         nv.add( new Tuplet( name, o, 5 ) );
 542                 }
 543             } // End of I for loop.
 544             if ( ( nv != null ) && ( nv.size() > 0 ) ) {
 545                 orthologs.append( "[seq name]\t\t[ortho]\t[st-n]\t[sup-o]\t[dist]" + ForesterUtil.LINE_SEPARATOR );
 546                 final Tuplet[] nv_array = new Tuplet[ nv.size() ];
 547                 for( int j = 0; j < nv.size(); ++j ) {
 548                     nv_array[ j ] = nv.get( j );
 549                 }
 550                 Arrays.sort( nv_array );
 551                 for( int i = 0; i < nv_array.length; ++i ) {
 552                     name = nv_array[ i ].getKey();
 553                     value1 = nv_array[ i ].getValue1();
 554                     value2 = nv_array[ i ].getValue2();
 555                     value3 = nv_array[ i ].getValue3();
 556                     value4 = nv_array[ i ].getValue4();
 557                     orthologs.append( addNameAndValues( name, value1, value2, value3, value4, sort ) );
 558                 }
 559             }
 560         }
 561         // No orthologs found.
 562         if ( ( orthologs == null ) || ( orthologs.length() < 1 ) ) {
 563             orthologs.append( "-" );
 564         }
 565         return orthologs;
 566     } // inferredOrthologsToString( String, int, double )
 567
 568     // Helper method for inferredOrthologTableToFile.
 569     // Returns individual rows for the table as String.
 570     private String inferredOrthologsToTableHelper( final String name2,
 571                                                    final String[] names,
 572                                                    final int j,
 573                                                    final boolean super_orthologs ) {
 574         HashMap<String, Integer> hashmap = null;
 575         String name = null, orthologs = new String( "" );
 576         int value = 0;
 577         if ( !super_orthologs ) {
 578             hashmap = getInferredOrthologs( name2 );
 579         }
 580         else {
 581             hashmap = getInferredSuperOrthologs( name2 );
 582         }
 583         if ( hashmap == null ) {
 584             throw new RuntimeException( "Unexpected failure in method inferredOrthologsToTableHelper" );
 585         }
 586         for( int i = 0; i < names.length; ++i ) {
 587             name = names[ i ];
 588             if ( !hashmap.containsKey( name ) ) {
 589                 value = 0;
 590             }
 591             else {
 592                 value = hashmap.get( name );
 593             }
 594             if ( i == j ) {
 595                 // Sanity check.
 596                 if ( value != 0 ) {
 597                     throw new RuntimeException( "Failed sanity check in method inferredOrthologsToTableHelper: value not 0." );
 598                 }
 599                 orthologs += ( " " + "\t" );
 600             }
 601             else {
 602                 orthologs += ( value + "\t" );
 603             }
 604         }
 605         return orthologs;
 606     }
 607
 608     /**
 609      * Writes the orthologs for each external node of the gene trees to outfile
 610      * in the form of a table. Orthologs are to be inferred by method
 611      * "inferOrthologs". Overwrites without asking! (Last modified: 12/07/00)
 612      *
 613      * @param outfile
 614      *            the File to write to
 615      */
 616     public void inferredOrthologTableToFile( final File outfile ) throws IOException {
 617         if ( _o_hash_maps == null ) {
 618             return;
 619         }
 620         inferredOrthologTableToFile( outfile, false );
 621     }
 622
 623     // Helper for inferredOrthologTableToFile(File).
 624     // (Last modified: 11/28/00)
 625     private void inferredOrthologTableToFile( final File outfile, final boolean super_orthologs ) throws IOException {
 626         String name = "", line = "";
 627         PrintWriter out = null;
 628         if ( _seq_names == null ) {
 629             throw new RuntimeException( "inferredOrthologTableToFile: seq_names_ is null." );
 630         }
 631         Arrays.sort( _seq_names );
 632         out = new PrintWriter( new FileWriter( outfile ), true );
 633         if ( out == null ) {
 634             throw new RuntimeException( "inferredOrthologTableToFile: failure to create PrintWriter." );
 635         }
 636         line = "\t\t\t\t";
 637         for( int i = 0; i < _seq_names.length; ++i ) {
 638             line += ( i + ")\t" );
 639         }
 640         line += "\n";
 641         out.println( line );
 642         for( int i = 0; i < _seq_names.length; ++i ) {
 643             name = _seq_names[ i ];
 644             if ( name.length() < 8 ) {
 645                 line = i + ")\t" + name + "\t\t\t";
 646             }
 647             else if ( name.length() < 16 ) {
 648                 line = i + ")\t" + name + "\t\t";
 649             }
 650             else {
 651                 line = i + ")\t" + name + "\t";
 652             }
 653             line += inferredOrthologsToTableHelper( name, _seq_names, i, super_orthologs );
 654             out.println( line );
 655         }
 656         out.close();
 657     }
 658
 659     /**
 660      * Writes the "super orthologs" for each external nodes of the gene trees to
 661      * outfile in the form of a table. Super orthologs are to be inferred by
 662      * method "inferOrthologs". Overwrites without asking!
 663      *
 664      * @param outfile
 665      *            the File to write to
 666      */
 667     public void inferredSuperOrthologTableToFile( final File outfile ) throws IOException {
 668         if ( _so_hash_maps == null ) {
 669             return;
 670         }
 671         inferredOrthologTableToFile( outfile, true );
 672     }
 673
 674     /**
 675      * Returns a String containg the names of orthologs of the PhylogenyNode
 676      * with seq name query_name. The String also contains how many times a
 677      * particular ortholog has been observed. Returns "-" if no putative
 678      * orthologs have been found (given threshold_orthologs).
 679      * <p>
 680      * Orthologs are to be inferred by method "inferOrthologs".
 681      *
 682      * @param query_name
 683      *            sequence name of a external node of the gene trees
 684      * @param return_dists
 685      * @param threshold_ultra_paralogs
 686      *            between 1 and 100
 687      * @return String containing the inferred orthologs, String containing "-"
 688      *         if no orthologs have been found null in case of error
 689      */
 690     public String inferredUltraParalogsToString( final String query_name,
 691                                                  final boolean return_dists,
 692                                                  double threshold_ultra_paralogs ) {
 693         HashMap<String, Integer> sp_hashmap = null;
 694         String name = "", ultra_paralogs = "";
 695         int sort = 0;
 696         double sp = 0.0, value1 = 0.0, value2 = 0.0, d = 0.0;
 697         final List<Tuplet> nv = new ArrayList<Tuplet>();
 698         if ( threshold_ultra_paralogs < 1.0 ) {
 699             threshold_ultra_paralogs = 1.0;
 700         }
 701         else if ( threshold_ultra_paralogs > 100.0 ) {
 702             threshold_ultra_paralogs = 100.0;
 703         }
 704         if ( _up_hash_maps == null ) {
 705             throw new RuntimeException( "Ultra paralogs have not been calculated (successfully)." );
 706         }
 707         if ( return_dists && ( _m == null ) && ( _l == null ) ) {
 708             throw new RuntimeException( "Distance list or matrix have not been read in (successfully)." );
 709         }
 710         sp_hashmap = getInferredUltraParalogs( query_name );
 711         if ( sp_hashmap == null ) {
 712             throw new RuntimeException( "Ultra paralogs for " + query_name + " were not established" );
 713         }
 714         if ( _seq_names.length > 0 ) {
 715             I: for( int i = 0; i < _seq_names.length; ++i ) {
 716                 name = _seq_names[ i ];
 717                 if ( name.equals( query_name ) ) {
 718                     continue I;
 719                 }
 720                 sp = getBootstrapValueFromHash( sp_hashmap, name );
 721                 if ( sp < threshold_ultra_paralogs ) {
 722                     continue I;
 723                 }
 724                 if ( return_dists ) {
 725                     if ( _m != null ) {
 726                         d = getDistance( query_name, name );
 727                     }
 728                     else {
 729                         d = getDistance( name );
 730                     }
 731                     nv.add( new Tuplet( name, sp, d, 1 ) );
 732                 }
 733                 else {
 734                     nv.add( new Tuplet( name, sp, 5 ) );
 735                 }
 736             } // End of I for loop.
 737             if ( ( nv != null ) && ( nv.size() > 0 ) ) {
 738                 final Tuplet[] nv_array = new Tuplet[ nv.size() ];
 739                 for( int j = 0; j < nv.size(); ++j ) {
 740                     nv_array[ j ] = nv.get( j );
 741                 }
 742                 Arrays.sort( nv_array );
 743                 if ( return_dists ) {
 744                     sort = 91;
 745                 }
 746                 else {
 747                     sort = 90;
 748                 }
 749                 for( int i = 0; i < nv_array.length; ++i ) {
 750                     name = nv_array[ i ].getKey();
 751                     value1 = nv_array[ i ].getValue1();
 752                     value2 = nv_array[ i ].getValue2();
 753                     ultra_paralogs += addNameAndValues( name, value1, value2, 0.0, 0.0, sort );
 754                 }
 755             }
 756         }
 757         // No ultra paralogs found.
 758         if ( ( ultra_paralogs == null ) || ( ultra_paralogs.length() < 1 ) ) {
 759             ultra_paralogs = "-";
 760         }
 761         return ultra_paralogs;
 762     }
 763
 764     public final void readDistanceMatrix( final File matrix_file ) throws IOException {
 765         DistanceMatrix[] matrices = null;
 766         final SymmetricalDistanceMatrixParser parser = SymmetricalDistanceMatrixParser.createInstance();
 767         matrices = parser.parse( matrix_file );
 768         if ( ( matrices == null ) || ( matrices.length == 0 ) ) {
 769             throw new IOException( "failed to parse distance matrix from [" + matrix_file + "]" );
 770         }
 771         if ( matrices.length > 1 ) {
 772             throw new IOException( "[" + matrix_file + "] contains more than once distance matrix" );
 773         }
 774         _m = matrices[ 0 ];
 775     }
 776
 777     /**
 778      * Brings this into the same state as immediately after construction.
 779      */
 780     private final void reset() {
 781         _o_hash_maps = null;
 782         _so_hash_maps = null;
 783         _up_hash_maps = null;
 784         _seq_names = null;
 785         _m = null;
 786         _l = null;
 787         _bootstraps = 1;
 788         _ext_nodes_ = 0;
 789         _time = 0;
 790     }
 791
 792     /**
 793      * Sets the numbers of trees analyzed.
 794      * @param the
 795      *            numbers of trees analyzed
 796      */
 797     private void setBootstraps( int i ) {
 798         if ( i < 1 ) {
 799             i = 1;
 800         }
 801         _bootstraps = i;
 802     }
 803
 804     /**
 805      * Sets number of ext nodes in gene trees analyzed (after stripping).
 806      * @param the
 807      *            number of ext nodes in gene trees analyzed (after stripping)
 808      */
 809     private void setExtNodesOfAnalyzedGeneTrees( int i ) {
 810         if ( i < 1 ) {
 811             i = 0;
 812         }
 813         _ext_nodes_ = i;
 814     }
 815
 816     // Helper for doInferOrthologs( Phylogeny, Phylogeny, String )
 817     // and doInferOrthologs( Phylogeny, Phylogeny ).
 818     private void updateHash( final HashMap<String, HashMap<String, Integer>> counter_map,
 819                              final String query_seq_name,
 820                              final List<PhylogenyNode> nodes ) {
 821         final HashMap<String, Integer> hash_map = counter_map.get( query_seq_name );
 822         if ( hash_map == null ) {
 823             throw new RuntimeException( "Unexpected failure in method updateHash." );
 824         }
 825         for( int j = 0; j < nodes.size(); ++j ) {
 826             final String seq_name = ( nodes.get( j ) ).getNodeData().getSequence().getName();
 827             if ( hash_map.containsKey( seq_name ) ) {
 828                 hash_map.put( seq_name, hash_map.get( seq_name ) + 1 );
 829             }
 830             else {
 831                 hash_map.put( seq_name, 1 );
 832             }
 833         }
 834     }
 835
 836     // Helper method for inferredOrthologsToString
 837     // and inferredUltraParalogsToString.
 838     private final static String addNameAndValues( final String name,
 839                                                   final double value1,
 840                                                   final double value2,
 841                                                   final double value3,
 842                                                   final double value4,
 843                                                   final int sort ) {
 844         final java.text.DecimalFormat df = new java.text.DecimalFormat( "0.#####" );
 845         df.setDecimalSeparatorAlwaysShown( false );
 846         String line = "";
 847         if ( name.length() < 8 ) {
 848             line += ( name + "\t\t\t" );
 849         }
 850         else if ( name.length() < 16 ) {
 851             line += ( name + "\t\t" );
 852         }
 853         else {
 854             line += ( name + "\t" );
 855         }
 856         switch ( sort ) {
 857             case 0:
 858                 line += addToLine( value1, df );
 859                 line += "-\t";
 860                 line += "-\t";
 861                 line += "-\t";
 862                 break;
 863             case 1:
 864                 line += addToLine( value1, df );
 865                 line += "-\t";
 866                 line += addToLine( value2, df );
 867                 line += "-\t";
 868                 break;
 869             case 2:
 870                 line += addToLine( value2, df );
 871                 line += "-\t";
 872                 line += addToLine( value1, df );
 873                 line += "-\t";
 874                 break;
 875             case 3:
 876                 line += addToLine( value1, df );
 877                 line += "-\t";
 878                 line += "-\t";
 879                 line += addToLine( value2, df );
 880                 break;
 881             case 4:
 882                 line += addToLine( value2, df );
 883                 line += "-\t";
 884                 line += "-\t";
 885                 line += addToLine( value1, df );
 886                 break;
 887             case 5:
 888                 line += addToLine( value1, df );
 889                 line += "-\t";
 890                 line += addToLine( value2, df );
 891                 line += addToLine( value3, df );
 892                 break;
 893             case 6:
 894                 line += addToLine( value1, df );
 895                 line += "-\t";
 896                 line += addToLine( value3, df );
 897                 line += addToLine( value2, df );
 898                 break;
 899             case 7:
 900                 line += addToLine( value2, df );
 901                 line += "-\t";
 902                 line += addToLine( value1, df );
 903                 line += addToLine( value3, df );
 904                 break;
 905             case 8:
 906                 line += addToLine( value3, df );
 907                 line += "-\t";
 908                 line += addToLine( value1, df );
 909                 line += addToLine( value2, df );
 910                 break;
 911             case 9:
 912                 line += addToLine( value2, df );
 913                 line += "-\t";
 914                 line += addToLine( value3, df );
 915                 line += addToLine( value1, df );
 916                 break;
 917             case 10:
 918                 line += addToLine( value3, df );
 919                 line += "-\t";
 920                 line += addToLine( value2, df );
 921                 line += addToLine( value1, df );
 922                 break;
 923             case 11:
 924                 line += addToLine( value1, df );
 925                 line += addToLine( value2, df );
 926                 line += "-\t";
 927                 line += addToLine( value3, df );
 928                 break;
 929             case 12:
 930                 line += addToLine( value1, df );
 931                 line += addToLine( value2, df );
 932                 line += addToLine( value3, df );
 933                 line += addToLine( value4, df );
 934                 break;
 935             case 13:
 936                 line += addToLine( value1, df );
 937                 line += addToLine( value3, df );
 938                 line += addToLine( value2, df );
 939                 line += addToLine( value4, df );
 940                 break;
 941             case 14:
 942                 line += addToLine( value2, df );
 943                 line += addToLine( value1, df );
 944                 line += addToLine( value3, df );
 945                 line += addToLine( value4, df );
 946                 break;
 947             case 15:
 948                 line += addToLine( value3, df );
 949                 line += addToLine( value1, df );
 950                 line += addToLine( value4, df );
 951                 line += addToLine( value2, df );
 952                 break;
 953             case 16:
 954                 line += addToLine( value1, df );
 955                 line += addToLine( value3, df );
 956                 line += addToLine( value4, df );
 957                 line += addToLine( value2, df );
 958                 break;
 959             case 17:
 960                 line += addToLine( value1, df );
 961                 line += addToLine( value2, df );
 962                 line += addToLine( value4, df );
 963                 line += addToLine( value3, df );
 964                 break;
 965             case 90:
 966                 line += addToLine( value1, df );
 967                 line += "-\t";
 968                 break;
 969             case 91:
 970                 line += addToLine( value1, df );
 971                 line += addToLine( value2, df );
 972                 break;
 973         }
 974         line += ForesterUtil.LINE_SEPARATOR;
 975         return line;
 976     }
 977
 978     // Helper for addNameAndValues.
 979     private final static String addToLine( final double value, final java.text.DecimalFormat df ) {
 980         String s = "";
 981         if ( value != Tuplet.DEFAULT ) {
 982             s = df.format( value ) + "\t";
 983         }
 984         else {
 985             s = "-\t";
 986         }
 987         return s;
 988     }
 989
 990     private static String[] getAllExternalSequenceNames( final Phylogeny phy ) {
 991         if ( phy.isEmpty() ) {
 992             return null;
 993         }
 994         int i = 0;
 995         final String[] names = new String[ phy.getNumberOfExternalNodes() ];
 996         for( final PhylogenyNodeIterator iter = phy.iteratorExternalForward(); iter.hasNext(); ) {
 997             names[ i++ ] = iter.next().getNodeData().getSequence().getName();
 998         }
 999         return names;
1000     }
1001
1002     /**
1003      * Returns the order in which ortholog (o), "super ortholog" (s) and
1004      * distance (d) are returned and sorted (priority of sort always goes from
1005      * left to right), given sort. For the meaning of sort
1006      *
1007      * @see #inferredOrthologsToString(String,int,double,double)
1008      *
1009      * @param sort
1010      *            determines order and sort priority
1011      * @return String indicating the order
1012      */
1013     public final static String getOrder( final int sort ) {
1014         String order = "";
1015         switch ( sort ) {
1016             case 0:
1017                 order = "orthologies";
1018                 break;
1019             case 1:
1020                 order = "orthologies > super orthologies";
1021                 break;
1022             case 2:
1023                 order = "super orthologies > orthologies";
1024                 break;
1025             case 3:
1026                 order = "orthologies > distance to query";
1027                 break;
1028             case 4:
1029                 order = "distance to query > orthologies";
1030                 break;
1031             case 5:
1032                 order = "orthologies > super orthologies > distance to query";
1033                 break;
1034             case 6:
1035                 order = "orthologies > distance to query > super orthologies";
1036                 break;
1037             case 7:
1038                 order = "super orthologies > orthologies > distance to query";
1039                 break;
1040             case 8:
1041                 order = "super orthologies > distance to query > orthologies";
1042                 break;
1043             case 9:
1044                 order = "distance to query > orthologies > super orthologies";
1045                 break;
1046             case 10:
1047                 order = "distance to query > super orthologies > orthologies";
1048                 break;
1049             case 11:
1050                 order = "orthologies > subtree neighbors > distance to query";
1051                 break;
1052             case 12:
1053                 order = "orthologies > subtree neighbors > super orthologies > distance to query";
1054                 break;
1055             case 13:
1056                 order = "orthologies > super orthologies > subtree neighbors > distance to query";
1057                 break;
1058             case 14:
1059                 order = "subtree neighbors > orthologies > super orthologies > distance to query";
1060                 break;
1061             case 15:
1062                 order = "subtree neighbors > distance to query > orthologies > super orthologies";
1063                 break;
1064             case 16:
1065                 order = "orthologies > distance to query > subtree neighbors > super orthologies";
1066                 break;
1067             case 17:
1068                 order = "orthologies > subtree neighbors > distance to query > super orthologies";
1069                 break;
1070             default:
1071                 order = "orthologies";
1072                 break;
1073         }
1074         return order;
1075     }
1076
1077     public final static StringBuffer getOrderHelp() {
1078         final StringBuffer sb = new StringBuffer();
1079         sb.append( "  0: orthologies" + ForesterUtil.LINE_SEPARATOR );
1080         sb.append( "  1: orthologies > super orthologies" + ForesterUtil.LINE_SEPARATOR );
1081         sb.append( "  2: super orthologies > orthologies" + ForesterUtil.LINE_SEPARATOR );
1082         sb.append( "  3: orthologies > distance to query" + ForesterUtil.LINE_SEPARATOR );
1083         sb.append( "  4: distance to query > orthologies" + ForesterUtil.LINE_SEPARATOR );
1084         sb.append( "  5: orthologies > super orthologies > distance to query" + ForesterUtil.LINE_SEPARATOR );
1085         sb.append( "  6: orthologies > distance to query > super orthologies" + ForesterUtil.LINE_SEPARATOR );
1086         sb.append( "  7: super orthologies > orthologies > distance to query" + ForesterUtil.LINE_SEPARATOR );
1087         sb.append( "  8: super orthologies > distance to query > orthologies" + ForesterUtil.LINE_SEPARATOR );
1088         sb.append( "  9: distance to query > orthologies > super orthologies" + ForesterUtil.LINE_SEPARATOR );
1089         sb.append( " 10: distance to query > super orthologies > orthologies" + ForesterUtil.LINE_SEPARATOR );
1090         sb.append( " 11: orthologies > subtree neighbors > distance to query" + ForesterUtil.LINE_SEPARATOR );
1091         sb.append( " 12: orthologies > subtree neighbors > super orthologies > distance to query"
1092                 + ForesterUtil.LINE_SEPARATOR );
1093         sb.append( " 13: orthologies > super orthologies > subtree neighbors > distance to query"
1094                 + ForesterUtil.LINE_SEPARATOR );
1095         sb.append( " 14: subtree neighbors > orthologies > super orthologies > distance to query"
1096                 + ForesterUtil.LINE_SEPARATOR );
1097         sb.append( " 15: subtree neighbors > distance to query > orthologies > super orthologies"
1098                 + ForesterUtil.LINE_SEPARATOR );
1099         sb.append( " 16: orthologies > distance to query > subtree neighbors > super orthologies"
1100                 + ForesterUtil.LINE_SEPARATOR );
1101         sb.append( " 17: orthologies > subtree neighbors > distance to query > super orthologies"
1102                 + ForesterUtil.LINE_SEPARATOR );
1103         return sb;
1104     }
1105
1106     private final static List<PhylogenyNode> getSubtreeNeighbors( final PhylogenyNode query, final int level ) {
1107         PhylogenyNode node = query;
1108         if ( !node.isExternal() ) {
1109             return null;
1110         }
1111         if ( !node.isRoot() ) {
1112             node = node.getParent();
1113         }
1114         if ( level == 2 ) {
1115             if ( !node.isRoot() ) {
1116                 node = node.getParent();
1117             }
1118         }
1119         else {
1120             throw new IllegalArgumentException( "currently only supporting level 2 subtree neighbors " );
1121         }
1122         final List<PhylogenyNode> sn = node.getAllExternalDescendants();
1123         sn.remove( query );
1124         return sn;
1125     }
1126 }