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