forester/java/src/org/forester/sdi/RIO.java

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