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