forester/java/src/org/forester/msa_compactor/MsaCompactor.java

   1 // $Id:
   2 // FORESTER -- software libraries and applications
   3 // for evolutionary biology research and applications.
   4 //
   5 // Copyright (C) 2014 Christian M. Zmasek
   6 // Copyright (C) 2014 Sanford-Burnham Medical Research Institute
   7 // All rights reserved
   8 //
   9 // This library is free software; you can redistribute it and/or
  10 // modify it under the terms of the GNU Lesser General Public
  11 // License as published by the Free Software Foundation; either
  12 // version 2.1 of the License, or (at your option) any later version.
  13 //
  14 // This library is distributed in the hope that it will be useful,
  15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  17 // Lesser General Public License for more details.
  18 //
  19 // You should have received a copy of the GNU Lesser General Public
  20 // License along with this library; if not, write to the Free Software
  21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  22 //
  23 // WWW: https://sites.google.com/site/cmzmasek/home/software/forester
  24
  25 package org.forester.msa_compactor;
  26
  27 import java.awt.Color;
  28 import java.io.File;
  29 import java.io.IOException;
  30 import java.io.Writer;
  31 import java.math.RoundingMode;
  32 import java.text.DecimalFormat;
  33 import java.text.NumberFormat;
  34 import java.util.ArrayList;
  35 import java.util.Arrays;
  36 import java.util.List;
  37 import java.util.SortedSet;
  38 import java.util.TreeSet;
  39
  40 import org.forester.archaeopteryx.Archaeopteryx;
  41 import org.forester.archaeopteryx.Configuration;
  42 import org.forester.evoinference.distance.NeighborJoiningF;
  43 import org.forester.evoinference.distance.PairwiseDistanceCalculator;
  44 import org.forester.evoinference.distance.PairwiseDistanceCalculator.PWD_DISTANCE_METHOD;
  45 import org.forester.evoinference.matrix.distance.BasicSymmetricalDistanceMatrix;
  46 import org.forester.evoinference.tools.BootstrapResampler;
  47 import org.forester.io.parsers.nhx.NHXParser.TAXONOMY_EXTRACTION;
  48 import org.forester.io.parsers.phyloxml.PhyloXmlDataFormatException;
  49 import org.forester.io.parsers.util.ParserUtils;
  50 import org.forester.io.writers.SequenceWriter;
  51 import org.forester.io.writers.SequenceWriter.SEQ_FORMAT;
  52 import org.forester.msa.DeleteableMsa;
  53 import org.forester.msa.Mafft;
  54 import org.forester.msa.Msa;
  55 import org.forester.msa.Msa.MSA_FORMAT;
  56 import org.forester.msa.MsaInferrer;
  57 import org.forester.msa.MsaMethods;
  58 import org.forester.msa.ResampleableMsa;
  59 import org.forester.phylogeny.Phylogeny;
  60 import org.forester.phylogeny.PhylogenyMethods;
  61 import org.forester.phylogeny.PhylogenyMethods.DESCENDANT_SORT_PRIORITY;
  62 import org.forester.phylogeny.PhylogenyNode;
  63 import org.forester.phylogeny.data.NodeVisualData;
  64 import org.forester.phylogeny.data.NodeVisualData.NodeFill;
  65 import org.forester.phylogeny.data.NodeVisualData.NodeShape;
  66 import org.forester.phylogeny.iterators.PhylogenyNodeIterator;
  67 import org.forester.sequence.Sequence;
  68 import org.forester.tools.ConfidenceAssessor;
  69 import org.forester.util.BasicDescriptiveStatistics;
  70 import org.forester.util.ForesterUtil;
  71
  72 public class MsaCompactor {
  73
  74     final private static NumberFormat NF_3                      = new DecimalFormat( "#.###" );
  75     final private static NumberFormat NF_4                      = new DecimalFormat( "#.####" );
  76     private double                    _gap_ratio                = -1;
  77     //
  78     private String                    _infile_name              = null;
  79     private final short               _longest_id_length;
  80     //
  81     private String                    _maffts_opts              = "--auto";
  82     private int                       _min_length               = -1;
  83     private DeleteableMsa             _msa                      = null;
  84     private boolean                   _norm                     = true;
  85     private File                      _out_file_base            = null;
  86     private MSA_FORMAT                _output_format            = MSA_FORMAT.FASTA;
  87     private String                    _path_to_mafft            = null;
  88     private boolean                   _phylogentic_inference    = false;
  89     //
  90     private boolean                   _realign                  = false;
  91     private final SortedSet<String>   _removed_seq_ids;
  92     private final ArrayList<Sequence> _removed_seqs;
  93     private File                      _removed_seqs_out_base    = null;
  94     private boolean                   _report_aln_mean_identity = false;
  95     private int                       _step                     = -1;
  96     private int                       _step_for_diagnostics     = -1;
  97     static {
  98         NF_4.setRoundingMode( RoundingMode.HALF_UP );
  99         NF_3.setRoundingMode( RoundingMode.HALF_UP );
 100     }
 101
 102     public MsaCompactor( final DeleteableMsa msa ) {
 103         _msa = msa;
 104         _removed_seq_ids = new TreeSet<String>();
 105         _longest_id_length = _msa.determineMaxIdLength();
 106         _removed_seqs = new ArrayList<Sequence>();
 107     }
 108
 109     public final Phylogeny calcTree() {
 110         final Phylogeny phy = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, _msa, false, "" );
 111         PhylogenyMethods.midpointRoot( phy );
 112         PhylogenyMethods.orderAppearance( phy.getRoot(), true, true, DESCENDANT_SORT_PRIORITY.NODE_NAME );
 113         final boolean x = PhylogenyMethods.extractFastaInformation( phy );
 114         if ( !x ) {
 115             final PhylogenyNodeIterator it = phy.iteratorExternalForward();
 116             while ( it.hasNext() ) {
 117                 final PhylogenyNode n = it.next();
 118                 final String name = n.getName().trim();
 119                 if ( !ForesterUtil.isEmpty( name ) ) {
 120                     try {
 121                         ParserUtils.extractTaxonomyDataFromNodeName( n, TAXONOMY_EXTRACTION.AGGRESSIVE );
 122                     }
 123                     catch ( final PhyloXmlDataFormatException e ) {
 124                         // Ignore.
 125                     }
 126                 }
 127             }
 128         }
 129         return phy;
 130     }
 131
 132     public final List<MsaProperties> chart( final int step, final boolean realign, final boolean norm )
 133             throws IOException, InterruptedException {
 134         final GapContribution stats[] = calcGapContribtionsStats( norm );
 135         final List<String> to_remove_ids = new ArrayList<String>();
 136         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
 137         for( final GapContribution gap_gontribution : stats ) {
 138             to_remove_ids.add( gap_gontribution.getId() );
 139         }
 140         Phylogeny phy = null;
 141         if ( _phylogentic_inference ) {
 142             System.out.println( "calculating phylogentic tree..." );
 143             System.out.println();
 144             phy = calcTree();
 145         }
 146         if ( !_realign ) {
 147             _step = -1;
 148         }
 149         int x = ForesterUtil.roundToInt( _msa.getNumberOfSequences() / 20.0 );
 150         if ( x < 1 ) {
 151             x = 1;
 152         }
 153         MsaProperties msa_prop = new MsaProperties( _msa, "", _report_aln_mean_identity );
 154         msa_props.add( msa_prop );
 155         printTableHeader();
 156         printMsaProperties( msa_prop );
 157         System.out.println();
 158         int i = 0;
 159         while ( _msa.getNumberOfSequences() > x ) {
 160             final String id = to_remove_ids.get( i );
 161             _msa.deleteRow( id, false );
 162             if ( realign && isPrintMsaStatsWriteOutfileAndRealign( i ) ) {
 163                 removeGapColumns();
 164                 realignWithMafft();
 165                 msa_prop = new MsaProperties( _msa, id, _report_aln_mean_identity );
 166                 msa_props.add( msa_prop );
 167                 printMsaProperties( msa_prop );
 168                 System.out.print( "(realigned)" );
 169                 System.out.println();
 170             }
 171             else if ( isPrintMsaStats( i ) ) {
 172                 removeGapColumns();
 173                 msa_prop = new MsaProperties( _msa, id, _report_aln_mean_identity );
 174                 msa_props.add( msa_prop );
 175                 printMsaProperties( msa_prop );
 176                 System.out.println();
 177             }
 178             ++i;
 179         }
 180         if ( _phylogentic_inference ) {
 181             decorateTree( phy, msa_props, true );
 182             displayTree( phy );
 183         }
 184         return msa_props;
 185     }
 186
 187     public final void decorateTree( final Phylogeny phy, final List<MsaProperties> msa_props, final boolean chart_only ) {
 188         final BasicDescriptiveStatistics length_stats = new BasicDescriptiveStatistics();
 189         for( int i = 0; i < msa_props.size(); ++i ) {
 190             final MsaProperties msa_prop = msa_props.get( i );
 191             final String id = msa_prop.getRemovedSeq();
 192             if ( !ForesterUtil.isEmpty( id ) ) {
 193                 length_stats.addValue( msa_prop.getLength() );
 194             }
 195         }
 196         final double mean = length_stats.arithmeticMean();
 197         final double min = length_stats.getMin();
 198         final double max = length_stats.getMax();
 199         final Color min_color = new Color( 0, 255, 0 );
 200         final Color max_color = new Color( 255, 0, 0 );
 201         final Color mean_color = new Color( 255, 255, 0 );
 202         final PhylogenyNodeIterator it = phy.iteratorExternalForward();
 203         if ( chart_only ) {
 204             while ( it.hasNext() ) {
 205                 final NodeVisualData vis = new NodeVisualData();
 206                 vis.setFillType( NodeFill.SOLID );
 207                 vis.setShape( NodeShape.RECTANGLE );
 208                 vis.setNodeColor( min_color );
 209                 it.next().getNodeData().setNodeVisualData( vis );
 210             }
 211         }
 212         for( int i = 0; i < msa_props.size(); ++i ) {
 213             final MsaProperties msa_prop = msa_props.get( i );
 214             final String id = msa_prop.getRemovedSeq();
 215             if ( !ForesterUtil.isEmpty( id ) ) {
 216                 final PhylogenyNode n = phy.getNode( id );
 217                 n.setName( n.getName() + " [" + i + "]" );
 218                 if ( !chart_only ) {
 219                     final NodeVisualData vis = new NodeVisualData();
 220                     vis.setFillType( NodeFill.SOLID );
 221                     vis.setShape( NodeShape.RECTANGLE );
 222                     vis.setNodeColor( ForesterUtil.calcColor( msa_prop.getLength(), min, max, mean_color, max_color ) );
 223                     n.getNodeData().setNodeVisualData( vis );
 224                 }
 225                 else {
 226                     n.getNodeData()
 227                             .getNodeVisualData()
 228                             .setNodeColor( ForesterUtil.calcColor( msa_prop.getLength(),
 229                                                                    min,
 230                                                                    max,
 231                                                                    mean,
 232                                                                    min_color,
 233                                                                    max_color,
 234                                                                    mean_color ) );
 235                 }
 236             }
 237         }
 238     }
 239
 240     final public void deleteGapColumns( final double max_allowed_gap_ratio ) {
 241         _msa.deleteGapColumns( max_allowed_gap_ratio );
 242     }
 243
 244     public final void displayTree( final Phylogeny phy ) {
 245         final Configuration config = new Configuration();
 246         config.setDisplayAsPhylogram( true );
 247         config.setUseStyle( true );
 248         config.setDisplayTaxonomyCode( false );
 249         config.setDisplayTaxonomyCommonNames( false );
 250         config.setDisplayTaxonomyScientificNames( false );
 251         config.setDisplaySequenceNames( false );
 252         config.setDisplaySequenceSymbols( false );
 253         config.setDisplayGeneNames( false );
 254         config.setShowScale( true );
 255         config.setAddTaxonomyImagesCB( false );
 256         config.setBaseFontSize( 9 );
 257         config.setBaseFontFamilyName( "Arial" );
 258         Archaeopteryx.createApplication( phy, config, _infile_name );
 259     }
 260
 261     final public Msa getMsa() {
 262         return _msa;
 263     }
 264
 265     public final void removeSequencesByMinimalLength( final int min_effective_length ) {
 266         printMsaProperties( new MsaProperties( _msa, "", _report_aln_mean_identity ) );
 267         System.out.println();
 268         _msa = DeleteableMsa.createInstance( MsaMethods.removeSequencesByMinimalLength( _msa, min_effective_length ) );
 269         removeGapColumns();
 270         printMsaProperties( new MsaProperties( _msa, "", _report_aln_mean_identity ) );
 271         System.out.println();
 272     }
 273
 274     public final List<MsaProperties> removeViaGapAverage( final double mean_gapiness ) throws IOException,
 275             InterruptedException {
 276         final GapContribution stats[] = calcGapContribtionsStats( _norm );
 277         final List<String> to_remove_ids = new ArrayList<String>();
 278         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
 279         for( final GapContribution gap_gontribution : stats ) {
 280             to_remove_ids.add( gap_gontribution.getId() );
 281         }
 282         Phylogeny phy = null;
 283         if ( _phylogentic_inference ) {
 284             System.out.println( "calculating phylogentic tree..." );
 285             System.out.println();
 286             phy = calcTree();
 287         }
 288         printTableHeader();
 289         MsaProperties msa_prop = new MsaProperties( _msa, "", _report_aln_mean_identity );
 290         msa_props.add( msa_prop );
 291         printMsaProperties( msa_prop );
 292         System.out.println();
 293         int i = 0;
 294         while ( MsaMethods.calcGapRatio( _msa ) > mean_gapiness ) {
 295             final String id = to_remove_ids.get( i );
 296             _removed_seq_ids.add( id );
 297             final Sequence deleted = _msa.deleteRow( id, true );
 298             _removed_seqs.add( deleted );
 299             removeGapColumns();
 300             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( MsaMethods.calcGapRatio( _msa ) <= mean_gapiness ) ) {
 301                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
 302                 msa_props.add( msa_prop );
 303                 System.out.println();
 304             }
 305             else if ( isPrintMsaStats( i ) ) {
 306                 msa_prop = new MsaProperties( _msa, id, _report_aln_mean_identity );
 307                 msa_props.add( msa_prop );
 308                 printMsaProperties( msa_prop );
 309                 System.out.println();
 310             }
 311             ++i;
 312         }
 313         if ( _removed_seqs_out_base != null ) {
 314             final String msg = writeAndAlignRemovedSeqs();
 315             System.out.println();
 316             System.out.println( msg );
 317         }
 318         if ( _phylogentic_inference ) {
 319             decorateTree( phy, msa_props, false );
 320             displayTree( phy );
 321         }
 322         return msa_props;
 323     }
 324
 325     public List<MsaProperties> removeViaLength( final int length ) throws IOException, InterruptedException {
 326         final GapContribution stats[] = calcGapContribtionsStats( _norm );
 327         final List<String> to_remove_ids = new ArrayList<String>();
 328         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
 329         for( final GapContribution gap_gontribution : stats ) {
 330             to_remove_ids.add( gap_gontribution.getId() );
 331         }
 332         Phylogeny phy = null;
 333         if ( _phylogentic_inference ) {
 334             System.out.println( "calculating phylogentic tree..." );
 335             System.out.println();
 336             phy = calcTree();
 337         }
 338         printTableHeader();
 339         MsaProperties msa_prop = new MsaProperties( _msa, "", _report_aln_mean_identity );
 340         msa_props.add( msa_prop );
 341         printMsaProperties( msa_prop );
 342         System.out.println();
 343         int i = 0;
 344         while ( _msa.getLength() > length ) {
 345             final String id = to_remove_ids.get( i );
 346             _removed_seq_ids.add( id );
 347             final Sequence deleted = _msa.deleteRow( id, true );
 348             _removed_seqs.add( deleted );
 349             removeGapColumns();
 350             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( _msa.getLength() <= length ) ) {
 351                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
 352                 msa_props.add( msa_prop );
 353                 System.out.println();
 354             }
 355             else if ( isPrintMsaStats( i ) ) {
 356                 msa_prop = new MsaProperties( _msa, id, _report_aln_mean_identity );
 357                 printMsaProperties( msa_prop );
 358                 msa_props.add( msa_prop );
 359                 System.out.println();
 360             }
 361             ++i;
 362         }
 363         if ( _removed_seqs_out_base != null ) {
 364             final String msg = writeAndAlignRemovedSeqs();
 365             System.out.println();
 366             System.out.println( msg );
 367         }
 368         if ( _phylogentic_inference ) {
 369             decorateTree( phy, msa_props, false );
 370             displayTree( phy );
 371         }
 372         return msa_props;
 373     }
 374
 375     public final List<MsaProperties> removeWorstOffenders( final int to_remove ) throws IOException,
 376             InterruptedException {
 377         final GapContribution stats[] = calcGapContribtionsStats( _norm );
 378         final List<String> to_remove_ids = new ArrayList<String>();
 379         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
 380         for( int j = 0; j < to_remove; ++j ) {
 381             to_remove_ids.add( stats[ j ].getId() );
 382         }
 383         Phylogeny phy = null;
 384         if ( _phylogentic_inference ) {
 385             System.out.println( "calculating phylogentic tree..." );
 386             System.out.println();
 387             phy = calcTree();
 388         }
 389         printTableHeader();
 390         MsaProperties msa_prop = new MsaProperties( _msa, "", _report_aln_mean_identity );
 391         msa_props.add( msa_prop );
 392         printMsaProperties( msa_prop );
 393         System.out.println();
 394         for( int i = 0; i < to_remove_ids.size(); ++i ) {
 395             final String id = to_remove_ids.get( i );
 396             _removed_seq_ids.add( id );
 397             final Sequence deleted = _msa.deleteRow( id, true );
 398             _removed_seqs.add( deleted );
 399             removeGapColumns();
 400             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( i == ( to_remove_ids.size() - 1 ) ) ) {
 401                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
 402                 msa_props.add( msa_prop );
 403                 System.out.println();
 404             }
 405             else if ( isPrintMsaStats( i ) ) {
 406                 msa_prop = new MsaProperties( _msa, id, _report_aln_mean_identity );
 407                 msa_props.add( msa_prop );
 408                 printMsaProperties( msa_prop );
 409                 System.out.println();
 410             }
 411         }
 412         if ( _removed_seqs_out_base != null ) {
 413             final String msg = writeAndAlignRemovedSeqs();
 414             System.out.println();
 415             System.out.println( msg );
 416         }
 417         if ( _phylogentic_inference ) {
 418             decorateTree( phy, msa_props, false );
 419             displayTree( phy );
 420         }
 421         return msa_props;
 422     }
 423
 424     public final void setGapRatio( final double gap_ratio ) {
 425         _gap_ratio = gap_ratio;
 426     }
 427
 428     public void setInfileName( final String infile_name ) {
 429         _infile_name = infile_name;
 430     }
 431
 432     public final void setMafftOptions( final String maffts_opts ) {
 433         _maffts_opts = maffts_opts;
 434     }
 435
 436     public final void setMinLength( final int min_length ) {
 437         _min_length = min_length;
 438     }
 439
 440     public final void setNorm( final boolean norm ) {
 441         _norm = norm;
 442     }
 443
 444     final public void setOutFileBase( final File out_file_base ) {
 445         _out_file_base = out_file_base;
 446     }
 447
 448     public final void setOutputFormat( final MSA_FORMAT output_format ) {
 449         _output_format = output_format;
 450     }
 451
 452     public void setPathToMafft( final String path_to_mafft ) {
 453         _path_to_mafft = path_to_mafft;
 454     }
 455
 456     public void setPeformPhylogenticInference( final boolean phylogentic_inference ) {
 457         _phylogentic_inference = phylogentic_inference;
 458     }
 459
 460     public final void setRealign( final boolean realign ) {
 461         _realign = realign;
 462     }
 463
 464     public final void setRemovedSeqsOutBase( final File removed_seqs_out_base ) {
 465         _removed_seqs_out_base = removed_seqs_out_base;
 466     }
 467
 468     public final void setReportAlnMeanIdentity( final boolean report_aln_mean_identity ) {
 469         _report_aln_mean_identity = report_aln_mean_identity;
 470     }
 471
 472     public final void setStep( final int step ) {
 473         _step = step;
 474     }
 475
 476     public final void setStepForDiagnostics( final int step_for_diagnostics ) {
 477         _step_for_diagnostics = step_for_diagnostics;
 478     }
 479
 480     final public String writeAndAlignRemovedSeqs() throws IOException, InterruptedException {
 481         final StringBuilder msg = new StringBuilder();
 482         final String n = _removed_seqs_out_base + "_" + _removed_seqs.size() + ".fasta";
 483         SequenceWriter.writeSeqs( _removed_seqs, new File( n ), SEQ_FORMAT.FASTA, 100 );
 484         msg.append( "wrote " + _removed_seqs.size() + " removed sequences to " + "\"" + n + "\"" );
 485         if ( _realign ) {
 486             final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
 487             final List<String> opts = new ArrayList<String>();
 488             for( final String o : _maffts_opts.split( "\\s" ) ) {
 489                 opts.add( o );
 490             }
 491             final Msa removed_msa = mafft.infer( _removed_seqs, opts );
 492             final Double gr = MsaMethods.calcGapRatio( removed_msa );
 493             String s = _removed_seqs_out_base + "_" + removed_msa.getNumberOfSequences() + "_"
 494                     + removed_msa.getLength() + "_" + ForesterUtil.roundToInt( gr * 100 );
 495             final String suffix = obtainSuffix();
 496             s += suffix;
 497             writeMsa( removed_msa, s, _output_format );
 498             msg.append( ", and as MSA of length " + removed_msa.getLength() + " to \"" + s + "\"" );
 499         }
 500         return msg.toString();
 501     }
 502
 503     final public String writeMsa( final File outfile ) throws IOException {
 504         final Double gr = MsaMethods.calcGapRatio( _msa );
 505         final String s = outfile + "_" + _msa.getNumberOfSequences() + "_" + _msa.getLength() + "_"
 506                 + ForesterUtil.roundToInt( gr * 100 );
 507         writeMsa( _msa, s + obtainSuffix(), _output_format );
 508         return s;
 509     }
 510
 511     final int calcNonGapResidues( final Sequence seq ) {
 512         int ng = 0;
 513         for( int i = 0; i < seq.getLength(); ++i ) {
 514             if ( !seq.isGapAt( i ) ) {
 515                 ++ng;
 516             }
 517         }
 518         return ng;
 519     }
 520
 521     private final GapContribution[] calcGapContribtions( final boolean normalize_for_effective_seq_length ) {
 522         final double gappiness[] = calcGappiness();
 523         final GapContribution stats[] = new GapContribution[ _msa.getNumberOfSequences() ];
 524         for( int row = 0; row < _msa.getNumberOfSequences(); ++row ) {
 525             stats[ row ] = new GapContribution( _msa.getIdentifier( row ) );
 526             for( int col = 0; col < _msa.getLength(); ++col ) {
 527                 if ( !_msa.isGapAt( row, col ) ) {
 528                     stats[ row ].addToValue( gappiness[ col ] );
 529                 }
 530             }
 531             if ( normalize_for_effective_seq_length ) {
 532                 stats[ row ].divideValue( calcNonGapResidues( _msa.getSequence( row ) ) );
 533             }
 534             else {
 535                 stats[ row ].divideValue( _msa.getLength() );
 536             }
 537         }
 538         return stats;
 539     }
 540
 541     final private GapContribution[] calcGapContribtionsStats( final boolean norm ) {
 542         final GapContribution stats[] = calcGapContribtions( norm );
 543         Arrays.sort( stats );
 544         return stats;
 545     }
 546
 547     private final double[] calcGappiness() {
 548         final int l = _msa.getLength();
 549         final double gappiness[] = new double[ l ];
 550         final int seqs = _msa.getNumberOfSequences();
 551         for( int i = 0; i < l; ++i ) {
 552             gappiness[ i ] = ( double ) MsaMethods.calcGapSumPerColumn( _msa, i ) / seqs;
 553         }
 554         return gappiness;
 555     }
 556
 557     private final Phylogeny inferNJphylogeny( final PWD_DISTANCE_METHOD pwd_distance_method,
 558                                               final Msa msa,
 559                                               final boolean write_matrix,
 560                                               final String matrix_name ) {
 561         BasicSymmetricalDistanceMatrix m = null;
 562         switch ( pwd_distance_method ) {
 563             case KIMURA_DISTANCE:
 564                 m = PairwiseDistanceCalculator.calcKimuraDistances( msa );
 565                 break;
 566             case POISSON_DISTANCE:
 567                 m = PairwiseDistanceCalculator.calcPoissonDistances( msa );
 568                 break;
 569             case FRACTIONAL_DISSIMILARITY:
 570                 m = PairwiseDistanceCalculator.calcFractionalDissimilarities( msa );
 571                 break;
 572             default:
 573                 throw new IllegalArgumentException( "invalid pwd method" );
 574         }
 575         if ( write_matrix ) {
 576             try {
 577                 m.write( ForesterUtil.createBufferedWriter( matrix_name ) );
 578             }
 579             catch ( final IOException e ) {
 580                 e.printStackTrace();
 581             }
 582         }
 583         final NeighborJoiningF nj = NeighborJoiningF.createInstance( false, 5 );
 584         final Phylogeny phy = nj.execute( m );
 585         return phy;
 586     }
 587
 588     private final boolean isPrintMsaStats( final int i ) {
 589         return ( ( ( _step == 1 ) && ( _step_for_diagnostics == 1 ) ) || ( ( _step_for_diagnostics > 0 ) && ( ( ( i + 1 ) % _step_for_diagnostics ) == 0 ) ) );
 590     }
 591
 592     private final boolean isPrintMsaStatsWriteOutfileAndRealign( final int i ) {
 593         return ( ( ( _step == 1 ) && ( _step_for_diagnostics == 1 ) ) || ( ( _step > 0 ) && ( ( ( i + 1 ) % _step ) == 0 ) ) );
 594     }
 595
 596     private final StringBuilder msaPropertiesAsSB( final MsaProperties msa_properties ) {
 597         final StringBuilder sb = new StringBuilder();
 598         sb.append( msa_properties.getNumberOfSequences() );
 599         sb.append( "\t" );
 600         sb.append( msa_properties.getLength() );
 601         sb.append( "\t" );
 602         sb.append( NF_4.format( msa_properties.getGapRatio() ) );
 603         if ( _report_aln_mean_identity ) {
 604             sb.append( "\t" );
 605             sb.append( NF_4.format( msa_properties.getAverageIdentityRatio() ) );
 606         }
 607         return sb;
 608     }
 609
 610     private String obtainSuffix() {
 611         if ( _output_format == MSA_FORMAT.FASTA ) {
 612             return ".fasta";
 613         }
 614         else if ( _output_format == MSA_FORMAT.PHYLIP ) {
 615             return ".aln";
 616         }
 617         return "";
 618     }
 619
 620     private final Phylogeny pi( final String matrix, final int boostrap ) {
 621         final Phylogeny master_phy = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, _msa, true, matrix );
 622         final int seed = 15;
 623         final int n = 100;
 624         final ResampleableMsa resampleable_msa = new ResampleableMsa( _msa );
 625         final int[][] resampled_column_positions = BootstrapResampler.createResampledColumnPositions( _msa.getLength(),
 626                                                                                                       n,
 627                                                                                                       seed );
 628         final Phylogeny[] eval_phys = new Phylogeny[ n ];
 629         for( int i = 0; i < n; ++i ) {
 630             resampleable_msa.resample( resampled_column_positions[ i ] );
 631             eval_phys[ i ] = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, resampleable_msa, false, null );
 632         }
 633         ConfidenceAssessor.evaluate( "bootstrap", eval_phys, master_phy, true, 1 );
 634         PhylogenyMethods.extractFastaInformation( master_phy );
 635         return master_phy;
 636     }
 637
 638     private final void printMsaProperties( final MsaProperties msa_properties ) {
 639         if ( ( _step == 1 ) || ( _step_for_diagnostics == 1 ) ) {
 640             System.out.print( ForesterUtil.pad( msa_properties.getRemovedSeq(), _longest_id_length, ' ', false ) );
 641             System.out.print( "\t" );
 642         }
 643         System.out.print( msaPropertiesAsSB( msa_properties ) );
 644         System.out.print( "\t" );
 645     }
 646
 647     final private MsaProperties printMsaStatsWriteOutfileAndRealign( final boolean realign, final String id )
 648             throws IOException, InterruptedException {
 649         if ( realign ) {
 650             realignWithMafft();
 651         }
 652         final MsaProperties msa_prop = new MsaProperties( _msa, id, _report_aln_mean_identity );
 653         printMsaProperties( msa_prop );
 654         final String s = writeOutfile();
 655         System.out.print( "-> " + s + ( realign ? "\t(realigned)" : "" ) );
 656         return msa_prop;
 657     }
 658
 659     private final void printTableHeader() {
 660         if ( ( _step == 1 ) || ( _step_for_diagnostics == 1 ) ) {
 661             System.out.print( ForesterUtil.pad( "Id", _longest_id_length, ' ', false ) );
 662             System.out.print( "\t" );
 663         }
 664         System.out.print( "Seqs" );
 665         System.out.print( "\t" );
 666         System.out.print( "Length" );
 667         System.out.print( "\t" );
 668         System.out.print( "Gaps" );
 669         System.out.print( "\t" );
 670         if ( _report_aln_mean_identity ) {
 671             System.out.print( "MSA qual" );
 672             System.out.print( "\t" );
 673         }
 674         System.out.println();
 675     }
 676
 677     final private void realignWithMafft() throws IOException, InterruptedException {
 678         final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
 679         final List<String> opts = new ArrayList<String>();
 680         for( final String o : _maffts_opts.split( "\\s" ) ) {
 681             opts.add( o );
 682         }
 683         _msa = DeleteableMsa.createInstance( mafft.infer( _msa.asSequenceList(), opts ) );
 684     }
 685
 686     final private void removeGapColumns() {
 687         _msa.deleteGapOnlyColumns();
 688     }
 689
 690     private final String writeOutfile() throws IOException {
 691         final String s = writeMsa( _out_file_base );
 692         return s;
 693     }
 694
 695     // Returns null if not path found.
 696     final public static String guessPathToMafft() {
 697         String path;
 698         if ( ForesterUtil.OS_NAME.toLowerCase().indexOf( "win" ) >= 0 ) {
 699             path = "C:\\Program Files\\mafft-win\\mafft.bat";
 700             if ( MsaInferrer.isInstalled( path ) ) {
 701                 return path;
 702             }
 703         }
 704         path = "/home/czmasek/SOFTWARE/MSA/MAFFT/mafft-7.130-without-extensions/scripts/mafft";
 705         if ( MsaInferrer.isInstalled( path ) ) {
 706             return path;
 707         }
 708         path = "/usr/local/bin/mafft";
 709         if ( MsaInferrer.isInstalled( path ) ) {
 710             return path;
 711         }
 712         path = "/usr/bin/mafft";
 713         if ( MsaInferrer.isInstalled( path ) ) {
 714             return path;
 715         }
 716         path = "/bin/mafft";
 717         if ( MsaInferrer.isInstalled( path ) ) {
 718             return path;
 719         }
 720         path = "mafft";
 721         if ( MsaInferrer.isInstalled( path ) ) {
 722             return path;
 723         }
 724         return null;
 725     }
 726
 727     final private static void writeMsa( final Msa msa, final String outfile, final MSA_FORMAT format )
 728             throws IOException {
 729         final Writer w = ForesterUtil.createBufferedWriter( outfile );
 730         msa.write( w, format );
 731         w.close();
 732     }
 733 }