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.NeighborJoining;
  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     static {
  43         NF_4.setRoundingMode( RoundingMode.HALF_UP );
  44         NF_3.setRoundingMode( RoundingMode.HALF_UP );
  45     }
  46
  47     private MsaCompactor( final Msa msa ) {
  48         _msa = msa;
  49         _removed_seq_ids = new TreeSet<String>();
  50     }
  51
  52     final public Msa getMsa() {
  53         return _msa;
  54     }
  55
  56     final public SortedSet<String> getRemovedSeqIds() {
  57         return _removed_seq_ids;
  58     }
  59
  60     final public void writeMsa( final File outfile, final MSA_FORMAT format, final String suffix ) throws IOException {
  61         final Double gr = MsaMethods.calcGapRatio( _msa );
  62         writeMsa( outfile + "_" + _msa.getNumberOfSequences() + "_" + _msa.getLength() + "_"
  63                           + ForesterUtil.roundToInt( gr * 100 ) + suffix,
  64                   format );
  65     }
  66
  67     final int calcNonGapResidues( final Sequence seq ) {
  68         int ng = 0;
  69         for( int i = 0; i < seq.getLength(); ++i ) {
  70             if ( !seq.isGapAt( i ) ) {
  71                 ++ng;
  72             }
  73         }
  74         return ng;
  75     }
  76
  77     private final GapContribution[] calcGapContribtions( final boolean normalize_for_effective_seq_length ) {
  78         final double gappiness[] = calcGappiness();
  79         final GapContribution stats[] = new GapContribution[ _msa.getNumberOfSequences() ];
  80         for( int row = 0; row < _msa.getNumberOfSequences(); ++row ) {
  81             stats[ row ] = new GapContribution( _msa.getIdentifier( row ) );
  82             for( int col = 0; col < _msa.getLength(); ++col ) {
  83                 if ( !_msa.isGapAt( row, col ) ) {
  84                     stats[ row ].addToValue( gappiness[ col ] );
  85                 }
  86             }
  87             if ( normalize_for_effective_seq_length ) {
  88                 stats[ row ].divideValue( calcNonGapResidues( _msa.getSequence( row ) ) );
  89             }
  90             else {
  91                 stats[ row ].divideValue( _msa.getLength() );
  92             }
  93         }
  94         return stats;
  95     }
  96
  97     final private GapContribution[] calcGapContribtionsStats( final boolean norm ) {
  98         final GapContribution stats[] = calcGapContribtions( norm );
  99         Arrays.sort( stats );
 100         for( final GapContribution stat : stats ) {
 101             final StringBuilder sb = new StringBuilder();
 102             sb.append( stat.getId() );
 103             sb.append( "\t" );
 104             sb.append( NF_4.format( stat.getValue() ) );
 105             sb.append( "\t" );
 106             //            sb.append( NF_4.format( stat.median() ) );
 107             //            sb.append( "\t" );
 108             //            sb.append( NF_4.format( stat.getMin() ) );
 109             //            sb.append( "\t" );
 110             //            sb.append( NF_4.format( stat.getMax() ) );
 111             //sb.append( "\t" );
 112             System.out.println( sb );
 113         }
 114         return stats;
 115     }
 116
 117     private final double[] calcGappiness() {
 118         final int l = _msa.getLength();
 119         final double gappiness[] = new double[ l ];
 120         final int seqs = _msa.getNumberOfSequences();
 121         for( int i = 0; i < l; ++i ) {
 122             gappiness[ i ] = ( double ) MsaMethods.calcGapSumPerColumn( _msa, i ) / seqs;
 123         }
 124         return gappiness;
 125     }
 126
 127     final private void mafft() throws IOException, InterruptedException {
 128         final MsaInferrer mafft = Mafft
 129                 .createInstance( "/home/czmasek/SOFTWARE/MSA/MAFFT/mafft-7.130-without-extensions/scripts/mafft" );
 130         final List<String> opts = new ArrayList<String>();
 131         // opts.add( "--maxiterate" );
 132         // opts.add( "1000" );
 133         // opts.add( "--localpair" );
 134         opts.add( "--quiet" );
 135         _msa = mafft.infer( _msa.asSequenceList(), opts );
 136     }
 137
 138     private StringBuilder msaStatsAsSB() {
 139         final StringBuilder sb = new StringBuilder();
 140         sb.append( _msa.getNumberOfSequences() );
 141         sb.append( "\t" );
 142         sb.append( _msa.getLength() );
 143         sb.append( "\t" );
 144         sb.append( NF_3.format( MsaMethods.calcGapRatio( _msa ) ) );
 145         return sb;
 146     }
 147
 148     final private void removeGapColumns() {
 149         _msa = MsaMethods.createInstance().removeGapColumns( 1, 0, _msa );
 150     }
 151
 152     final private void removeViaGapAverage( final double mean_gapiness,
 153                                             final int step,
 154                                             final boolean realign,
 155                                             final File outfile,
 156                                             final int minimal_effective_length ) throws IOException,
 157             InterruptedException {
 158         if ( step < 1 ) {
 159             throw new IllegalArgumentException( "step cannot be less than 1" );
 160         }
 161         if ( mean_gapiness < 0 ) {
 162             throw new IllegalArgumentException( "target average gap ratio cannot be less than 0" );
 163         }
 164         if ( VERBOSE ) {
 165             System.out.println( "orig: " + msaStatsAsSB() );
 166         }
 167         if ( minimal_effective_length > 1 ) {
 168             _msa = MsaMethods.removeSequencesByMinimalLength( _msa, minimal_effective_length );
 169             if ( VERBOSE ) {
 170                 System.out.println( "short seq removal: " + msaStatsAsSB() );
 171             }
 172         }
 173         int counter = step;
 174         double gr;
 175         do {
 176             removeWorstOffenders( step, 1, false, false );
 177             if ( realign ) {
 178                 mafft();
 179             }
 180             gr = MsaMethods.calcGapRatio( _msa );
 181             if ( VERBOSE ) {
 182                 System.out.println( counter + ": " + msaStatsAsSB() );
 183             }
 184             //   write( outfile, gr );
 185             counter += step;
 186         } while ( gr > mean_gapiness );
 187         if ( VERBOSE ) {
 188             System.out.println( "final: " + msaStatsAsSB() );
 189         }
 190     }
 191
 192     final private void removeViaLength( final int length, final int step, final boolean realign ) throws IOException,
 193             InterruptedException {
 194         if ( step < 1 ) {
 195             throw new IllegalArgumentException( "step cannot be less than 1" );
 196         }
 197         if ( length < 11 ) {
 198             throw new IllegalArgumentException( "target length cannot be less than 1" );
 199         }
 200         if ( VERBOSE ) {
 201             System.out.println( "orig: " + msaStatsAsSB() );
 202         }
 203         int counter = step;
 204         while ( _msa.getLength() > length ) {
 205             removeWorstOffenders( step, 1, false, false );
 206             if ( realign ) {
 207                 mafft();
 208             }
 209             if ( VERBOSE ) {
 210                 System.out.println( counter + ": " + msaStatsAsSB() );
 211             }
 212             counter += step;
 213         }
 214     }
 215
 216     Phylogeny pi() {
 217         final Phylogeny master_phy = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, _msa );
 218         final int seed = 15;
 219         final int n = 100;
 220         final ResampleableMsa resampleable_msa = new ResampleableMsa( ( BasicMsa ) _msa );
 221         final int[][] resampled_column_positions = BootstrapResampler.createResampledColumnPositions( _msa.getLength(),
 222                                                                                                       n,
 223                                                                                                       seed );
 224         final Phylogeny[] eval_phys = new Phylogeny[ n ];
 225         for( int i = 0; i < n; ++i ) {
 226             resampleable_msa.resample( resampled_column_positions[ i ] );
 227             eval_phys[ i ] = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, resampleable_msa );
 228         }
 229         ConfidenceAssessor.evaluate( "bootstrap", eval_phys, master_phy, true, 1 );
 230         PhylogenyMethods.extractFastaInformation( master_phy );
 231         return master_phy;
 232     }
 233
 234     private Phylogeny inferNJphylogeny( PWD_DISTANCE_METHOD pwd_distance_method, final Msa msa ) {
 235         BasicSymmetricalDistanceMatrix m = null;
 236         switch ( pwd_distance_method ) {
 237             case KIMURA_DISTANCE:
 238                 m = PairwiseDistanceCalculator.calcKimuraDistances( msa );
 239                 break;
 240             case POISSON_DISTANCE:
 241                 m = PairwiseDistanceCalculator.calcPoissonDistances( msa );
 242                 break;
 243             case FRACTIONAL_DISSIMILARITY:
 244                 m = PairwiseDistanceCalculator.calcFractionalDissimilarities( msa );
 245                 break;
 246             default:
 247                 throw new IllegalArgumentException( "invalid pwd method" );
 248         }
 249         final NeighborJoining nj = NeighborJoining.createInstance();
 250         final Phylogeny phy = nj.execute( m );
 251         return phy;
 252     }
 253
 254     final private void removeWorstOffenders( final int to_remove,
 255                                              final int step,
 256                                              final boolean realign,
 257                                              final boolean norm ) throws IOException, InterruptedException {
 258         final Phylogeny a = pi();
 259         Archaeopteryx.createApplication( a );
 260         final GapContribution stats[] = calcGapContribtionsStats( norm );
 261         final List<String> to_remove_ids = new ArrayList<String>();
 262         for( int j = 0; j < to_remove; ++j ) {
 263             to_remove_ids.add( stats[ j ].getId() );
 264             _removed_seq_ids.add( stats[ j ].getId() );
 265         }
 266         //TODO if verbose/interactve
 267         for( final String id : to_remove_ids ) {
 268             _msa = MsaMethods.removeSequence( _msa, id );
 269             removeGapColumns();
 270             System.out.print( id );
 271             System.out.print( "\t" );
 272             final StringBuilder sb = msaStatsAsSB();
 273             System.out.println( sb );
 274         }
 275         //TODO else:
 276         //_msa = MsaMethods.removeSequences( _msa, to_remove_ids );
 277         //removeGapColumns();
 278         if ( realign ) {
 279             mafft();
 280         }
 281         final Phylogeny b = pi();
 282         Archaeopteryx.createApplication( b );
 283     }
 284
 285     final private void writeMsa( final String outfile, final MSA_FORMAT format ) throws IOException {
 286         final Writer w = ForesterUtil.createBufferedWriter( outfile );
 287         _msa.write( w, format );
 288         w.close();
 289     }
 290
 291     public final static MsaCompactor reduceGapAverage( final Msa msa,
 292                                                        final double max_gap_average,
 293                                                        final int step,
 294                                                        final boolean realign,
 295                                                        final File out,
 296                                                        final int minimal_effective_length ) throws IOException,
 297             InterruptedException {
 298         final MsaCompactor mc = new MsaCompactor( msa );
 299         mc.removeViaGapAverage( max_gap_average, step, realign, out, minimal_effective_length );
 300         return mc;
 301     }
 302
 303     public final static MsaCompactor reduceLength( final Msa msa,
 304                                                    final int length,
 305                                                    final int step,
 306                                                    final boolean realign ) throws IOException, InterruptedException {
 307         final MsaCompactor mc = new MsaCompactor( msa );
 308         mc.removeViaLength( length, step, realign );
 309         return mc;
 310     }
 311
 312     public final static MsaCompactor removeWorstOffenders( final Msa msa,
 313                                                            final int worst_offenders_to_remove,
 314                                                            final boolean realign,
 315                                                            final boolean norm ) throws IOException,
 316             InterruptedException {
 317         final MsaCompactor mc = new MsaCompactor( msa );
 318         mc.removeWorstOffenders( worst_offenders_to_remove, 1, realign, norm );
 319         return mc;
 320     }
 321 }