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