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.Comparator;
  13 import java.util.List;
  14 import java.util.SortedSet;
  15 import java.util.TreeSet;
  16
  17 import org.forester.msa.Mafft;
  18 import org.forester.msa.Msa;
  19 import org.forester.msa.Msa.MSA_FORMAT;
  20 import org.forester.msa.MsaInferrer;
  21 import org.forester.msa.MsaMethods;
  22 import org.forester.sequence.Sequence;
  23 import org.forester.util.BasicDescriptiveStatistics;
  24 import org.forester.util.DescriptiveStatistics;
  25 import org.forester.util.ForesterUtil;
  26
  27 public class MsaCompactor {
  28
  29     final private static NumberFormat NF_3    = new DecimalFormat( "#.###" );
  30     final private static NumberFormat NF_4    = new DecimalFormat( "#.####" );
  31     private static final boolean      VERBOSE = true;
  32     private Msa                       _msa;
  33     private final SortedSet<String>   _removed_seq_ids;
  34     static {
  35         NF_4.setRoundingMode( RoundingMode.HALF_UP );
  36         NF_3.setRoundingMode( RoundingMode.HALF_UP );
  37     }
  38
  39     private MsaCompactor( final Msa msa ) {
  40         _msa = msa;
  41         _removed_seq_ids = new TreeSet<String>();
  42     }
  43
  44     final public Msa getMsa() {
  45         return _msa;
  46     }
  47
  48     final public SortedSet<String> getRemovedSeqIds() {
  49         return _removed_seq_ids;
  50     }
  51
  52     final public void writeMsa( final File outfile, final MSA_FORMAT format, final String suffix ) throws IOException {
  53         final Double gr = MsaMethods.calcGapRatio( _msa );
  54         writeMsa( outfile + "_" + _msa.getNumberOfSequences() + "_" + _msa.getLength() + "_"
  55                           + ForesterUtil.roundToInt( gr * 100 ) + suffix,
  56                   format );
  57     }
  58
  59     final int calcNonGapResidues( final Sequence seq ) {
  60         int ng = 0;
  61         for( int i = 0; i < seq.getLength(); ++i ) {
  62             if ( !seq.isGapAt( i ) ) {
  63                 ++ng;
  64             }
  65         }
  66         return ng;
  67     }
  68
  69     private final DescriptiveStatistics[] calcGapContribtionsX( final boolean normalize_for_effective_seq_length ) {
  70         final double gappiness[] = calcGappiness();
  71         final DescriptiveStatistics stats[] = new DescriptiveStatistics[ _msa.getNumberOfSequences() ];
  72         for( int row = 0; row < _msa.getNumberOfSequences(); ++row ) {
  73             stats[ row ] = new BasicDescriptiveStatistics( _msa.getIdentifier( row ) );
  74             final double l = calculateEffectiveLengthRatio( row );
  75             for( int col = 0; col < _msa.getLength(); ++col ) {
  76                 if ( !_msa.isGapAt( row, col ) ) {
  77                     if ( normalize_for_effective_seq_length ) {
  78                         stats[ row ].addValue( gappiness[ col ] / l );
  79                     }
  80                     else {
  81                         stats[ row ].addValue( gappiness[ col ] );
  82                     }
  83                 }
  84             }
  85         }
  86         return stats;
  87     }
  88
  89     private final GapContribution[] calcGapContribtions( final boolean normalize_for_effective_seq_length ) {
  90         final double gappiness[] = calcGappiness();
  91         final GapContribution stats[] = new GapContribution[ _msa.getNumberOfSequences() ];
  92         for( int row = 0; row < _msa.getNumberOfSequences(); ++row ) {
  93             stats[ row ] = new GapContribution( _msa.getIdentifier( row ) );
  94             for( int col = 0; col < _msa.getLength(); ++col ) {
  95                 if ( !_msa.isGapAt( row, col ) ) {
  96                     stats[ row ].addToValue( gappiness[ col ] );
  97                 }
  98             }
  99             if ( normalize_for_effective_seq_length ) {
 100                 stats[ row ].divideValue( calculateEffectiveLengthRatio( row ) );
 101             }
 102             else {
 103                 //
 104             }
 105         }
 106         return stats;
 107     }
 108
 109     final private GapContribution[] calcGapContribtionsStats( final boolean norm ) {
 110         final GapContribution stats[] = calcGapContribtions( norm );
 111         Arrays.sort( stats );
 112         for( final GapContribution stat : stats ) {
 113             final StringBuilder sb = new StringBuilder();
 114             sb.append( stat.getId() );
 115             sb.append( "\t" );
 116             sb.append( NF_4.format( stat.getValue() ) );
 117             sb.append( "\t" );
 118             //            sb.append( NF_4.format( stat.median() ) );
 119             //            sb.append( "\t" );
 120             //            sb.append( NF_4.format( stat.getMin() ) );
 121             //            sb.append( "\t" );
 122             //            sb.append( NF_4.format( stat.getMax() ) );
 123             //sb.append( "\t" );
 124             System.out.println( sb );
 125         }
 126         return stats;
 127     }
 128
 129     private final double[] calcGappiness() {
 130         final int l = _msa.getLength();
 131         final double gappiness[] = new double[ l ];
 132         final int seqs = _msa.getNumberOfSequences();
 133         for( int i = 0; i < l; ++i ) {
 134             gappiness[ i ] = ( double ) MsaMethods.calcGapSumPerColumn( _msa, i ) / seqs;
 135         }
 136         return gappiness;
 137     }
 138
 139     private double calculateEffectiveLengthRatio( final int row ) {
 140         return ( double ) calcNonGapResidues( _msa.getSequence( row ) ) / _msa.getLength();
 141     }
 142
 143     final private void mafft() throws IOException, InterruptedException {
 144         final MsaInferrer mafft = Mafft
 145                 .createInstance( "/home/czmasek/SOFTWARE/MSA/MAFFT/mafft-7.130-without-extensions/scripts/mafft" );
 146         final List<String> opts = new ArrayList<String>();
 147         // opts.add( "--maxiterate" );
 148         // opts.add( "1000" );
 149         // opts.add( "--localpair" );
 150         opts.add( "--quiet" );
 151         _msa = mafft.infer( _msa.asSequenceList(), opts );
 152     }
 153
 154     private StringBuilder msaStatsAsSB() {
 155         final StringBuilder sb = new StringBuilder();
 156         sb.append( _msa.getNumberOfSequences() );
 157         sb.append( "\t" );
 158         sb.append( _msa.getLength() );
 159         sb.append( "\t" );
 160         sb.append( NF_3.format( MsaMethods.calcGapRatio( _msa ) ) );
 161         return sb;
 162     }
 163
 164     final private void removeGapColumns() {
 165         _msa = MsaMethods.createInstance().removeGapColumns( 1, 0, _msa );
 166     }
 167
 168     final private void removeViaGapAverage( final double mean_gapiness,
 169                                             final int step,
 170                                             final boolean realign,
 171                                             final File outfile,
 172                                             final int minimal_effective_length ) throws IOException,
 173             InterruptedException {
 174         if ( step < 1 ) {
 175             throw new IllegalArgumentException( "step cannot be less than 1" );
 176         }
 177         if ( mean_gapiness < 0 ) {
 178             throw new IllegalArgumentException( "target average gap ratio cannot be less than 0" );
 179         }
 180         if ( VERBOSE ) {
 181             System.out.println( "orig: " + msaStatsAsSB() );
 182         }
 183         if ( minimal_effective_length > 1 ) {
 184             _msa = MsaMethods.removeSequencesByMinimalLength( _msa, minimal_effective_length );
 185             if ( VERBOSE ) {
 186                 System.out.println( "short seq removal: " + msaStatsAsSB() );
 187             }
 188         }
 189         int counter = step;
 190         double gr;
 191         do {
 192             removeWorstOffenders( step, 1, false, false );
 193             if ( realign ) {
 194                 mafft();
 195             }
 196             gr = MsaMethods.calcGapRatio( _msa );
 197             if ( VERBOSE ) {
 198                 System.out.println( counter + ": " + msaStatsAsSB() );
 199             }
 200             //   write( outfile, gr );
 201             counter += step;
 202         } while ( gr > mean_gapiness );
 203         if ( VERBOSE ) {
 204             System.out.println( "final: " + msaStatsAsSB() );
 205         }
 206     }
 207
 208     final private void removeViaLength( final int length, final int step, final boolean realign ) throws IOException,
 209             InterruptedException {
 210         if ( step < 1 ) {
 211             throw new IllegalArgumentException( "step cannot be less than 1" );
 212         }
 213         if ( length < 11 ) {
 214             throw new IllegalArgumentException( "target length cannot be less than 1" );
 215         }
 216         if ( VERBOSE ) {
 217             System.out.println( "orig: " + msaStatsAsSB() );
 218         }
 219         int counter = step;
 220         while ( _msa.getLength() > length ) {
 221             removeWorstOffenders( step, 1, false, false );
 222             if ( realign ) {
 223                 mafft();
 224             }
 225             if ( VERBOSE ) {
 226                 System.out.println( counter + ": " + msaStatsAsSB() );
 227             }
 228             counter += step;
 229         }
 230     }
 231
 232     final private void removeWorstOffenders( final int to_remove,
 233                                              final int step,
 234                                              final boolean realign,
 235                                              final boolean norm ) throws IOException, InterruptedException {
 236         final GapContribution stats[] = calcGapContribtionsStats( norm );
 237         final List<String> to_remove_ids = new ArrayList<String>();
 238         for( int j = 0; j < to_remove; ++j ) {
 239             to_remove_ids.add( stats[ j ].getId() );
 240             _removed_seq_ids.add( stats[ j ].getId() );
 241         }
 242         //TODO if verbose/interactve
 243         for( final String id : to_remove_ids ) {
 244             _msa = MsaMethods.removeSequence( _msa, id );
 245             removeGapColumns();
 246             System.out.print( id );
 247             System.out.print( "\t" );
 248             final StringBuilder sb = msaStatsAsSB();
 249             System.out.println( sb );
 250         }
 251         //TODO else:
 252         //_msa = MsaMethods.removeSequences( _msa, to_remove_ids );
 253         //removeGapColumns();
 254         if ( realign ) {
 255             mafft();
 256         }
 257     }
 258
 259     final private void writeMsa( final String outfile, final MSA_FORMAT format ) throws IOException {
 260         final Writer w = ForesterUtil.createBufferedWriter( outfile );
 261         _msa.write( w, format );
 262         w.close();
 263     }
 264
 265     public final static MsaCompactor reduceGapAverage( final Msa msa,
 266                                                        final double max_gap_average,
 267                                                        final int step,
 268                                                        final boolean realign,
 269                                                        final File out,
 270                                                        final int minimal_effective_length ) throws IOException,
 271             InterruptedException {
 272         final MsaCompactor mc = new MsaCompactor( msa );
 273         mc.removeViaGapAverage( max_gap_average, step, realign, out, minimal_effective_length );
 274         return mc;
 275     }
 276
 277     public final static MsaCompactor reduceLength( final Msa msa,
 278                                                    final int length,
 279                                                    final int step,
 280                                                    final boolean realign ) throws IOException, InterruptedException {
 281         final MsaCompactor mc = new MsaCompactor( msa );
 282         mc.removeViaLength( length, step, realign );
 283         return mc;
 284     }
 285
 286     public final static MsaCompactor removeWorstOffenders( final Msa msa,
 287                                                            final int worst_offenders_to_remove,
 288                                                            final boolean realign,
 289                                                            final boolean norm ) throws IOException,
 290             InterruptedException {
 291         final MsaCompactor mc = new MsaCompactor( msa );
 292         mc.removeWorstOffenders( worst_offenders_to_remove, 1, realign, norm );
 293         return mc;
 294     }
 295
 296     public static enum SORT_BY {
 297         MAX, MEAN, MEDIAN;
 298     }
 299
 300     final static class DescriptiveStatisticsComparator implements Comparator<DescriptiveStatistics> {
 301
 302         final private boolean _ascending;
 303         final private SORT_BY _sort_by;
 304
 305         public DescriptiveStatisticsComparator( final boolean ascending, final SORT_BY sort_by ) {
 306             _ascending = ascending;
 307             _sort_by = sort_by;
 308         }
 309
 310         @Override
 311         public final int compare( final DescriptiveStatistics s0, final DescriptiveStatistics s1 ) {
 312             switch ( _sort_by ) {
 313                 case MAX:
 314                     if ( s0.getMax() < s1.getMax() ) {
 315                         return _ascending ? -1 : 1;
 316                     }
 317                     else if ( s0.getMax() > s1.getMax() ) {
 318                         return _ascending ? 1 : -1;
 319                     }
 320                     return 0;
 321                 case MEAN:
 322                     if ( s0.arithmeticMean() < s1.arithmeticMean() ) {
 323                         return _ascending ? -1 : 1;
 324                     }
 325                     else if ( s0.arithmeticMean() > s1.arithmeticMean() ) {
 326                         return _ascending ? 1 : -1;
 327                     }
 328                     return 0;
 329                 case MEDIAN:
 330                     if ( s0.median() < s1.median() ) {
 331                         return _ascending ? -1 : 1;
 332                     }
 333                     else if ( s0.median() > s1.median() ) {
 334                         return _ascending ? 1 : -1;
 335                     }
 336                     return 0;
 337                 default:
 338                     return 0;
 339             }
 340         }
 341     }
 342 }