forester/java/src/org/forester/msa_compactor/MsaCompactor.java

   1
   2 package org.forester.msa_compactor;
   3
   4 import java.io.File;
   5 import java.io.IOException;
   6 import java.io.Writer;
   7 import java.math.RoundingMode;
   8 import java.text.DecimalFormat;
   9 import java.text.NumberFormat;
  10 import java.util.ArrayList;
  11 import java.util.Arrays;
  12 import java.util.List;
  13 import java.util.SortedSet;
  14 import java.util.TreeSet;
  15
  16 import org.forester.evoinference.distance.NeighborJoiningF;
  17 import org.forester.evoinference.distance.PairwiseDistanceCalculator;
  18 import org.forester.evoinference.distance.PairwiseDistanceCalculator.PWD_DISTANCE_METHOD;
  19 import org.forester.evoinference.matrix.distance.BasicSymmetricalDistanceMatrix;
  20 import org.forester.evoinference.tools.BootstrapResampler;
  21 import org.forester.msa.BasicMsa;
  22 import org.forester.msa.Mafft;
  23 import org.forester.msa.Msa;
  24 import org.forester.msa.Msa.MSA_FORMAT;
  25 import org.forester.msa.MsaInferrer;
  26 import org.forester.msa.MsaMethods;
  27 import org.forester.msa.ResampleableMsa;
  28 import org.forester.phylogeny.Phylogeny;
  29 import org.forester.phylogeny.PhylogenyMethods;
  30 import org.forester.sequence.Sequence;
  31 import org.forester.tools.ConfidenceAssessor;
  32 import org.forester.util.BasicDescriptiveStatistics;
  33 import org.forester.util.DescriptiveStatistics;
  34 import org.forester.util.ForesterUtil;
  35
  36 public class MsaCompactor {
  37
  38     final private static NumberFormat NF_3    = new DecimalFormat( "#.###" );
  39     final private static NumberFormat NF_4    = new DecimalFormat( "#.####" );
  40     private static final boolean      VERBOSE = true;
  41     private Msa                       _msa;
  42     private final SortedSet<String>   _removed_seq_ids;
  43     private String                    _path_to_mafft;
  44     static {
  45         NF_4.setRoundingMode( RoundingMode.HALF_UP );
  46         NF_3.setRoundingMode( RoundingMode.HALF_UP );
  47     }
  48
  49     private MsaCompactor( final Msa msa ) {
  50         _msa = msa;
  51         _removed_seq_ids = new TreeSet<String>();
  52     }
  53
  54     final public Msa getMsa() {
  55         return _msa;
  56     }
  57
  58     final public SortedSet<String> getRemovedSeqIds() {
  59         return _removed_seq_ids;
  60     }
  61
  62     final public void writeMsa( final File outfile, final MSA_FORMAT format, final String suffix ) throws IOException {
  63         final Double gr = MsaMethods.calcGapRatio( _msa );
  64         writeMsa( outfile + "_" + _msa.getNumberOfSequences() + "_" + _msa.getLength() + "_"
  65                           + ForesterUtil.roundToInt( gr * 100 ) + suffix,
  66                   format );
  67     }
  68
  69     final int calcNonGapResidues( final Sequence seq ) {
  70         int ng = 0;
  71         for( int i = 0; i < seq.getLength(); ++i ) {
  72             if ( !seq.isGapAt( i ) ) {
  73                 ++ng;
  74             }
  75         }
  76         return ng;
  77     }
  78
  79     private final GapContribution[] calcGapContribtions( final boolean normalize_for_effective_seq_length ) {
  80         final double gappiness[] = calcGappiness();
  81         final GapContribution stats[] = new GapContribution[ _msa.getNumberOfSequences() ];
  82         for( int row = 0; row < _msa.getNumberOfSequences(); ++row ) {
  83             stats[ row ] = new GapContribution( _msa.getIdentifier( row ) );
  84             for( int col = 0; col < _msa.getLength(); ++col ) {
  85                 if ( !_msa.isGapAt( row, col ) ) {
  86                     stats[ row ].addToValue( gappiness[ col ] );
  87                 }
  88             }
  89             if ( normalize_for_effective_seq_length ) {
  90                 stats[ row ].divideValue( calcNonGapResidues( _msa.getSequence( row ) ) );
  91             }
  92             else {
  93                 stats[ row ].divideValue( _msa.getLength() );
  94             }
  95         }
  96         return stats;
  97     }
  98
  99     final private GapContribution[] calcGapContribtionsStats( final boolean norm ) {
 100         final GapContribution stats[] = calcGapContribtions( norm );
 101         Arrays.sort( stats );
 102         // for( final GapContribution stat : stats ) {
 103         //  final StringBuilder sb = new StringBuilder();
 104         //  sb.append( stat.getId() );
 105         //  sb.append( "\t" );
 106         //  sb.append( NF_4.format( stat.getValue() ) );
 107         //  sb.append( "\t" );
 108         //            sb.append( NF_4.format( stat.median() ) );
 109         //            sb.append( "\t" );
 110         //            sb.append( NF_4.format( stat.getMin() ) );
 111         //            sb.append( "\t" );
 112         //            sb.append( NF_4.format( stat.getMax() ) );
 113         //sb.append( "\t" );
 114         //System.out.println( sb );
 115         // }
 116         return stats;
 117     }
 118
 119     private static DescriptiveStatistics calculateIdentityRatio( final int from, final int to, final Msa msa ) {
 120         final DescriptiveStatistics stats = new BasicDescriptiveStatistics();
 121         for( int c = from; c <= to; ++c ) {
 122             stats.addValue( MsaMethods.calculateIdentityRatio( msa, c ) );
 123         }
 124         return stats;
 125     }
 126
 127     private final double[] calcGappiness() {
 128         final int l = _msa.getLength();
 129         final double gappiness[] = new double[ l ];
 130         final int seqs = _msa.getNumberOfSequences();
 131         for( int i = 0; i < l; ++i ) {
 132             gappiness[ i ] = ( double ) MsaMethods.calcGapSumPerColumn( _msa, i ) / seqs;
 133         }
 134         return gappiness;
 135     }
 136
 137     // Returns null if not path found.
 138     final public static String guessPathToMafft() {
 139         String path;
 140         if ( ForesterUtil.OS_NAME.toLowerCase().indexOf( "win" ) >= 0 ) {
 141             path = "C:\\Program Files\\mafft-win\\mafft.bat";
 142             if ( MsaInferrer.isInstalled( path ) ) {
 143                 return path;
 144             }
 145         }
 146         path = "/usr/local/bin/mafft";
 147         if ( MsaInferrer.isInstalled( path ) ) {
 148             return path;
 149         }
 150         path = "/usr/bin/mafft";
 151         if ( MsaInferrer.isInstalled( path ) ) {
 152             return path;
 153         }
 154         path = "/bin/mafft";
 155         if ( MsaInferrer.isInstalled( path ) ) {
 156             return path;
 157         }
 158         path = "mafft";
 159         if ( MsaInferrer.isInstalled( path ) ) {
 160             return path;
 161         }
 162         return null;
 163     }
 164
 165     final private void mafft() throws IOException, InterruptedException {
 166         //  final MsaInferrer mafft = Mafft
 167         //       .createInstance( "/home/czmasek/SOFTWARE/MSA/MAFFT/mafft-7.130-without-extensions/scripts/mafft" );
 168         final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
 169         final List<String> opts = new ArrayList<String>();
 170         opts.add( "--maxiterate" );
 171         opts.add( "1000" );
 172         opts.add( "--localpair" );
 173         opts.add( "--quiet" );
 174         _msa = mafft.infer( _msa.asSequenceList(), opts );
 175     }
 176
 177     private StringBuilder msaStatsAsSB() {
 178         final StringBuilder sb = new StringBuilder();
 179         sb.append( _msa.getNumberOfSequences() );
 180         sb.append( "\t" );
 181         sb.append( _msa.getLength() );
 182         sb.append( "\t" );
 183         sb.append( NF_3.format( MsaMethods.calcGapRatio( _msa ) ) );
 184         sb.append( "\t" );
 185         sb.append( NF_3.format( calculateIdentityRatio( 0, _msa.getLength() - 1, _msa ).arithmeticMean() ) );
 186         return sb;
 187     }
 188
 189     final private void removeGapColumns() {
 190         _msa = MsaMethods.createInstance().removeGapColumns( 1, 0, _msa );
 191     }
 192
 193     final private void removeViaGapAverage( final double mean_gapiness,
 194                                             final int step,
 195                                             final boolean realign,
 196                                             final File outfile,
 197                                             final int minimal_effective_length ) throws IOException,
 198             InterruptedException {
 199         if ( step < 1 ) {
 200             throw new IllegalArgumentException( "step cannot be less than 1" );
 201         }
 202         if ( mean_gapiness < 0 ) {
 203             throw new IllegalArgumentException( "target average gap ratio cannot be less than 0" );
 204         }
 205         if ( VERBOSE ) {
 206             System.out.println( "orig: " + msaStatsAsSB() );
 207         }
 208         if ( minimal_effective_length > 1 ) {
 209             _msa = MsaMethods.removeSequencesByMinimalLength( _msa, minimal_effective_length );
 210             if ( VERBOSE ) {
 211                 System.out.println( "short seq removal: " + msaStatsAsSB() );
 212             }
 213         }
 214         int counter = step;
 215         double gr;
 216         do {
 217             removeWorstOffenders( step, 1, false, false );
 218             if ( realign ) {
 219                 mafft();
 220             }
 221             gr = MsaMethods.calcGapRatio( _msa );
 222             if ( VERBOSE ) {
 223                 System.out.println( counter + ": " + msaStatsAsSB() );
 224             }
 225             //   write( outfile, gr );
 226             counter += step;
 227         } while ( gr > mean_gapiness );
 228         if ( VERBOSE ) {
 229             System.out.println( "final: " + msaStatsAsSB() );
 230         }
 231     }
 232
 233     final private void removeViaLength( final int length, final int step, final boolean realign ) throws IOException,
 234             InterruptedException {
 235         if ( step < 1 ) {
 236             throw new IllegalArgumentException( "step cannot be less than 1" );
 237         }
 238         if ( length < 11 ) {
 239             throw new IllegalArgumentException( "target length cannot be less than 1" );
 240         }
 241         if ( VERBOSE ) {
 242             System.out.println( "orig: " + msaStatsAsSB() );
 243         }
 244         int counter = step;
 245         while ( _msa.getLength() > length ) {
 246             removeWorstOffenders( step, 1, false, false );
 247             if ( realign ) {
 248                 mafft();
 249             }
 250             if ( VERBOSE ) {
 251                 System.out.println( counter + ": " + msaStatsAsSB() );
 252             }
 253             counter += step;
 254         }
 255     }
 256
 257     Phylogeny pi( final String matrix ) {
 258         final Phylogeny master_phy = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, _msa, true, matrix );
 259         final int seed = 15;
 260         final int n = 100;
 261         final ResampleableMsa resampleable_msa = new ResampleableMsa( ( BasicMsa ) _msa );
 262         final int[][] resampled_column_positions = BootstrapResampler.createResampledColumnPositions( _msa.getLength(),
 263                                                                                                       n,
 264                                                                                                       seed );
 265         final Phylogeny[] eval_phys = new Phylogeny[ n ];
 266         for( int i = 0; i < n; ++i ) {
 267             resampleable_msa.resample( resampled_column_positions[ i ] );
 268             eval_phys[ i ] = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, resampleable_msa, false, null );
 269         }
 270         ConfidenceAssessor.evaluate( "bootstrap", eval_phys, master_phy, true, 1 );
 271         PhylogenyMethods.extractFastaInformation( master_phy );
 272         return master_phy;
 273     }
 274
 275     private Phylogeny inferNJphylogeny( final PWD_DISTANCE_METHOD pwd_distance_method,
 276                                         final Msa msa,
 277                                         final boolean write_matrix,
 278                                         final String matrix_name ) {
 279         BasicSymmetricalDistanceMatrix m = null;
 280         switch ( pwd_distance_method ) {
 281             case KIMURA_DISTANCE:
 282                 m = PairwiseDistanceCalculator.calcKimuraDistances( msa );
 283                 break;
 284             case POISSON_DISTANCE:
 285                 m = PairwiseDistanceCalculator.calcPoissonDistances( msa );
 286                 break;
 287             case FRACTIONAL_DISSIMILARITY:
 288                 m = PairwiseDistanceCalculator.calcFractionalDissimilarities( msa );
 289                 break;
 290             default:
 291                 throw new IllegalArgumentException( "invalid pwd method" );
 292         }
 293         if ( write_matrix ) {
 294             try {
 295                 m.write( ForesterUtil.createBufferedWriter( matrix_name ) );
 296             }
 297             catch ( final IOException e ) {
 298                 // TODO Auto-generated catch block
 299                 e.printStackTrace();
 300             }
 301         }
 302         final NeighborJoiningF nj = NeighborJoiningF.createInstance( false, 5 );
 303         final Phylogeny phy = nj.execute( m );
 304         return phy;
 305     }
 306
 307     final private void removeWorstOffenders( final int to_remove,
 308                                              final int step,
 309                                              final boolean realign,
 310                                              final boolean norm ) throws IOException, InterruptedException {
 311         //final Phylogeny a = pi( "a.pwd" );
 312         //Archaeopteryx.createApplication( a );
 313         final GapContribution stats[] = calcGapContribtionsStats( norm );
 314         final List<String> to_remove_ids = new ArrayList<String>();
 315         for( int j = 0; j < to_remove; ++j ) {
 316             to_remove_ids.add( stats[ j ].getId() );
 317             _removed_seq_ids.add( stats[ j ].getId() );
 318         }
 319         //TODO if verbose/interactive
 320         for( final String id : to_remove_ids ) {
 321             _msa = MsaMethods.removeSequence( _msa, id );
 322             removeGapColumns();
 323             //System.out.print( id );
 324             System.out.print( ForesterUtil.pad( id, 20, ' ', false ) );
 325             System.out.print( "\t" );
 326             final StringBuilder sb = msaStatsAsSB();
 327             System.out.println( sb );
 328         }
 329         //TODO else:
 330         //_msa = MsaMethods.removeSequences( _msa, to_remove_ids );
 331         //removeGapColumns();
 332         if ( realign ) {
 333             mafft();
 334         }
 335         //final Phylogeny b = pi( "b.pwd" );
 336         //Archaeopteryx.createApplication( b );
 337     }
 338
 339     final private void writeMsa( final String outfile, final MSA_FORMAT format ) throws IOException {
 340         final Writer w = ForesterUtil.createBufferedWriter( outfile );
 341         _msa.write( w, format );
 342         w.close();
 343     }
 344
 345     public final static MsaCompactor reduceGapAverage( final Msa msa,
 346                                                        final double max_gap_average,
 347                                                        final int step,
 348                                                        final boolean realign,
 349                                                        final File out,
 350                                                        final int minimal_effective_length,
 351                                                        final String path_to_mafft ) throws IOException,
 352             InterruptedException {
 353         final MsaCompactor mc = new MsaCompactor( msa );
 354         if ( realign ) {
 355             mc.setPathToMafft( path_to_mafft );
 356         }
 357         mc.removeViaGapAverage( max_gap_average, step, realign, out, minimal_effective_length );
 358         return mc;
 359     }
 360
 361     public final static MsaCompactor reduceLength( final Msa msa,
 362                                                    final int length,
 363                                                    final int step,
 364                                                    final boolean realign,
 365                                                    final String path_to_mafft ) throws IOException,
 366             InterruptedException {
 367         final MsaCompactor mc = new MsaCompactor( msa );
 368         if ( realign ) {
 369             mc.setPathToMafft( path_to_mafft );
 370         }
 371         mc.removeViaLength( length, step, realign );
 372         return mc;
 373     }
 374
 375     public final static MsaCompactor removeWorstOffenders( final Msa msa,
 376                                                            final int worst_offenders_to_remove,
 377                                                            final boolean realign,
 378                                                            final boolean norm,
 379                                                            final String path_to_mafft ) throws IOException,
 380             InterruptedException {
 381         final MsaCompactor mc = new MsaCompactor( msa );
 382         if ( realign ) {
 383             mc.setPathToMafft( path_to_mafft );
 384         }
 385         mc.removeWorstOffenders( worst_offenders_to_remove, 1, realign, norm );
 386         return mc;
 387     }
 388
 389     private void setPathToMafft( final String path_to_mafft ) {
 390         _path_to_mafft = path_to_mafft;
 391     }
 392 }