in progress
[jalview.git] / forester / java / src / org / forester / application / msa_compactor.java
index 19d5caf..e555593 100644 (file)
@@ -46,11 +46,19 @@ import org.forester.util.CommandLineArguments;
 import org.forester.util.DescriptiveStatistics;
 import org.forester.util.ForesterUtil;
 
+
+/*
+java -cp C:\Users\czmasek\SOFTWARE_DEV\ECLIPSE\forester\java\fo
+rester.jar org.forester.application.msa_compactor Bcl-2_e1_20_mafft -t
+*/
+
 public class msa_compactor {
 
-    final private static NumberFormat NF_1                                   = new DecimalFormat( "#.0" );
+    final private static NumberFormat NF_1                                   = new DecimalFormat( "0.#" );
+    final private static NumberFormat NF_4                                   = new DecimalFormat( "0.####" );
     static {
         NF_1.setRoundingMode( RoundingMode.HALF_UP );
+        NF_4.setRoundingMode( RoundingMode.HALF_UP );
     }
     final static private String       HELP_OPTION_1                          = "help";
     final static private String       HELP_OPTION_2                          = "h";
@@ -59,16 +67,17 @@ public class msa_compactor {
     final static private String       STEP_OPTION                            = "s";
     final static private String       LENGTH_OPTION                          = "l";
     final static private String       REALIGN_OPTION                         = "a";
+    final static private String       INFO_ONLY_OPTION                       = "i";
     //
     final static private String       STEP_FOR_DIAGNOSTICS_OPTION            = "sd";
     final static private String       MIN_LENGTH_OPTION                      = "ml";
     final static private String       GAP_RATIO_LENGTH_OPTION                = "gr";
     final static private String       REPORT_ENTROPY                         = "e";
-    final static private String       OUTPUT_FORMAT_PHYLIP_OPTION            = "p";
+    final static private String       OUTPUT_FORMAT_OPTION                   = "f";
     final static private String       OUTPUT_REMOVED_SEQS_OPTION             = "ro";
     final static private String       MAFFT_OPTIONS                          = "mo";
     final static private String       PERFORM_PHYLOGENETIC_INFERENCE         = "t";
-    //        
+    //
     final static private String       PATH_TO_MAFFT_OPTION                   = "mafft";
     final static private String       DO_NOT_NORMALIZE_FOR_EFF_LENGTH_OPTION = "nn";
     final static private String       PRG_NAME                               = "msa_compactor";
@@ -96,7 +105,7 @@ public class msa_compactor {
             int length = -1;
             int step = 1;
             boolean realign = false;
-            boolean norm = true;
+            boolean normalize_for_effective_seq_length = true;
             String path_to_mafft = null;
             int step_for_diagnostics = 1;
             int min_length = -1;
@@ -118,10 +127,11 @@ public class msa_compactor {
             allowed_options.add( MIN_LENGTH_OPTION );
             allowed_options.add( GAP_RATIO_LENGTH_OPTION );
             allowed_options.add( REPORT_ENTROPY );
-            allowed_options.add( OUTPUT_FORMAT_PHYLIP_OPTION );
+            allowed_options.add( OUTPUT_FORMAT_OPTION );
             allowed_options.add( OUTPUT_REMOVED_SEQS_OPTION );
             allowed_options.add( MAFFT_OPTIONS );
             allowed_options.add( PERFORM_PHYLOGENETIC_INFERENCE );
+            allowed_options.add( INFO_ONLY_OPTION );
             final String dissallowed_options = cla.validateAllowedOptionsAsString( allowed_options );
             if ( dissallowed_options.length() > 0 ) {
                 ForesterUtil.fatalError( PRG_NAME, "unknown option(s): " + dissallowed_options );
@@ -135,15 +145,20 @@ public class msa_compactor {
                 msa = DeleteableMsa.createInstance( GeneralMsaParser.parse( is ) );
             }
             final DescriptiveStatistics initial_msa_stats = MsaMethods.calculateEffectiveLengthStatistics( msa );
+            if (cla.isOptionSet( INFO_ONLY_OPTION ) ) {
+                printInfo( in, msa, initial_msa_stats );
+                System.exit( 0 );
+            }
+            
             final boolean chart_only = ( !cla.isOptionSet( LENGTH_OPTION ) )
                     && ( !cla.isOptionSet( REMOVE_WORST_OFFENDERS_OPTION ) )
-                    && ( !cla.isOptionSet( AV_GAPINESS_OPTION ) );
+                    && ( !cla.isOptionSet( AV_GAPINESS_OPTION ) && ( !cla.isOptionSet( MIN_LENGTH_OPTION ) ) );
             if ( !chart_only && ( out == null ) ) {
                 ForesterUtil.fatalError( PRG_NAME, "outfile file missing" );
             }
             if ( cla.isOptionSet( REMOVE_WORST_OFFENDERS_OPTION ) ) {
                 worst_remove = cla.getOptionValueAsInt( REMOVE_WORST_OFFENDERS_OPTION );
-                if ( ( worst_remove < 1 ) || ( worst_remove >= msa.getNumberOfSequences() - 1 ) ) {
+                if ( ( worst_remove < 1 ) || ( worst_remove >= ( msa.getNumberOfSequences() - 1 ) ) ) {
                     ForesterUtil.fatalError( PRG_NAME, "number of worst offender sequences to remove is out of range: "
                             + worst_remove );
                 }
@@ -167,7 +182,7 @@ public class msa_compactor {
                 if ( length >= msa.getLength() ) {
                     ForesterUtil.fatalError( PRG_NAME,
                                              "target length is out of range [longer than MSA (" + msa.getLength()
-                                                     + ")]: " + length );
+                                             + ")]: " + length );
                 }
                 else if ( length < initial_msa_stats.getMin() ) {
                     ForesterUtil.fatalError( PRG_NAME,
@@ -175,6 +190,22 @@ public class msa_compactor {
                                                      + initial_msa_stats.getMin() + ") ]: " + length );
                 }
             }
+            if ( cla.isOptionSet( MIN_LENGTH_OPTION ) ) {
+                if ( cla.isOptionSet( LENGTH_OPTION ) || cla.isOptionSet( REMOVE_WORST_OFFENDERS_OPTION )
+                        || cla.isOptionSet( AV_GAPINESS_OPTION ) || cla.isOptionSet( STEP_OPTION )
+                        || cla.isOptionSet( REALIGN_OPTION ) || cla.isOptionSet( PATH_TO_MAFFT_OPTION )
+                        || cla.isOptionSet( STEP_FOR_DIAGNOSTICS_OPTION ) || cla.isOptionSet( REPORT_ENTROPY )
+                        || cla.isOptionSet( OUTPUT_REMOVED_SEQS_OPTION )
+                        || cla.isOptionSet( PERFORM_PHYLOGENETIC_INFERENCE ) ) {
+                    printHelp();
+                    System.exit( 0 );
+                }
+                min_length = cla.getOptionValueAsInt( MIN_LENGTH_OPTION );
+                if ( ( min_length < 2 ) || ( min_length > initial_msa_stats.getMax() ) ) {
+                    ForesterUtil.fatalError( PRG_NAME, "value for minimal sequence length is out of range: "
+                            + min_length );
+                }
+            }
             if ( cla.isOptionSet( STEP_OPTION ) ) {
                 step = cla.getOptionValueAsInt( STEP_OPTION );
                 if ( ( step < 1 )
@@ -192,7 +223,7 @@ public class msa_compactor {
                 path_to_mafft = cla.getOptionValueAsCleanString( PATH_TO_MAFFT_OPTION );
             }
             if ( cla.isOptionSet( DO_NOT_NORMALIZE_FOR_EFF_LENGTH_OPTION ) ) {
-                norm = false;
+                normalize_for_effective_seq_length = false;
             }
             if ( cla.isOptionSet( STEP_FOR_DIAGNOSTICS_OPTION ) ) {
                 step_for_diagnostics = cla.getOptionValueAsInt( STEP_FOR_DIAGNOSTICS_OPTION );
@@ -202,13 +233,6 @@ public class msa_compactor {
                             + step_for_diagnostics );
                 }
             }
-            if ( cla.isOptionSet( MIN_LENGTH_OPTION ) ) {
-                min_length = cla.getOptionValueAsInt( MIN_LENGTH_OPTION );
-                if ( ( min_length < 2 ) || ( min_length > initial_msa_stats.getMax() ) ) {
-                    ForesterUtil.fatalError( PRG_NAME, "value for minimal sequence length is out of range: "
-                            + min_length );
-                }
-            }
             if ( cla.isOptionSet( GAP_RATIO_LENGTH_OPTION ) ) {
                 gap_ratio = cla.getOptionValueAsDouble( GAP_RATIO_LENGTH_OPTION );
                 if ( ( gap_ratio < 0 ) || ( gap_ratio > 1 ) ) {
@@ -218,8 +242,20 @@ public class msa_compactor {
             if ( cla.isOptionSet( REPORT_ENTROPY ) ) {
                 report_entropy = true;
             }
-            if ( cla.isOptionSet( OUTPUT_FORMAT_PHYLIP_OPTION ) ) {
-                output_format = MSA_FORMAT.PHYLIP;
+            if ( cla.isOptionSet( OUTPUT_FORMAT_OPTION ) ) {
+                final String fs = cla.getOptionValueAsCleanString( OUTPUT_FORMAT_OPTION );
+                if ( fs.equalsIgnoreCase( "p" ) ) {
+                    output_format = MSA_FORMAT.PHYLIP;
+                }
+                else if ( fs.equalsIgnoreCase( "f" ) ) {
+                    output_format = MSA_FORMAT.FASTA;
+                }
+                else if ( fs.equalsIgnoreCase( "n" ) ) {
+                    output_format = MSA_FORMAT.NEXUS;
+                }
+                else {
+                    ForesterUtil.fatalError( PRG_NAME, "illegal or empty output format option: " + fs );
+                }
             }
             if ( cla.isOptionSet( OUTPUT_REMOVED_SEQS_OPTION ) ) {
                 final String s = cla.getOptionValueAsCleanString( OUTPUT_REMOVED_SEQS_OPTION );
@@ -246,8 +282,8 @@ public class msa_compactor {
             if ( chart_only ) {
                 if ( ( out != null ) || ( removed_seqs_out_base != null ) ) {
                     ForesterUtil
-                            .fatalError( PRG_NAME,
-                                         "chart only, no outfile(s) produced, thus no need to indicate output file(s)" );
+                    .fatalError( PRG_NAME,
+                                 "chart only, no outfile(s) produced, thus no need to indicate output file(s)" );
                 }
                 if ( !realign && cla.isOptionSet( STEP_OPTION ) ) {
                     ForesterUtil.fatalError( PRG_NAME,
@@ -258,30 +294,14 @@ public class msa_compactor {
             if ( perform_phylogenetic_inference ) {
                 if ( step_for_diagnostics != 1 ) {
                     ForesterUtil.fatalError( PRG_NAME,
-                                             "step for diagnostics reports needs to be set to 1 for tree calculation" );
-                }
-            }
-            ForesterUtil.printProgramInformation( PRG_NAME,
-                                                  PRG_DESC,
-                                                  PRG_VERSION,
-                                                  PRG_DATE,
-                                                  E_MAIL,
-                                                  WWW,
-                                                  ForesterUtil.getForesterLibraryInformation() );
-            System.out.println( "Input MSA                            : " + in );
-            System.out.println( "  MSA length                         : " + msa.getLength() );
-            System.out.println( "  Number of sequences                : " + msa.getNumberOfSequences() );
-            System.out.println( "  Median sequence length             : " + NF_1.format( initial_msa_stats.median() ) );
-            System.out.println( "  Mean sequence length               : "
-                    + NF_1.format( initial_msa_stats.arithmeticMean() ) );
-            System.out.println( "  Max sequence length                : " + ( ( int ) initial_msa_stats.getMax() ) );
-            System.out.println( "  Min sequence length                : " + ( ( int ) initial_msa_stats.getMin() ) );
+                            "step for diagnostics reports needs to be set to 1 for tree calculation" );
+                }
+            }
+            printInfo( in, msa, initial_msa_stats );
             if ( !chart_only ) {
                 System.out.println( "Output                               : " + out );
             }
-            else {
-                System.out.println( "Output                               : n/a" );
-            }
+            
             if ( removed_seqs_out_base != null ) {
                 System.out.println( "Write removed sequences to           : " + removed_seqs_out_base );
             }
@@ -301,69 +321,86 @@ public class msa_compactor {
                 System.out.println( "Maximum allowed gap ratio per column : " + gap_ratio );
             }
             if ( ( out != null ) || ( removed_seqs_out_base != null ) ) {
-                System.out.println( "Output format                        : "
-                        + ( output_format == MSA_FORMAT.FASTA ? "fasta" : "phylip" ) );
-            }
-            if ( chart_only && !realign ) {
-                System.out.println( "Step for output and re-aligning      : n/a" );
+                System.out.print( "Output format                        : " );
+                if ( output_format == MSA_FORMAT.FASTA ) {
+                    System.out.println( "fasta" );
+                }
+                else if ( output_format == MSA_FORMAT.PHYLIP ) {
+                    System.out.println( "phylip" );
+                }
+                else if ( output_format == MSA_FORMAT.NEXUS ) {
+                    System.out.println( "nexus" );
+                }
             }
-            else {
-                if ( chart_only ) {
-                    System.out.println( "Step for re-aligning                 : " + step );
+            if ( min_length == -1 ) {
+                if ( chart_only && !realign ) {
+                    System.out.println( "Step for output and re-aligning      : n/a" );
+                }
+                else {
+                    if ( chart_only ) {
+                        System.out.println( "Step for re-aligning                 : " + step );
+                    }
+                    else {
+                        System.out.println( "Step for output and re-aligning      : " + step );
+                    }
+                }
+                System.out.println( "Step for diagnostics reports         : " + step_for_diagnostics );
+                System.out.println( "Calculate normalized Shannon Entropy : " + report_entropy );
+                if ( normalize_for_effective_seq_length ) {
+                    System.out.println( "Normalize                            : with individual, effective sequence lenghts"  );
                 }
                 else {
-                    System.out.println( "Step for output and re-aligning      : " + step );
+                    System.out.println( "Normalize                            : with MSA length" );
                 }
+                System.out.println( "Realign with MAFFT                   : " + realign );
+                if ( realign ) {
+                    System.out.println( "MAFFT options                        : " + mafft_options );
+                }
+                System.out.println( "Simple tree (Kimura distances, NJ)   : " + perform_phylogenetic_inference );
             }
-            System.out.println( "Step for diagnostics reports         : " + step_for_diagnostics );
-            System.out.println( "Calculate normalized Shannon Entropy : " + report_entropy );
-            if ( !norm ) {
-                System.out.println( "Normalize                            : " + norm );
-            }
-            System.out.println( "Realign with MAFFT                   : " + realign );
-            if ( realign ) {
-                System.out.println( "MAFFT options                        : " + mafft_options );
-            }
-            System.out.println( "Simple tree (Kimura distances, NJ)   : " + perform_phylogenetic_inference );
             System.out.println();
             final int initial_number_of_seqs = msa.getNumberOfSequences();
             List<MsaProperties> msa_props = null;
             final MsaCompactor mc = new MsaCompactor( msa );
             mc.setInfileName( in.getName() );
-            mc.setNorm( norm );
-            mc.setRealign( realign );
-            if ( realign ) {
-                mc.setPathToMafft( path_to_mafft );
-                mc.setMafftOptions( mafft_options );
-            }
-            mc.setStep( step );
-            mc.setStepForDiagnostics( step_for_diagnostics );
-            mc.setCalculateNormalizedShannonEntropy( report_entropy );
-            mc.setPeformPhylogenticInference( perform_phylogenetic_inference );
-            if ( ( worst_remove > 0 ) || ( av_gap > 0 ) || ( length > 0 ) ) {
+            if ( ( worst_remove > 0 ) || ( av_gap > 0 ) || ( length > 0 ) || ( min_length != -1 ) ) {
                 mc.setOutputFormat( output_format );
                 mc.setOutFileBase( out );
-                if ( removed_seqs_out_base != null ) {
-                    mc.setRemovedSeqsOutBase( removed_seqs_out_base );
-                }
             }
-            if ( min_length > 1 ) {
+            if ( min_length != -1 ) {
                 mc.removeSequencesByMinimalLength( min_length );
-                mc.writeMsa( new File( "removed" ) );
-            }
-            if ( worst_remove > 0 ) {
-                msa_props = mc.removeWorstOffenders( worst_remove );
-            }
-            else if ( av_gap > 0 ) {
-                msa_props = mc.removeViaGapAverage( av_gap );
-            }
-            else if ( length > 0 ) {
-                msa_props = mc.removeViaLength( length );
             }
             else {
-                msa_props = mc.chart( step, realign, norm );
+                mc.setPeformPhylogenticInference( perform_phylogenetic_inference );
+                if ( removed_seqs_out_base != null ) {
+                    mc.setRemovedSeqsOutBase( removed_seqs_out_base );
+                }
+                mc.setNorm( normalize_for_effective_seq_length );
+                mc.setRealign( realign );
+                if ( realign ) {
+                    mc.setPathToMafft( path_to_mafft );
+                    mc.setMafftOptions( mafft_options );
+                }
+                mc.setStep( step );
+                mc.setStepForDiagnostics( step_for_diagnostics );
+                mc.setCalculateNormalizedShannonEntropy( report_entropy );
+                if ( worst_remove > 0 ) {
+                    msa_props = mc.removeWorstOffenders( worst_remove );
+                }
+                else if ( av_gap > 0 ) {
+                    msa_props = mc.removeViaGapAverage( av_gap );
+                }
+                else if ( length > 0 ) {
+                    msa_props = mc.removeViaLength( length );
+                }
+                else {
+                    msa_props = mc.chart( step, realign, normalize_for_effective_seq_length );
+                }
+                Chart.display( msa_props, initial_number_of_seqs, report_entropy, in.getName() );
+                System.out.println();
+                System.out.println( "Final MSA properties" );
+                printMsaInfo( msa,  MsaMethods.calculateEffectiveLengthStatistics( msa ));
             }
-            Chart.display( msa_props, initial_number_of_seqs, report_entropy, in.getName() );
         }
         catch ( final IllegalArgumentException iae ) {
             //  iae.printStackTrace(); //TODO remove me
@@ -378,13 +415,43 @@ public class msa_compactor {
         }
     }
 
+    private static void printInfo( final File in, DeleteableMsa msa, final DescriptiveStatistics initial_msa_stats ) {
+        ForesterUtil.printProgramInformation( PRG_NAME,
+                                              PRG_DESC,
+                                              PRG_VERSION,
+                                              PRG_DATE,
+                                              E_MAIL,
+                                              WWW,
+                                              ForesterUtil.getForesterLibraryInformation() );
+        System.out.println( "Input MSA                            : " + in );
+        printMsaInfo( msa, initial_msa_stats );
+    }
+
+    private static void printMsaInfo( DeleteableMsa msa, final DescriptiveStatistics msa_stats ) {
+        System.out.println( "MSA length                           : " + msa.getLength() );
+        System.out.println( "Number of sequences                  : " + msa.getNumberOfSequences() );
+        System.out.println( "Median sequence length               : " + NF_1.format( msa_stats.median() ) );
+        System.out.println( "Mean sequence length                 : "
+                + NF_1.format( msa_stats.arithmeticMean() ) );
+        System.out.println( "Max sequence length                  : " + ( ( int ) msa_stats.getMax() ) );
+        System.out.println( "Min sequence length                  : " + ( ( int ) msa_stats.getMin() ) );
+        System.out.println( "Gap ratio                            : "
+                + NF_4.format( MsaMethods.calcGapRatio( msa ) ) );
+        System.out.println( "Mean gap count per sequence          : "
+                + NF_1.format( MsaMethods.calcNumberOfGapsStats( msa ).arithmeticMean() ) );
+        System.out.println( "Normalized Shannon Entropy (entn7)   : "
+                + NF_4.format( MsaMethods.calcNormalizedShannonsEntropy( 7, msa ) ) );
+        System.out.println( "Normalized Shannon Entropy (entn21)  : "
+                + NF_4.format( MsaMethods.calcNormalizedShannonsEntropy( 21, msa ) ) );
+    }
+
     private static void checkPathToMafft( final String path_to_mafft ) {
         if ( !ForesterUtil.isEmpty( path_to_mafft ) && MsaInferrer.isInstalled( path_to_mafft ) ) {
         }
         else {
             if ( ForesterUtil.isEmpty( path_to_mafft ) ) {
                 ForesterUtil.fatalError( PRG_NAME, "no MAFFT executable found, use -\"" + PATH_TO_MAFFT_OPTION
-                        + "=<path to MAFFT>\" option" );
+                                         + "=<path to MAFFT>\" option" );
             }
             else {
                 ForesterUtil.fatalError( PRG_NAME, "no MAFFT executable at \"" + path_to_mafft + "\"" );
@@ -410,30 +477,35 @@ public class msa_compactor {
         }
         System.out.println( "Usage:" );
         System.out.println();
-        System.out.println( PRG_NAME + " <options> <msa input file> <output file base>" );
+        System.out.println( PRG_NAME + " [options] <msa input file> [output file base]" );
         System.out.println();
         System.out.println( " options: " );
         System.out.println();
+        System.out.println( "   -" + INFO_ONLY_OPTION
+                            + "             to only display same basic information about the MSA" );
         System.out.println( "   -" + REMOVE_WORST_OFFENDERS_OPTION
-                + "=<integer>   number of worst offender sequences to remove" );
+                            + "=<integer>   number of worst offender sequences to remove" );
         System.out.println( "   -" + LENGTH_OPTION + "=<integer>   target MSA length" );
         System.out.println( "   -" + AV_GAPINESS_OPTION + "=<decimal>   target gap-ratio (0.0-1.0)" );
         System.out.println( "   -" + REALIGN_OPTION + "             to realign using MAFFT" + mafft_comment );
         System.out.println( "   -" + MAFFT_OPTIONS + "=<string>   options for MAFFT (default: --auto)" );
         System.out.println( "   -" + STEP_OPTION + "=<integer>   step for output and re-aligning (default: 1)" );
         System.out.println( "   -" + STEP_FOR_DIAGNOSTICS_OPTION
-                + "=<integer>  step for diagnostics reports (default: 1)" );
+                            + "=<integer>  step for diagnostics reports (default: 1)" );
         System.out.println( "   -" + REPORT_ENTROPY
-                + "             to calculate normalized Shannon Entropy (not recommended for very large alignments)" );
-        System.out.println( "   -" + OUTPUT_FORMAT_PHYLIP_OPTION
-                + "             to write output alignments in phylip format instead of fasta" );
+                            + "             to calculate normalized Shannon Entropy (not recommended for very large alignments)" );
+        System.out.println( "   -" + OUTPUT_FORMAT_OPTION
+                            + "=<f|p|n>     format for output alignments: f for fasta (default), p for phylip, or n for nexus" );
         System.out.println( "   -" + OUTPUT_REMOVED_SEQS_OPTION + "=<file>     to output the removed sequences" );
         System.out.println( "   -" + MIN_LENGTH_OPTION
-                + "=<integer>  minimal effecive sequence length (for deleting of shorter sequences)" );
+                            + "=<integer>  minimal effecive sequence length (for deleting of shorter sequences)" );
         System.out.println( "   -" + GAP_RATIO_LENGTH_OPTION
-                + "=<decimal>  maximal allowed gap ratio per column (for deleting of columms) (0.0-1.0)" );
+                            + "=<decimal>  maximal allowed gap ratio per column (for deleting of columms) (0.0-1.0)" );
         System.out.println( "   -" + PERFORM_PHYLOGENETIC_INFERENCE
-                + "             to calculate a simple phylogenetic tree (Kimura distances, NJ)" );
+                            + "             to calculate a simple phylogenetic tree (Kimura distances, NJ)" );
+        System.out.println( "   -" + DO_NOT_NORMALIZE_FOR_EFF_LENGTH_OPTION
+                            + "            to normalize gap-contributions with MSA length, instead of individual effective sequence lenghts" );
+       
         System.out.println();
         System.out.println();
         System.out.println();