2 // FORESTER -- software libraries and applications
3 // for evolutionary biology research and applications.
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
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.
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.
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
25 // Contact: phylosoft @ gmail . com
26 // WWW: www.phylosoft.org/forester
28 package org.forester.application;
31 import java.io.IOException;
32 import java.util.ArrayList;
33 import java.util.List;
35 import org.forester.datastructures.IntMatrix;
36 import org.forester.io.parsers.phyloxml.PhyloXmlParser;
37 import org.forester.phylogeny.Phylogeny;
38 import org.forester.phylogeny.PhylogenyMethods;
39 import org.forester.phylogeny.factories.ParserBasedPhylogenyFactory;
40 import org.forester.phylogeny.factories.PhylogenyFactory;
41 import org.forester.rio.RIO;
42 import org.forester.rio.RIO.REROOTING;
43 import org.forester.rio.RIOException;
44 import org.forester.sdi.SDIException;
45 import org.forester.sdi.SDIutil.ALGORITHM;
46 import org.forester.util.BasicDescriptiveStatistics;
47 import org.forester.util.CommandLineArguments;
48 import org.forester.util.EasyWriter;
49 import org.forester.util.ForesterUtil;
53 final static private String PRG_NAME = "rio";
54 final static private String PRG_VERSION = "4.000 beta 3";
55 final static private String PRG_DATE = "2012.12.17";
56 final static private String E_MAIL = "czmasek@burnham.org";
57 final static private String WWW = "www.phylosoft.org/forester/";
58 final static private String HELP_OPTION_1 = "help";
59 final static private String HELP_OPTION_2 = "h";
60 final static private String USE_SDIR = "b";
62 public static void main( final String[] args ) {
63 ForesterUtil.printProgramInformation( PRG_NAME,
64 "resampled inference of orthologs",
69 ForesterUtil.getForesterLibraryInformation() );
70 CommandLineArguments cla = null;
72 cla = new CommandLineArguments( args );
74 catch ( final Exception e ) {
75 ForesterUtil.fatalError( PRG_NAME, e.getMessage() );
77 if ( cla.isOptionSet( HELP_OPTION_1 ) || cla.isOptionSet( HELP_OPTION_2 ) || ( args.length == 0 ) ) {
80 if ( ( args.length < 3 ) || ( args.length > 8 ) ) {
82 System.out.println( "[" + PRG_NAME + "] incorrect number of arguments" );
86 final List<String> allowed_options = new ArrayList<String>();
87 allowed_options.add( USE_SDIR );
88 final String dissallowed_options = cla.validateAllowedOptionsAsString( allowed_options );
89 if ( dissallowed_options.length() > 0 ) {
90 ForesterUtil.fatalError( PRG_NAME, "unknown option(s): " + dissallowed_options );
92 final File gene_trees_file = cla.getFile( 0 );
93 final File species_tree_file = cla.getFile( 1 );
94 final File orthology_outtable = cla.getFile( 2 );
96 if ( cla.getNumberOfNames() > 3 ) {
97 logfile = cla.getFile( 3 );
98 if ( logfile.exists() ) {
99 ForesterUtil.fatalError( PRG_NAME, "\"" + logfile + "\" already exists" );
105 final String outgroup = "";
106 ForesterUtil.fatalErrorIfFileNotReadable( PRG_NAME, gene_trees_file );
107 ForesterUtil.fatalErrorIfFileNotReadable( PRG_NAME, species_tree_file );
108 if ( orthology_outtable.exists() ) {
109 ForesterUtil.fatalError( PRG_NAME, "\"" + orthology_outtable + "\" already exists" );
111 boolean sdir = false;
112 if ( cla.isOptionSet( USE_SDIR ) ) {
114 if ( logfile != null ) {
115 ForesterUtil.fatalError( PRG_NAME, "logfile output only for GSDIR algorithm" );
119 System.out.println( "Gene trees : " + gene_trees_file );
120 System.out.println( "Species tree : " + species_tree_file );
121 System.out.println( "All vs all orthology table: " + orthology_outtable );
123 if ( logfile != null ) {
124 System.out.println( "Logfile : " + logfile );
126 System.out.println( "Non binary species tree : allowed (GSDIR algorithm)" );
129 System.out.println( "Non binary species tree : disallowed (SDIR algorithm)" );
131 time = System.currentTimeMillis();
132 Phylogeny species_tree = null;
134 final PhylogenyFactory factory = ParserBasedPhylogenyFactory.getInstance();
135 species_tree = factory.create( species_tree_file, new PhyloXmlParser() )[ 0 ];
137 catch ( final Exception e ) {
141 if ( !species_tree.isRooted() ) {
142 ForesterUtil.fatalError( PRG_NAME, "species tree is not rooted" );
144 final int o = PhylogenyMethods.countNumberOfOneDescendantNodes( species_tree );
146 ForesterUtil.printWarningMessage( PRG_NAME, "species tree has " + o
147 + " internal nodes with only one descendent! Going to strip them." );
148 PhylogenyMethods.deleteInternalNodesWithOnlyOneDescendent( species_tree );
149 if ( PhylogenyMethods.countNumberOfOneDescendantNodes( species_tree ) > 0 ) {
150 ForesterUtil.unexpectedFatalError( PRG_NAME, "stripping of one-desc nodes failed" );
153 final ALGORITHM algorithm;
155 algorithm = ALGORITHM.SDIR;
158 algorithm = ALGORITHM.GSDIR;
161 final RIO rio = RIO.executeAnalysis( gene_trees_file,
164 REROOTING.BY_ALGORITHM,
168 if ( algorithm == ALGORITHM.GSDIR ) {
169 ForesterUtil.programMessage( PRG_NAME, "taxonomy linking based on: " + rio.getGSDIRtaxCompBase() );
171 tableOutput( orthology_outtable, rio );
172 if ( ( algorithm != ALGORITHM.SDIR ) && ( logfile != null ) ) {
173 writeLogFile( logfile,
181 ForesterUtil.getForesterLibraryInformation() );
183 final BasicDescriptiveStatistics stats = rio.getDuplicationsStatistics();
184 final java.text.DecimalFormat df = new java.text.DecimalFormat( "0.#" );
185 ForesterUtil.programMessage( PRG_NAME,
186 "Mean number of duplications : " + df.format( stats.arithmeticMean() )
187 + " (sd: " + df.format( stats.sampleStandardDeviation() ) + ")" );
188 if ( stats.getN() > 3 ) {
189 ForesterUtil.programMessage( PRG_NAME, "Median number of duplications: " + df.format( stats.median() ) );
191 ForesterUtil.programMessage( PRG_NAME, "Minimum duplications : " + ( int ) stats.getMin() );
192 ForesterUtil.programMessage( PRG_NAME, "Maximum duplications : " + ( int ) stats.getMax() );
194 catch ( final RIOException e ) {
195 ForesterUtil.fatalError( PRG_NAME, e.getLocalizedMessage() );
197 catch ( final SDIException e ) {
198 ForesterUtil.fatalError( PRG_NAME, e.getLocalizedMessage() );
200 catch ( final IOException e ) {
201 ForesterUtil.fatalError( PRG_NAME, e.getLocalizedMessage() );
203 catch ( final Exception e ) {
204 ForesterUtil.unexpectedFatalError( PRG_NAME, e );
206 time = System.currentTimeMillis() - time;
207 ForesterUtil.programMessage( PRG_NAME, "time: " + time + "ms" );
208 ForesterUtil.programMessage( PRG_NAME, "OK" );
212 private final static void printHelp() {
213 System.out.println( "Usage" );
214 System.out.println();
217 + " [options] <gene trees infile> <species tree infile> <all vs all orthology table outfile> [logfile]" );
218 System.out.println();
219 System.out.println( " Options" );
220 System.out.println( " -" + USE_SDIR
221 + " : to use SDIR instead of GSDIR (faster, but non-binary species trees are disallowed)" );
222 System.out.println();
223 System.out.println( " Formats" );
224 System.out.println( " The species tree is expected to be in phyloXML format." );
226 .println( " The gene trees ideally are in phyloXML as well, but can also be in New Hamphshire (Newick)" );
227 System.out.println( " or Nexus format as long as species information can be extracted from the gene names" );
228 System.out.println( " (e.g. \"HUMAN\" from \"BCL2_HUMAN\")." );
229 System.out.println();
230 System.out.println( " Examples" );
231 System.out.println( " \"rio gene_trees.nh species.xml outtable.tsv log.txt\"" );
232 System.out.println();
233 System.out.println( " More information: http://code.google.com/p/forester/wiki/RIO" );
234 System.out.println();
238 private static void tableOutput( final File table_outfile, final RIO rio ) throws IOException, RIOException {
239 final IntMatrix m = RIO.calculateOrthologTable( rio.getAnalyzedGeneTrees(), true );
240 writeTable( table_outfile, rio, m );
243 private static void writeLogFile( final File logfile,
245 final File species_tree_file,
246 final File gene_trees_file,
248 final String prg_name,
250 final String prg_date,
251 final String f ) throws IOException {
252 final EasyWriter out = ForesterUtil.createEasyWriter( logfile );
253 out.println( prg_name );
254 out.println( "version : " + prg_v );
255 out.println( "date : " + prg_date );
256 out.println( "based on: " + f );
257 out.println( "----------------------------------" );
258 out.println( "Gene trees : " + gene_trees_file );
259 out.println( "Species tree : " + species_tree_file );
260 out.println( "All vs all orthology table : " + outtable );
262 out.println( rio.getLog().toString() );
264 ForesterUtil.programMessage( PRG_NAME, "wrote log to \"" + logfile + "\"" );
267 private static void writeTable( final File table_outfile, final RIO rio, final IntMatrix m ) throws IOException {
268 final EasyWriter w = ForesterUtil.createEasyWriter( table_outfile );
269 final java.text.DecimalFormat df = new java.text.DecimalFormat( "0.###" );
270 df.setDecimalSeparatorAlwaysShown( false );
271 for( int i = 0; i < m.size(); ++i ) {
273 w.print( m.getLabel( i ) );
276 for( int x = 0; x < m.size(); ++x ) {
277 w.print( m.getLabel( x ) );
278 for( int y = 0; y < m.size(); ++y ) {
281 if ( m.get( x, y ) != rio.getAnalyzedGeneTrees().length ) {
282 ForesterUtil.unexpectedFatalError( PRG_NAME, "diagonal value is off" );
287 w.print( df.format( ( ( double ) m.get( x, y ) ) / rio.getAnalyzedGeneTrees().length ) );
293 ForesterUtil.programMessage( PRG_NAME, "wrote table to \"" + table_outfile + "\"" );