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