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