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