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