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