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.io.File;
  28 import java.io.IOException;
  29 import java.io.Writer;
  30 import java.math.RoundingMode;
  31 import java.text.DecimalFormat;
  32 import java.text.NumberFormat;
  33 import java.util.ArrayList;
  34 import java.util.Arrays;
  35 import java.util.List;
  36 import java.util.SortedSet;
  37 import java.util.TreeSet;
  38
  39 import org.forester.evoinference.distance.NeighborJoiningF;
  40 import org.forester.evoinference.distance.PairwiseDistanceCalculator;
  41 import org.forester.evoinference.distance.PairwiseDistanceCalculator.PWD_DISTANCE_METHOD;
  42 import org.forester.evoinference.matrix.distance.BasicSymmetricalDistanceMatrix;
  43 import org.forester.evoinference.tools.BootstrapResampler;
  44 import org.forester.io.writers.SequenceWriter;
  45 import org.forester.io.writers.SequenceWriter.SEQ_FORMAT;
  46 import org.forester.msa.DeleteableMsa;
  47 import org.forester.msa.Mafft;
  48 import org.forester.msa.Msa;
  49 import org.forester.msa.Msa.MSA_FORMAT;
  50 import org.forester.msa.MsaInferrer;
  51 import org.forester.msa.MsaMethods;
  52 import org.forester.msa.ResampleableMsa;
  53 import org.forester.phylogeny.Phylogeny;
  54 import org.forester.phylogeny.PhylogenyMethods;
  55 import org.forester.sequence.Sequence;
  56 import org.forester.tools.ConfidenceAssessor;
  57 import org.forester.util.ForesterUtil;
  58
  59 public class MsaCompactor {
  60
  61     final private static NumberFormat NF_3                   = new DecimalFormat( "#.###" );
  62     final private static NumberFormat NF_4                   = new DecimalFormat( "#.####" );
  63     private double                    _gap_ratio             = -1;
  64     //
  65     private final String              _maffts_opts           = "--auto";
  66     private int                       _min_length            = -1;
  67     //
  68     private DeleteableMsa             _msa                   = null;
  69     private boolean                   _norm                  = true;
  70     private File                      _out_file_base         = null;
  71     private MSA_FORMAT                _output_format         = MSA_FORMAT.FASTA;
  72     private String                    _path_to_mafft         = null;
  73     //
  74     private boolean                   _realign               = false;
  75     private final SortedSet<String>   _removed_seq_ids;
  76     private File                      _removed_seqs_out_base = null;
  77
  78     public final void setRemovedSeqsOutBase( final File removed_seqs_out_base ) {
  79         _removed_seqs_out_base = removed_seqs_out_base;
  80     }
  81     private boolean                   _report_aln_mean_identity = false;
  82     private int                       _step                     = -1;
  83     private int                       _step_for_diagnostics     = -1;
  84     private final short               _longest_id_length;
  85     private final ArrayList<Sequence> _removed_seqs;
  86     static {
  87         NF_4.setRoundingMode( RoundingMode.HALF_UP );
  88         NF_3.setRoundingMode( RoundingMode.HALF_UP );
  89     }
  90
  91     public MsaCompactor( final DeleteableMsa msa ) {
  92         _msa = msa;
  93         _removed_seq_ids = new TreeSet<String>();
  94         _longest_id_length = _msa.determineMaxIdLength();
  95         _removed_seqs = new ArrayList<Sequence>();
  96     }
  97
  98     final public Msa getMsa() {
  99         return _msa;
 100     }
 101
 102     final public SortedSet<String> getRemovedSeqIds() {
 103         return _removed_seq_ids;
 104     }
 105
 106     public final void removeViaGapAverage( final double mean_gapiness ) throws IOException, InterruptedException {
 107         final GapContribution stats[] = calcGapContribtionsStats( _norm );
 108         final List<String> to_remove_ids = new ArrayList<String>();
 109         for( final GapContribution gap_gontribution : stats ) {
 110             to_remove_ids.add( gap_gontribution.getId() );
 111         }
 112         printTableHeader();
 113         int i = 0;
 114         while ( MsaMethods.calcGapRatio( _msa ) > mean_gapiness ) {
 115             final String id = to_remove_ids.get( i );
 116             _removed_seq_ids.add( id );
 117             final Sequence deleted = _msa.deleteRow( id, true );
 118             _removed_seqs.add( deleted );
 119             removeGapColumns();
 120             if ( ( ( _step > 0 ) && ( ( ( i + 1 ) % _step ) == 0 ) )
 121                     || ( MsaMethods.calcGapRatio( _msa ) <= mean_gapiness ) ) {
 122                 printMsaStatsWriteOutfileAndRealign( _realign, id );
 123             }
 124             else {
 125                 final MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
 126                 printMsaProperties( id, msa_prop );
 127             }
 128             System.out.println();
 129             ++i;
 130         }
 131         if ( _removed_seqs_out_base != null ) {
 132             final String msg = writeAndAlignRemovedSeqs();
 133             System.out.println( msg );
 134         }
 135     }
 136
 137     public void removeViaLength( final int length ) throws IOException, InterruptedException {
 138         final GapContribution stats[] = calcGapContribtionsStats( _norm );
 139         final List<String> to_remove_ids = new ArrayList<String>();
 140         for( final GapContribution gap_gontribution : stats ) {
 141             to_remove_ids.add( gap_gontribution.getId() );
 142         }
 143         printTableHeader();
 144         int i = 0;
 145         while ( _msa.getLength() > length ) {
 146             final String id = to_remove_ids.get( i );
 147             _removed_seq_ids.add( id );
 148             final Sequence deleted = _msa.deleteRow( id, true );
 149             _removed_seqs.add( deleted );
 150             removeGapColumns();
 151             if ( ( ( _step > 0 ) && ( ( ( i + 1 ) % _step ) == 0 ) ) || ( _msa.getLength() <= length ) ) {
 152                 printMsaStatsWriteOutfileAndRealign( _realign, id );
 153             }
 154             final MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
 155             printMsaProperties( id, msa_prop );
 156             System.out.println();
 157             ++i;
 158         }
 159         if ( _removed_seqs_out_base != null ) {
 160             final String msg = writeAndAlignRemovedSeqs();
 161             System.out.println( msg );
 162         }
 163     }
 164
 165     public final void removeWorstOffenders( final int to_remove ) throws IOException, InterruptedException {
 166         final GapContribution stats[] = calcGapContribtionsStats( _norm );
 167         final List<String> to_remove_ids = new ArrayList<String>();
 168         for( int j = 0; j < to_remove; ++j ) {
 169             to_remove_ids.add( stats[ j ].getId() );
 170             _removed_seq_ids.add( stats[ j ].getId() );
 171         }
 172         printTableHeader();
 173         for( int i = 0; i < to_remove_ids.size(); ++i ) {
 174             final String id = to_remove_ids.get( i );
 175             _removed_seq_ids.add( id );
 176             final Sequence deleted = _msa.deleteRow( id, true );
 177             _removed_seqs.add( deleted );
 178             removeGapColumns();
 179             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( i == ( to_remove_ids.size() - 1 ) ) ) {
 180                 printMsaStatsWriteOutfileAndRealign( _realign, id );
 181                 System.out.println();
 182             }
 183             else if ( isPrintMsaStats( i ) ) {
 184                 final MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
 185                 printMsaProperties( id, msa_prop );
 186                 System.out.println();
 187             }
 188         }
 189         if ( _removed_seqs_out_base != null ) {
 190             final String msg = writeAndAlignRemovedSeqs();
 191             System.out.println( msg );
 192         }
 193     }
 194
 195     public final List<MsaProperties> chart( final int step, final boolean realign, final boolean norm )
 196             throws IOException, InterruptedException {
 197         final GapContribution stats[] = calcGapContribtionsStats( norm );
 198         final List<String> to_remove_ids = new ArrayList<String>();
 199         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
 200         for( final GapContribution gap_gontribution : stats ) {
 201             to_remove_ids.add( gap_gontribution.getId() );
 202         }
 203         printTableHeader();
 204         int i = 0;
 205         final int x = ForesterUtil.roundToInt( _msa.getNumberOfSequences() / 20.0 );
 206         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
 207         msa_props.add( msa_prop );
 208         printMsaProperties( "", msa_prop );
 209         System.out.println();
 210         while ( _msa.getNumberOfSequences() > x ) {
 211             final String id = to_remove_ids.get( i );
 212             _msa.deleteRow( id, false );
 213             if ( realign && isPrintMsaStatsWriteOutfileAndRealign( i ) ) {
 214                 removeGapColumns();
 215                 realignWithMafft();
 216                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
 217                 msa_props.add( msa_prop );
 218                 printMsaProperties( id, msa_prop );
 219                 System.out.print( "(realigned)" );
 220                 System.out.println();
 221             }
 222             else if ( isPrintMsaStats( i ) ) {
 223                 removeGapColumns();
 224                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
 225                 msa_props.add( msa_prop );
 226                 printMsaProperties( id, msa_prop );
 227                 System.out.println();
 228             }
 229             ++i;
 230         }
 231         return msa_props;
 232     }
 233
 234     private final boolean isPrintMsaStats( final int i ) {
 235         return ( ( _step_for_diagnostics < 2 ) || ( ( ( i + 1 ) % _step_for_diagnostics ) == 0 ) );
 236     }
 237
 238     private final boolean isPrintMsaStatsWriteOutfileAndRealign( final int i ) {
 239         return ( ( _step < 2 ) || ( ( ( i + 1 ) % _step ) == 0 ) );
 240     }
 241
 242     public final void setGapRatio( final double gap_ratio ) {
 243         _gap_ratio = gap_ratio;
 244     }
 245
 246     public final void setMinLength( final int min_length ) {
 247         _min_length = min_length;
 248     }
 249
 250     public final void setNorm( final boolean norm ) {
 251         _norm = norm;
 252     }
 253
 254     final public void setOutFileBase( final File out_file_base ) {
 255         _out_file_base = out_file_base;
 256     }
 257
 258     public final void setOutputFormat( final MSA_FORMAT output_format ) {
 259         _output_format = output_format;
 260     }
 261
 262     public void setPathToMafft( final String path_to_mafft ) {
 263         _path_to_mafft = path_to_mafft;
 264     }
 265
 266     public final void setRealign( final boolean realign ) {
 267         _realign = realign;
 268     }
 269
 270     public final void setStep( final int step ) {
 271         _step = step;
 272     }
 273
 274     public final void setStepForDiagnostics( final int step_for_diagnostics ) {
 275         _step_for_diagnostics = step_for_diagnostics;
 276     }
 277
 278     public final void setReportAlnMeanIdentity( final boolean report_aln_mean_identity ) {
 279         _report_aln_mean_identity = report_aln_mean_identity;
 280     }
 281
 282     final public String writeMsa( final File outfile, final MSA_FORMAT format, final String suffix ) throws IOException {
 283         final Double gr = MsaMethods.calcGapRatio( _msa );
 284         final String s = outfile + "_" + _msa.getNumberOfSequences() + "_" + _msa.getLength() + "_"
 285                 + ForesterUtil.roundToInt( gr * 100 );
 286         writeMsa( _msa, s + suffix, format );
 287         return s;
 288     }
 289
 290     final public String writeAndAlignRemovedSeqs() throws IOException, InterruptedException {
 291         final StringBuilder msg = new StringBuilder();
 292         final String n = _removed_seqs_out_base + "_" + _removed_seqs.size() + ".fasta";
 293         SequenceWriter.writeSeqs( _removed_seqs, new File( n ), SEQ_FORMAT.FASTA, 100 );
 294         msg.append( "wrote " + _removed_seqs.size() + " removed sequences to " + " to \"" + n + "\"" );
 295         if ( _realign ) {
 296             final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
 297             final List<String> opts = new ArrayList<String>();
 298             for( final String o : _maffts_opts.split( "\\s" ) ) {
 299                 opts.add( o );
 300             }
 301             final Msa removed_msa = mafft.infer( _removed_seqs, opts );
 302             final Double gr = MsaMethods.calcGapRatio( removed_msa );
 303             String s = _removed_seqs_out_base + "_" + removed_msa.getNumberOfSequences() + "_"
 304                     + removed_msa.getLength() + "_" + ForesterUtil.roundToInt( gr * 100 );
 305             String suffix = "";
 306             if ( _output_format == MSA_FORMAT.FASTA ) {
 307                 suffix = ".fasta";
 308             }
 309             else if ( _output_format == MSA_FORMAT.PHYLIP ) {
 310                 suffix = ".aln";
 311             }
 312             s += suffix;
 313             writeMsa( removed_msa, s, _output_format );
 314             msg.append( ", and as MSA of length " + removed_msa.getLength() + " to \"" + s + "\"" );
 315         }
 316         return msg.toString();
 317     }
 318
 319     final int calcNonGapResidues( final Sequence seq ) {
 320         int ng = 0;
 321         for( int i = 0; i < seq.getLength(); ++i ) {
 322             if ( !seq.isGapAt( i ) ) {
 323                 ++ng;
 324             }
 325         }
 326         return ng;
 327     }
 328
 329     private final GapContribution[] calcGapContribtions( final boolean normalize_for_effective_seq_length ) {
 330         final double gappiness[] = calcGappiness();
 331         final GapContribution stats[] = new GapContribution[ _msa.getNumberOfSequences() ];
 332         for( int row = 0; row < _msa.getNumberOfSequences(); ++row ) {
 333             stats[ row ] = new GapContribution( _msa.getIdentifier( row ) );
 334             for( int col = 0; col < _msa.getLength(); ++col ) {
 335                 if ( !_msa.isGapAt( row, col ) ) {
 336                     stats[ row ].addToValue( gappiness[ col ] );
 337                 }
 338             }
 339             if ( normalize_for_effective_seq_length ) {
 340                 stats[ row ].divideValue( calcNonGapResidues( _msa.getSequence( row ) ) );
 341             }
 342             else {
 343                 stats[ row ].divideValue( _msa.getLength() );
 344             }
 345         }
 346         return stats;
 347     }
 348
 349     final private GapContribution[] calcGapContribtionsStats( final boolean norm ) {
 350         final GapContribution stats[] = calcGapContribtions( norm );
 351         Arrays.sort( stats );
 352         return stats;
 353     }
 354
 355     private final double[] calcGappiness() {
 356         final int l = _msa.getLength();
 357         final double gappiness[] = new double[ l ];
 358         final int seqs = _msa.getNumberOfSequences();
 359         for( int i = 0; i < l; ++i ) {
 360             gappiness[ i ] = ( double ) MsaMethods.calcGapSumPerColumn( _msa, i ) / seqs;
 361         }
 362         return gappiness;
 363     }
 364
 365     private final Phylogeny inferNJphylogeny( final PWD_DISTANCE_METHOD pwd_distance_method,
 366                                               final Msa msa,
 367                                               final boolean write_matrix,
 368                                               final String matrix_name ) {
 369         BasicSymmetricalDistanceMatrix m = null;
 370         switch ( pwd_distance_method ) {
 371             case KIMURA_DISTANCE:
 372                 m = PairwiseDistanceCalculator.calcKimuraDistances( msa );
 373                 break;
 374             case POISSON_DISTANCE:
 375                 m = PairwiseDistanceCalculator.calcPoissonDistances( msa );
 376                 break;
 377             case FRACTIONAL_DISSIMILARITY:
 378                 m = PairwiseDistanceCalculator.calcFractionalDissimilarities( msa );
 379                 break;
 380             default:
 381                 throw new IllegalArgumentException( "invalid pwd method" );
 382         }
 383         if ( write_matrix ) {
 384             try {
 385                 m.write( ForesterUtil.createBufferedWriter( matrix_name ) );
 386             }
 387             catch ( final IOException e ) {
 388                 e.printStackTrace();
 389             }
 390         }
 391         final NeighborJoiningF nj = NeighborJoiningF.createInstance( false, 5 );
 392         final Phylogeny phy = nj.execute( m );
 393         return phy;
 394     }
 395
 396     private final Phylogeny pi( final String matrix ) {
 397         final Phylogeny master_phy = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, _msa, true, matrix );
 398         final int seed = 15;
 399         final int n = 100;
 400         final ResampleableMsa resampleable_msa = new ResampleableMsa( _msa );
 401         final int[][] resampled_column_positions = BootstrapResampler.createResampledColumnPositions( _msa.getLength(),
 402                                                                                                       n,
 403                                                                                                       seed );
 404         final Phylogeny[] eval_phys = new Phylogeny[ n ];
 405         for( int i = 0; i < n; ++i ) {
 406             resampleable_msa.resample( resampled_column_positions[ i ] );
 407             eval_phys[ i ] = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, resampleable_msa, false, null );
 408         }
 409         ConfidenceAssessor.evaluate( "bootstrap", eval_phys, master_phy, true, 1 );
 410         PhylogenyMethods.extractFastaInformation( master_phy );
 411         return master_phy;
 412     }
 413
 414     private final void printMsaProperties( final String id, final MsaProperties msa_properties ) {
 415         System.out.print( ForesterUtil.pad( id, _longest_id_length, ' ', false ) );
 416         System.out.print( "\t" );
 417         final StringBuilder sb = msaPropertiesAsSB( msa_properties );
 418         System.out.print( sb );
 419         System.out.print( "\t" );
 420     }
 421
 422     private final StringBuilder msaPropertiesAsSB( final MsaProperties msa_properties ) {
 423         final StringBuilder sb = new StringBuilder();
 424         sb.append( msa_properties.getNumberOfSequences() );
 425         sb.append( "\t" );
 426         sb.append( msa_properties.getLength() );
 427         sb.append( "\t" );
 428         sb.append( NF_4.format( msa_properties.getGapRatio() ) );
 429         if ( _report_aln_mean_identity /*msa_properties.getAverageIdentityRatio() >= 0*/) {
 430             sb.append( "\t" );
 431             sb.append( NF_4.format( msa_properties.getAverageIdentityRatio() ) );
 432         }
 433         return sb;
 434     }
 435
 436     final private void printMsaStatsWriteOutfileAndRealign( final boolean realign, final String id )
 437             throws IOException, InterruptedException {
 438         if ( realign ) {
 439             realignWithMafft();
 440         }
 441         final MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
 442         printMsaProperties( id, msa_prop );
 443         final String s = writeOutfile();
 444         System.out.print( "-> " + s + ( realign ? "\t(realigned)" : "" ) );
 445     }
 446
 447     final private void realignWithMafft() throws IOException, InterruptedException {
 448         final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
 449         final List<String> opts = new ArrayList<String>();
 450         for( final String o : _maffts_opts.split( "\\s" ) ) {
 451             opts.add( o );
 452         }
 453         _msa = DeleteableMsa.createInstance( mafft.infer( _msa.asSequenceList(), opts ) );
 454     }
 455
 456     final private void removeGapColumns() {
 457         _msa.deleteGapOnlyColumns();
 458     }
 459
 460     final private static void writeMsa( final Msa msa, final String outfile, final MSA_FORMAT format )
 461             throws IOException {
 462         final Writer w = ForesterUtil.createBufferedWriter( outfile );
 463         msa.write( w, format );
 464         w.close();
 465     }
 466
 467     private final String writeOutfile() throws IOException {
 468         final String s = writeMsa( _out_file_base, MSA_FORMAT.PHYLIP, ".aln" );
 469         //writeMsa( _out_file_base, MSA_FORMAT.FASTA, ".fasta" );
 470         return s;
 471     }
 472
 473     // Returns null if not path found.
 474     final public static String guessPathToMafft() {
 475         String path;
 476         if ( ForesterUtil.OS_NAME.toLowerCase().indexOf( "win" ) >= 0 ) {
 477             path = "C:\\Program Files\\mafft-win\\mafft.bat";
 478             if ( MsaInferrer.isInstalled( path ) ) {
 479                 return path;
 480             }
 481         }
 482         path = "/home/czmasek/SOFTWARE/MSA/MAFFT/mafft-7.130-without-extensions/scripts/mafft";
 483         if ( MsaInferrer.isInstalled( path ) ) {
 484             return path;
 485         }
 486         path = "/usr/local/bin/mafft";
 487         if ( MsaInferrer.isInstalled( path ) ) {
 488             return path;
 489         }
 490         path = "/usr/bin/mafft";
 491         if ( MsaInferrer.isInstalled( path ) ) {
 492             return path;
 493         }
 494         path = "/bin/mafft";
 495         if ( MsaInferrer.isInstalled( path ) ) {
 496             return path;
 497         }
 498         path = "mafft";
 499         if ( MsaInferrer.isInstalled( path ) ) {
 500             return path;
 501         }
 502         return null;
 503     }
 504
 505     private final void printTableHeader() {
 506         System.out.print( ForesterUtil.pad( "Id", _longest_id_length, ' ', false ) );
 507         System.out.print( "\t" );
 508         System.out.print( "Seqs" );
 509         System.out.print( "\t" );
 510         System.out.print( "Length" );
 511         System.out.print( "\t" );
 512         System.out.print( "Gaps" );
 513         System.out.print( "\t" );
 514         if ( _report_aln_mean_identity ) {
 515             System.out.print( "MSA qual" );
 516             System.out.print( "\t" );
 517         }
 518         System.out.println();
 519     }
 520 }