forester/java/src/org/forester/sdi/ORcount.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.util.ArrayList;
  32 import java.util.HashMap;
  33 import java.util.List;
  34
  35 import org.forester.io.parsers.PhylogenyParser;
  36 import org.forester.io.parsers.util.ParserUtils;
  37 import org.forester.phylogeny.Phylogeny;
  38 import org.forester.phylogeny.PhylogenyMethods;
  39 import org.forester.phylogeny.PhylogenyNode;
  40 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
  41 import org.forester.phylogeny.factories.PhylogenyFactory;
  42 import org.forester.phylogeny.iterators.PhylogenyNodeIterator;
  43 import org.forester.util.ForesterUtil;
  44
  45 /*
  46  * Allows to <ul> <li> <li> <li> </ul>
  47  *
  48  * @see SDIse
  49  *
  50  * @see SDI
  51  *
  52  * @author Christian M. Zmasek
  53  *
  54  * @version 1.400 -- last modified: 10/29/2005
  55  */
  56 public class ORcount {
  57
  58     private static final String[]                     group_1              = { "ANOGA", "DROME", "CAEBR", "CAEEL" };
  59     private static final String[]                     group_2              = { "CIOIN", "FUGRU", "MOUSE", "RAT",
  60             "HUMAN"                                                       };
  61     private static final String[]                     all_species          = { "ANOGA", "DROME", "CAEBR", "CAEEL",
  62             "CIOIN", "FUGRU", "MOUSE", "RAT", "HUMAN"                     };
  63     private final Phylogeny[]                         _trees;
  64     private HashMap<String, HashMap<Object, Integer>> _species             = null;
  65     private ArrayList<String>                         _names               = null;
  66     private int                                       _group1_vs_2_counter = 0;
  67
  68     /**
  69      * Default contructor which
  70      */
  71     public ORcount( final Phylogeny[] trees ) {
  72         _trees = trees;
  73     } // ORcount( final Phylogeny[] trees )
  74
  75     private void count( final PhylogenyNode node ) {
  76         final List<PhylogenyNode> external_nodes = node.getAllExternalDescendants();
  77         for( int i = 1; i < external_nodes.size(); ++i ) {
  78             for( int j = 0; j < i; ++j ) {
  79                 final PhylogenyNode node_i = external_nodes.get( i );
  80                 final PhylogenyNode node_j = external_nodes.get( j );
  81                 final String si = PhylogenyMethods.getSpecies( node_i );
  82                 final String sj = PhylogenyMethods.getSpecies( node_j );
  83                 count( si, sj, node_i.getName(), node_j.getName() );
  84             }
  85         }
  86     } // count( PhylogenyNode )
  87
  88     private void count( final String a, final String b, final String seq_name_a, final String seq_name_b ) {
  89         HashMap<Object, Integer> h1 = _species.get( a );
  90         if ( h1 == null ) {
  91             throw new RuntimeException( "Unexpected error: Species \"" + a + "\" not present in species matrix." );
  92         }
  93         Object h2 = h1.get( b );
  94         String species_in_h1 = b;
  95         // We only look at the half matrix, and we do not know/care about the
  96         // order
  97         // of the keys (species).
  98         if ( h2 == null ) {
  99             h1 = _species.get( b );
 100             if ( h1 == null ) {
 101                 throw new RuntimeException( "Unexpected error: Species \"" + b + "\" not present in species matrix." );
 102             }
 103             h2 = h1.get( a );
 104             species_in_h1 = a;
 105         }
 106         if ( h2 == null ) {
 107             throw new RuntimeException( "Unexpected error: Species \"" + a + "\" not present in species matrix." );
 108         }
 109         h1.put( species_in_h1, new Integer( ( ( Integer ) h2 ).intValue() + 1 ) );
 110         _names.add( a + "-" + seq_name_a + " = " + b + "-" + seq_name_b );
 111     } // count( String, String )
 112
 113     public void countSharedAncestralClades( final Phylogeny tree,
 114                                             final int bootstrap_threshold,
 115                                             final String[] group_1,
 116                                             final String[] group_2 ) {
 117         if ( ( group_1 == null ) || ( group_2 == null ) ) {
 118             throw new IllegalArgumentException( "String[](s) in arguments to method \"ORcount.countSharedAncestralClades\" is (are) null." );
 119         }
 120         if ( !tree.isRooted() ) {
 121             throw new IllegalArgumentException( "Phylogeny must be rooted in order to count shared ancestral clades." );
 122         }
 123         final PhylogenyNodeIterator it = tree.iteratorPostorder();
 124         tree.setIndicatorsToZero();
 125         while ( it.hasNext() ) {
 126             final PhylogenyNode current_node = it.next();
 127             if ( current_node.getNumberOfDescendants() != 2 ) {
 128                 throw new IllegalArgumentException( "Phylogeny can not contain multifurcations in order to count shared ancestral clades." );
 129             }
 130             if ( !current_node.isExternal() ) {
 131                 final PhylogenyNode child1 = current_node.getChildNode1();
 132                 final PhylogenyNode child2 = current_node.getChildNode2();
 133                 if ( ( child1.getIndicator() == 1 ) || ( child2.getIndicator() == 1 ) ) {
 134                     current_node.setIndicator( ( byte ) 1 );
 135                 }
 136                 else {
 137                     final List<PhylogenyNode> external_nodes = current_node.getAllExternalDescendants();
 138                     final String[] external_species = new String[ external_nodes.size() ];
 139                     for( int i = 0; i < external_nodes.size(); ++i ) {
 140                         final PhylogenyNode n = external_nodes.get( i );
 141                         external_species[ i ] = PhylogenyMethods.getSpecies( n ).trim().toUpperCase();
 142                     }
 143                     if ( ForesterUtil.isIntersecting( external_species, group_1 )
 144                             && ForesterUtil.isIntersecting( external_species, group_2 ) ) {
 145                         current_node.setIndicator( ( byte ) 1 );
 146                         if ( ( group_1.length == 1 ) && ( group_2.length == 1 ) ) {
 147                             count( group_1[ 0 ], group_2[ 0 ], "name a", "name b" );
 148                         }
 149                         else {
 150                             increaseGroup1Vs2Counter();
 151                         }
 152                     }
 153                 }
 154             }
 155         } // while
 156     } // countSharedAncestralClades( Phylogeny, int )
 157
 158     public void countSharedAncestralClades( final Phylogeny[] trees, final int bootstrap_threshold ) {
 159         for( int i = 1; i < ORcount.all_species.length; ++i ) {
 160             for( int j = 0; j < i; ++j ) {
 161                 final String all_i = ORcount.all_species[ i ].trim().toUpperCase();
 162                 final String all_j = ORcount.all_species[ j ].trim().toUpperCase();
 163                 final String[] a = { all_i };
 164                 final String[] b = { all_j };
 165                 for( int k = 0; k < trees.length; ++k ) {
 166                     countSharedAncestralClades( trees[ k ], bootstrap_threshold, a, b );
 167                 }
 168             }
 169         }
 170         // print();
 171         if ( ( ORcount.group_1 != null ) && ( ORcount.group_2 != null ) && ( ORcount.group_1.length > 0 )
 172                 && ( ORcount.group_2.length > 0 ) ) {
 173             setGroup1Vs2Counter( 0 );
 174             for( int k = 0; k < trees.length; ++k ) {
 175                 countSharedAncestralClades( trees[ k ], bootstrap_threshold, ORcount.group_1, ORcount.group_2 );
 176             }
 177             System.out.println( "\nCount [(" + ForesterUtil.stringArrayToString( ORcount.group_1 ) + ") vs ("
 178                     + ForesterUtil.stringArrayToString( ORcount.group_2 ) + ")] = " + getGroup1Vs2Counter() );
 179         }
 180     }
 181
 182     public void countSuperOrthologousRelations( final int bootstrap_threshold ) {
 183         reset();
 184         for( int i = 0; i < _trees.length; ++i ) {
 185             countSuperOrthologousRelations( _trees[ i ], bootstrap_threshold );
 186         }
 187     }
 188
 189     private void countSuperOrthologousRelations( final Phylogeny tree, final int bootstrap_threshold ) {
 190         final PhylogenyNodeIterator it = tree.iteratorPostorder();
 191         if ( !tree.isRooted() ) {
 192             throw new IllegalArgumentException( "Phylogeny must be rooted in order to count 1:1 orthologous relationships." );
 193         }
 194         // The purpose of this is to find all substrees
 195         // which contain only speciation events on all their nodes.
 196         // All nodes in these subtrees are "painted" with 0's, wheres
 197         // the rest od the nodes in painted with 1's.
 198         tree.setIndicatorsToZero();
 199         it.reset();
 200         while ( it.hasNext() ) {
 201             final PhylogenyNode current_node = it.next();
 202             if ( current_node.getNumberOfDescendants() != 2 ) {
 203                 throw new IllegalArgumentException( "Phylogeny can not contain multifurcations in order to count 1:1 orthologous relationships." );
 204             }
 205             if ( !current_node.isExternal() && !current_node.isHasAssignedEvent() ) {
 206                 throw new IllegalArgumentException( "All nodes must have duplication or speciation assigned in order to count 1:1 orthologous relationships." );
 207             }
 208             if ( !current_node.isExternal()
 209                     && ( current_node.isDuplication() || ( current_node.getChildNode1().getIndicator() == 1 ) || ( current_node
 210                             .getChildNode2().getIndicator() == 1 ) ) ) {
 211                 current_node.setIndicator( ( byte ) 1 );
 212             }
 213         }
 214         // These find the largest subtrees containing only speciations
 215         // and uses their largest nodes to count all possible species
 216         // combinations
 217         // in their extant external nodes.
 218         // ~~~ this could possibly be combined with the first iteration ~~
 219         // <<<<<<<<<<<~~~~~~~~~~~~~~~<<<<<<<<<<<<<<<
 220         it.reset();
 221         while ( it.hasNext() ) {
 222             final PhylogenyNode current_node = it.next();
 223             if ( !current_node.isExternal()
 224                     && ( current_node.getIndicator() == 0 )
 225                     && ( current_node.isRoot() || ( current_node.getParent().getIndicator() == 1 ) )
 226                     && ( ( bootstrap_threshold < 1 ) || ( ( PhylogenyMethods.getConfidenceValue( current_node ) >= bootstrap_threshold )
 227                             && ( PhylogenyMethods.getConfidenceValue( current_node.getChildNode1() ) >= bootstrap_threshold ) && ( PhylogenyMethods
 228                             .getConfidenceValue( current_node.getChildNode2() ) >= bootstrap_threshold ) ) ) ) {
 229                 count( current_node );
 230             }
 231         }
 232     } // countOneToOneOrthologs( Phylogeny, int )
 233
 234     // This puts all the species found in Phylogeny array _trees into
 235     // species HashMap.
 236     private void getAllSpecies() {
 237         if ( ( getTrees() == null ) || ( getTrees().length < 1 ) ) {
 238             throw new RuntimeException( "Phylogeny array in method \"getAllSpecies( HashMap hash )\" is null or empty." );
 239         }
 240         setSpecies( new HashMap<String, HashMap<Object, Integer>>() );
 241         for( int i = 0; i < getTrees().length; ++i ) {
 242             PhylogenyNode node = getTrees()[ i ].getFirstExternalNode();
 243             while ( node != null ) {
 244                 getSpecies().put( PhylogenyMethods.getSpecies( node ), null );
 245                 node = node.getNextExternalNode();
 246             }
 247         }
 248     } // void getAllSpecies( HashMap hash )
 249
 250     private int getGroup1Vs2Counter() {
 251         return _group1_vs_2_counter;
 252     }
 253
 254     private HashMap<String, HashMap<Object, Integer>> getSpecies() {
 255         return _species;
 256     }
 257
 258     private Phylogeny[] getTrees() {
 259         return _trees;
 260     }
 261
 262     private void increaseGroup1Vs2Counter() {
 263         _group1_vs_2_counter++;
 264     }
 265
 266     private void printCount() {
 267         if ( ( _species == null ) || ( _species.size() < 2 ) ) {
 268             throw new RuntimeException( "Species HashMap in method \"setUpCountingMatrix()\" is null or contains less than two species." );
 269         }
 270         final Object[] species_array = _species.keySet().toArray();
 271         final int s = species_array.length;
 272         for( int i = 0; i < s - 1; ++i ) {
 273             final String species = ( String ) species_array[ i ];
 274             System.out.println();
 275             System.out.println( species + ":" );
 276             final HashMap<?, ?> h = _species.get( species );
 277             // Setting up HashMaps linked to by hash (=_species)
 278             // Diagonals are ignored, only half the matrix is needed.
 279             for( int j = 1 + i; j < s; ++j ) {
 280                 final String sp = ( String ) species_array[ j ];
 281                 final int c = ( ( Integer ) h.get( sp ) ).intValue();
 282                 System.out.println( species + "-" + sp + ": " + c );
 283             }
 284         }
 285     }
 286
 287     private void printNames() {
 288         for( int i = 0; i < _names.size(); ++i ) {
 289             System.out.println( i + ": " + _names.get( i ) );
 290         }
 291     }
 292
 293     public void reset() {
 294         getAllSpecies();
 295         setUpCountingMatrix();
 296         setGroup1Vs2Counter( 0 );
 297         _names = new ArrayList<String>();
 298     }
 299
 300     private void setGroup1Vs2Counter( final int group1_vs_2_counter ) {
 301         _group1_vs_2_counter = group1_vs_2_counter;
 302     }
 303
 304     private void setSpecies( final HashMap<String, HashMap<Object, Integer>> species ) {
 305         _species = species;
 306     }
 307
 308     private void setUpCountingMatrix() {
 309         if ( ( getSpecies() == null ) || ( getSpecies().size() < 2 ) ) {
 310             throw new RuntimeException( "Species HashMap in method \"setUpCountingMatrix()\" is null or contains less than two species." );
 311         }
 312         final Object[] species_array = getSpecies().keySet().toArray();
 313         final int s = species_array.length;
 314         for( int i = 0; i < s; ++i ) {
 315             final String species = ( String ) species_array[ i ];
 316             final HashMap<Object, Integer> h = new HashMap<Object, Integer>();
 317             // Setting up HashMaps linked to by hash (=_species)
 318             // Diagonals are ignored, only half the matrix is needed.
 319             for( int j = 1 + i; j < s; ++j ) {
 320                 h.put( species_array[ j ], new Integer( 0 ) );
 321             }
 322             getSpecies().put( species, h );
 323         }
 324     }
 325
 326     private static void errorInCommandLine() {
 327         System.out.println( "\nORcount: Error in command line.\n" );
 328         System.out.println( "Usage: \"\"" );
 329         System.out.println( "\nOptions:" );
 330         System.out.println( " -" );
 331         System.out.println( "" );
 332         System.exit( -1 );
 333     } // errorInCommandLine()
 334
 335     /**
 336      * Main method for this class.
 337      * <p>
 338      * (Last modified: 11/26/03)
 339      *
 340      * @param args[1or2]
 341      *            gene tree file name (in NHX format with species names in
 342      *            species name fields and sequence names in sequence name
 343      *            fields; unless -n option is used)
 344      */
 345     public static void main( final String args[] ) {
 346         if ( args.length == 0 ) {
 347             ORcount.errorInCommandLine();
 348         }
 349         final Phylogeny[] trees = new Phylogeny[ args.length ];
 350         final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
 351         for( int i = 0; i < trees.length; ++i ) {
 352             try {
 353                 System.out.println( "Reading tree #" + i + "  [" + args[ i ] + "]" );
 354                 final PhylogenyParser pp = ParserUtils.createParserDependingOnFileType( new File( args[ i ] ), true );
 355                 trees[ i ] = factory.create( new File( args[ i ] ), pp )[ 0 ];
 356             }
 357             catch ( final Exception e ) {
 358                 System.out.println( "\nFailed to read \"" + args[ i ] + "\". Terminating.\n" );
 359                 System.exit( -1 );
 360             }
 361         }
 362         System.out.println( "Finished reading in trees.\n\n" );
 363         final ORcount or_count = new ORcount( trees );
 364         try {
 365             System.out.println( "\n\n\n\"1:1 ORTHOLOGOUS GENE PAIRS\":\n" );
 366             System.out.println( "\n\n\n\"SUPER ORTHOLOGOUS GENE PAIRS\":\n" );
 367             or_count.countSuperOrthologousRelations( 0 );
 368             or_count.printNames();
 369             or_count.printCount();
 370             // System.out.println( "\n\n\n\"SHARED ANCESTRAL CLADES\":\n");
 371             // or_count.reset();
 372             // or_count.countSharedAncestralClades( trees, 0 );
 373         }
 374         catch ( final Exception e ) {
 375             System.out.println( "\nException. Terminating.\n" );
 376             System.out.println( "\nException is: " + e + "\n" );
 377             e.printStackTrace();
 378             System.exit( -1 );
 379         }
 380         System.out.println( "\nDone." );
 381         System.exit( 0 );
 382     } // main ( String )
 383 } // End of class ORcount.