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