938273dd57f0fc086fe3b6499ec7af79418766e2
[jalview.git] / forester / java / src / org / forester / msa_compactor / MsaCompactor.java
1 // $Id:
2 // FORESTER -- software libraries and applications
3 // for evolutionary biology research and applications.
4 //
5 // Copyright (C) 2014 Christian M. Zmasek
6 // Copyright (C) 2014 Sanford-Burnham Medical Research Institute
7 // All rights reserved
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU Lesser General Public
11 // License as published by the Free Software Foundation; either
12 // version 2.1 of the License, or (at your option) any later version.
13 //
14 // This library is distributed in the hope that it will be useful,
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // Lesser General Public License for more details.
18 //
19 // You should have received a copy of the GNU Lesser General Public
20 // License along with this library; if not, write to the Free Software
21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
22 //
23 // WWW: https://sites.google.com/site/cmzmasek/home/software/forester
24
25 package org.forester.msa_compactor;
26
27 import java.io.File;
28 import java.io.IOException;
29 import java.io.Writer;
30 import java.math.RoundingMode;
31 import java.text.DecimalFormat;
32 import java.text.NumberFormat;
33 import java.util.ArrayList;
34 import java.util.Arrays;
35 import java.util.List;
36 import java.util.SortedSet;
37 import java.util.TreeSet;
38
39 import org.forester.evoinference.distance.NeighborJoiningF;
40 import org.forester.evoinference.distance.PairwiseDistanceCalculator;
41 import org.forester.evoinference.distance.PairwiseDistanceCalculator.PWD_DISTANCE_METHOD;
42 import org.forester.evoinference.matrix.distance.BasicSymmetricalDistanceMatrix;
43 import org.forester.evoinference.tools.BootstrapResampler;
44 import org.forester.io.writers.SequenceWriter;
45 import org.forester.io.writers.SequenceWriter.SEQ_FORMAT;
46 import org.forester.msa.DeleteableMsa;
47 import org.forester.msa.Mafft;
48 import org.forester.msa.Msa;
49 import org.forester.msa.Msa.MSA_FORMAT;
50 import org.forester.msa.MsaInferrer;
51 import org.forester.msa.MsaMethods;
52 import org.forester.msa.ResampleableMsa;
53 import org.forester.phylogeny.Phylogeny;
54 import org.forester.phylogeny.PhylogenyMethods;
55 import org.forester.sequence.Sequence;
56 import org.forester.tools.ConfidenceAssessor;
57 import org.forester.util.ForesterUtil;
58
59 public class MsaCompactor {
60
61     final private static NumberFormat NF_3                   = new DecimalFormat( "#.###" );
62     final private static NumberFormat NF_4                   = new DecimalFormat( "#.####" );
63     private double                    _gap_ratio             = -1;
64     //
65     private String                    _maffts_opts           = "--auto";
66     private int                       _min_length            = -1;
67     //
68     private DeleteableMsa             _msa                   = null;
69     private boolean                   _norm                  = true;
70     private File                      _out_file_base         = null;
71     private MSA_FORMAT                _output_format         = MSA_FORMAT.FASTA;
72     private String                    _path_to_mafft         = null;
73     //
74     private boolean                   _realign               = false;
75     private final SortedSet<String>   _removed_seq_ids;
76     private File                      _removed_seqs_out_base = null;
77
78     public final void setRemovedSeqsOutBase( final File removed_seqs_out_base ) {
79         _removed_seqs_out_base = removed_seqs_out_base;
80     }
81     private boolean                   _report_aln_mean_identity = false;
82     private int                       _step                     = -1;
83     private int                       _step_for_diagnostics     = -1;
84     private final short               _longest_id_length;
85     private final ArrayList<Sequence> _removed_seqs;
86     static {
87         NF_4.setRoundingMode( RoundingMode.HALF_UP );
88         NF_3.setRoundingMode( RoundingMode.HALF_UP );
89     }
90
91     public MsaCompactor( final DeleteableMsa msa ) {
92         _msa = msa;
93         _removed_seq_ids = new TreeSet<String>();
94         _longest_id_length = _msa.determineMaxIdLength();
95         _removed_seqs = new ArrayList<Sequence>();
96     }
97
98     final public Msa getMsa() {
99         return _msa;
100     }
101
102     final public SortedSet<String> getRemovedSeqIds() {
103         return _removed_seq_ids;
104     }
105
106     public final List<MsaProperties> removeViaGapAverage( final double mean_gapiness ) throws IOException,
107             InterruptedException {
108         final GapContribution stats[] = calcGapContribtionsStats( _norm );
109         final List<String> to_remove_ids = new ArrayList<String>();
110         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
111         for( final GapContribution gap_gontribution : stats ) {
112             to_remove_ids.add( gap_gontribution.getId() );
113         }
114         printTableHeader();
115         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
116         msa_props.add( msa_prop );
117         printMsaProperties( "", msa_prop );
118         System.out.println();
119         int i = 0;
120         while ( MsaMethods.calcGapRatio( _msa ) > mean_gapiness ) {
121             final String id = to_remove_ids.get( i );
122             _removed_seq_ids.add( id );
123             final Sequence deleted = _msa.deleteRow( id, true );
124             _removed_seqs.add( deleted );
125             removeGapColumns();
126             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( MsaMethods.calcGapRatio( _msa ) <= mean_gapiness ) ) {
127                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
128                 msa_props.add( msa_prop );
129                 System.out.println();
130             }
131             else if ( isPrintMsaStats( i ) ) {
132                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
133                 msa_props.add( msa_prop );
134                 printMsaProperties( id, msa_prop );
135                 System.out.println();
136             }
137             ++i;
138         }
139         if ( _removed_seqs_out_base != null ) {
140             final String msg = writeAndAlignRemovedSeqs();
141             System.out.println();
142             System.out.println( msg );
143         }
144         return msa_props;
145     }
146
147     public List<MsaProperties> removeViaLength( final int length ) throws IOException, InterruptedException {
148         final GapContribution stats[] = calcGapContribtionsStats( _norm );
149         final List<String> to_remove_ids = new ArrayList<String>();
150         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
151         for( final GapContribution gap_gontribution : stats ) {
152             to_remove_ids.add( gap_gontribution.getId() );
153         }
154         printTableHeader();
155         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
156         msa_props.add( msa_prop );
157         printMsaProperties( "", msa_prop );
158         System.out.println();
159         int i = 0;
160         while ( _msa.getLength() > length ) {
161             final String id = to_remove_ids.get( i );
162             _removed_seq_ids.add( id );
163             final Sequence deleted = _msa.deleteRow( id, true );
164             _removed_seqs.add( deleted );
165             removeGapColumns();
166             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( _msa.getLength() <= length ) ) {
167                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
168                 msa_props.add( msa_prop );
169                 System.out.println();
170             }
171             else if ( isPrintMsaStats( i ) ) {
172                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
173                 printMsaProperties( id, msa_prop );
174                 msa_props.add( msa_prop );
175                 System.out.println();
176             }
177             ++i;
178         }
179         if ( _removed_seqs_out_base != null ) {
180             final String msg = writeAndAlignRemovedSeqs();
181             System.out.println();
182             System.out.println( msg );
183         }
184         return msa_props;
185     }
186
187     public final List<MsaProperties> removeWorstOffenders( final int to_remove ) throws IOException,
188             InterruptedException {
189         final GapContribution stats[] = calcGapContribtionsStats( _norm );
190         final List<String> to_remove_ids = new ArrayList<String>();
191         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
192         for( int j = 0; j < to_remove; ++j ) {
193             to_remove_ids.add( stats[ j ].getId() );
194             _removed_seq_ids.add( stats[ j ].getId() );
195         }
196         printTableHeader();
197         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
198         msa_props.add( msa_prop );
199         printMsaProperties( "", msa_prop );
200         System.out.println();
201         for( int i = 0; i < to_remove_ids.size(); ++i ) {
202             final String id = to_remove_ids.get( i );
203             _removed_seq_ids.add( id );
204             final Sequence deleted = _msa.deleteRow( id, true );
205             _removed_seqs.add( deleted );
206             removeGapColumns();
207             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( i == ( to_remove_ids.size() - 1 ) ) ) {
208                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
209                 msa_props.add( msa_prop );
210                 System.out.println();
211             }
212             else if ( isPrintMsaStats( i ) ) {
213                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
214                 msa_props.add( msa_prop );
215                 printMsaProperties( id, msa_prop );
216                 System.out.println();
217             }
218         }
219         if ( _removed_seqs_out_base != null ) {
220             final String msg = writeAndAlignRemovedSeqs();
221             System.out.println();
222             System.out.println( msg );
223         }
224         return msa_props;
225     }
226
227     public final List<MsaProperties> chart( final int step, final boolean realign, final boolean norm )
228             throws IOException, InterruptedException {
229         final GapContribution stats[] = calcGapContribtionsStats( norm );
230         final List<String> to_remove_ids = new ArrayList<String>();
231         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
232         for( final GapContribution gap_gontribution : stats ) {
233             to_remove_ids.add( gap_gontribution.getId() );
234         }
235         printTableHeader();
236         final int x = ForesterUtil.roundToInt( _msa.getNumberOfSequences() / 20.0 );
237         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
238         msa_props.add( msa_prop );
239         printMsaProperties( "", msa_prop );
240         System.out.println();
241         int i = 0;
242         while ( _msa.getNumberOfSequences() > x ) {
243             final String id = to_remove_ids.get( i );
244             _msa.deleteRow( id, false );
245             if ( realign && isPrintMsaStatsWriteOutfileAndRealign( i ) ) {
246                 removeGapColumns();
247                 realignWithMafft();
248                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
249                 msa_props.add( msa_prop );
250                 printMsaProperties( id, msa_prop );
251                 System.out.print( "(realigned)" );
252                 System.out.println();
253             }
254             else if ( isPrintMsaStats( i ) ) {
255                 removeGapColumns();
256                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
257                 msa_props.add( msa_prop );
258                 printMsaProperties( id, msa_prop );
259                 System.out.println();
260             }
261             ++i;
262         }
263         return msa_props;
264     }
265
266     private final boolean isPrintMsaStats( final int i ) {
267         return ( ( ( _step < 2 ) && ( _step_for_diagnostics < 2 ) ) || ( ( _step_for_diagnostics > 0 ) && ( ( ( i + 1 ) % _step_for_diagnostics ) == 0 ) ) );
268     }
269
270     private final boolean isPrintMsaStatsWriteOutfileAndRealign( final int i ) {
271         return ( ( ( _step < 2 ) && ( _step_for_diagnostics < 2 ) ) || ( ( _step > 0 ) && ( ( ( i + 1 ) % _step ) == 0 ) ) );
272     }
273
274     public final void setGapRatio( final double gap_ratio ) {
275         _gap_ratio = gap_ratio;
276     }
277
278     public final void setMinLength( final int min_length ) {
279         _min_length = min_length;
280     }
281
282     public final void setNorm( final boolean norm ) {
283         _norm = norm;
284     }
285
286     final public void setOutFileBase( final File out_file_base ) {
287         _out_file_base = out_file_base;
288     }
289
290     public final void setOutputFormat( final MSA_FORMAT output_format ) {
291         _output_format = output_format;
292     }
293
294     public void setPathToMafft( final String path_to_mafft ) {
295         _path_to_mafft = path_to_mafft;
296     }
297
298     public final void setRealign( final boolean realign ) {
299         _realign = realign;
300     }
301
302     public final void setStep( final int step ) {
303         _step = step;
304     }
305
306     public final void setStepForDiagnostics( final int step_for_diagnostics ) {
307         _step_for_diagnostics = step_for_diagnostics;
308     }
309
310     public final void setReportAlnMeanIdentity( final boolean report_aln_mean_identity ) {
311         _report_aln_mean_identity = report_aln_mean_identity;
312     }
313
314     final public String writeMsa( final File outfile ) throws IOException {
315         final Double gr = MsaMethods.calcGapRatio( _msa );
316         final String s = outfile + "_" + _msa.getNumberOfSequences() + "_" + _msa.getLength() + "_"
317                 + ForesterUtil.roundToInt( gr * 100 );
318         writeMsa( _msa, s + obtainSuffix(), _output_format );
319         return s;
320     }
321
322     final public String writeAndAlignRemovedSeqs() throws IOException, InterruptedException {
323         final StringBuilder msg = new StringBuilder();
324         final String n = _removed_seqs_out_base + "_" + _removed_seqs.size() + ".fasta";
325         SequenceWriter.writeSeqs( _removed_seqs, new File( n ), SEQ_FORMAT.FASTA, 100 );
326         msg.append( "wrote " + _removed_seqs.size() + " removed sequences to " + "\"" + n + "\"" );
327         if ( _realign ) {
328             final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
329             final List<String> opts = new ArrayList<String>();
330             for( final String o : _maffts_opts.split( "\\s" ) ) {
331                 opts.add( o );
332             }
333             final Msa removed_msa = mafft.infer( _removed_seqs, opts );
334             final Double gr = MsaMethods.calcGapRatio( removed_msa );
335             String s = _removed_seqs_out_base + "_" + removed_msa.getNumberOfSequences() + "_"
336                     + removed_msa.getLength() + "_" + ForesterUtil.roundToInt( gr * 100 );
337             final String suffix = obtainSuffix();
338             s += suffix;
339             writeMsa( removed_msa, s, _output_format );
340             msg.append( ", and as MSA of length " + removed_msa.getLength() + " to \"" + s + "\"" );
341         }
342         return msg.toString();
343     }
344
345     private String obtainSuffix() {
346         if ( _output_format == MSA_FORMAT.FASTA ) {
347             return ".fasta";
348         }
349         else if ( _output_format == MSA_FORMAT.PHYLIP ) {
350             return ".aln";
351         }
352         return "";
353     }
354
355     final int calcNonGapResidues( final Sequence seq ) {
356         int ng = 0;
357         for( int i = 0; i < seq.getLength(); ++i ) {
358             if ( !seq.isGapAt( i ) ) {
359                 ++ng;
360             }
361         }
362         return ng;
363     }
364
365     private final GapContribution[] calcGapContribtions( final boolean normalize_for_effective_seq_length ) {
366         final double gappiness[] = calcGappiness();
367         final GapContribution stats[] = new GapContribution[ _msa.getNumberOfSequences() ];
368         for( int row = 0; row < _msa.getNumberOfSequences(); ++row ) {
369             stats[ row ] = new GapContribution( _msa.getIdentifier( row ) );
370             for( int col = 0; col < _msa.getLength(); ++col ) {
371                 if ( !_msa.isGapAt( row, col ) ) {
372                     stats[ row ].addToValue( gappiness[ col ] );
373                 }
374             }
375             if ( normalize_for_effective_seq_length ) {
376                 stats[ row ].divideValue( calcNonGapResidues( _msa.getSequence( row ) ) );
377             }
378             else {
379                 stats[ row ].divideValue( _msa.getLength() );
380             }
381         }
382         return stats;
383     }
384
385     final private GapContribution[] calcGapContribtionsStats( final boolean norm ) {
386         final GapContribution stats[] = calcGapContribtions( norm );
387         Arrays.sort( stats );
388         return stats;
389     }
390
391     private final double[] calcGappiness() {
392         final int l = _msa.getLength();
393         final double gappiness[] = new double[ l ];
394         final int seqs = _msa.getNumberOfSequences();
395         for( int i = 0; i < l; ++i ) {
396             gappiness[ i ] = ( double ) MsaMethods.calcGapSumPerColumn( _msa, i ) / seqs;
397         }
398         return gappiness;
399     }
400
401     private final Phylogeny inferNJphylogeny( final PWD_DISTANCE_METHOD pwd_distance_method,
402                                               final Msa msa,
403                                               final boolean write_matrix,
404                                               final String matrix_name ) {
405         BasicSymmetricalDistanceMatrix m = null;
406         switch ( pwd_distance_method ) {
407             case KIMURA_DISTANCE:
408                 m = PairwiseDistanceCalculator.calcKimuraDistances( msa );
409                 break;
410             case POISSON_DISTANCE:
411                 m = PairwiseDistanceCalculator.calcPoissonDistances( msa );
412                 break;
413             case FRACTIONAL_DISSIMILARITY:
414                 m = PairwiseDistanceCalculator.calcFractionalDissimilarities( msa );
415                 break;
416             default:
417                 throw new IllegalArgumentException( "invalid pwd method" );
418         }
419         if ( write_matrix ) {
420             try {
421                 m.write( ForesterUtil.createBufferedWriter( matrix_name ) );
422             }
423             catch ( final IOException e ) {
424                 e.printStackTrace();
425             }
426         }
427         final NeighborJoiningF nj = NeighborJoiningF.createInstance( false, 5 );
428         final Phylogeny phy = nj.execute( m );
429         return phy;
430     }
431
432     private final Phylogeny pi( final String matrix ) {
433         final Phylogeny master_phy = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, _msa, true, matrix );
434         final int seed = 15;
435         final int n = 100;
436         final ResampleableMsa resampleable_msa = new ResampleableMsa( _msa );
437         final int[][] resampled_column_positions = BootstrapResampler.createResampledColumnPositions( _msa.getLength(),
438                                                                                                       n,
439                                                                                                       seed );
440         final Phylogeny[] eval_phys = new Phylogeny[ n ];
441         for( int i = 0; i < n; ++i ) {
442             resampleable_msa.resample( resampled_column_positions[ i ] );
443             eval_phys[ i ] = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, resampleable_msa, false, null );
444         }
445         ConfidenceAssessor.evaluate( "bootstrap", eval_phys, master_phy, true, 1 );
446         PhylogenyMethods.extractFastaInformation( master_phy );
447         return master_phy;
448     }
449
450     private final void printMsaProperties( final String id, final MsaProperties msa_properties ) {
451         if ( ( _step < 2 ) || ( _step_for_diagnostics < 2 ) ) {
452             System.out.print( ForesterUtil.pad( id, _longest_id_length, ' ', false ) );
453             System.out.print( "\t" );
454         }
455         System.out.print( msaPropertiesAsSB( msa_properties ) );
456         System.out.print( "\t" );
457     }
458
459     private final StringBuilder msaPropertiesAsSB( final MsaProperties msa_properties ) {
460         final StringBuilder sb = new StringBuilder();
461         sb.append( msa_properties.getNumberOfSequences() );
462         sb.append( "\t" );
463         sb.append( msa_properties.getLength() );
464         sb.append( "\t" );
465         sb.append( NF_4.format( msa_properties.getGapRatio() ) );
466         if ( _report_aln_mean_identity /*msa_properties.getAverageIdentityRatio() >= 0*/) {
467             sb.append( "\t" );
468             sb.append( NF_4.format( msa_properties.getAverageIdentityRatio() ) );
469         }
470         return sb;
471     }
472
473     final private MsaProperties printMsaStatsWriteOutfileAndRealign( final boolean realign, final String id )
474             throws IOException, InterruptedException {
475         if ( realign ) {
476             realignWithMafft();
477         }
478         final MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
479         printMsaProperties( id, msa_prop );
480         final String s = writeOutfile();
481         System.out.print( "-> " + s + ( realign ? "\t(realigned)" : "" ) );
482         return msa_prop;
483     }
484
485     final private void realignWithMafft() throws IOException, InterruptedException {
486         final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
487         final List<String> opts = new ArrayList<String>();
488         for( final String o : _maffts_opts.split( "\\s" ) ) {
489             opts.add( o );
490         }
491         _msa = DeleteableMsa.createInstance( mafft.infer( _msa.asSequenceList(), opts ) );
492     }
493
494     public final void setMafftOptions( final String maffts_opts ) {
495         _maffts_opts = maffts_opts;
496     }
497
498     final private void removeGapColumns() {
499         _msa.deleteGapOnlyColumns();
500     }
501
502     final private static void writeMsa( final Msa msa, final String outfile, final MSA_FORMAT format )
503             throws IOException {
504         final Writer w = ForesterUtil.createBufferedWriter( outfile );
505         msa.write( w, format );
506         w.close();
507     }
508
509     private final String writeOutfile() throws IOException {
510         final String s = writeMsa( _out_file_base );
511         return s;
512     }
513
514     // Returns null if not path found.
515     final public static String guessPathToMafft() {
516         String path;
517         if ( ForesterUtil.OS_NAME.toLowerCase().indexOf( "win" ) >= 0 ) {
518             path = "C:\\Program Files\\mafft-win\\mafft.bat";
519             if ( MsaInferrer.isInstalled( path ) ) {
520                 return path;
521             }
522         }
523         path = "/home/czmasek/SOFTWARE/MSA/MAFFT/mafft-7.130-without-extensions/scripts/mafft";
524         if ( MsaInferrer.isInstalled( path ) ) {
525             return path;
526         }
527         path = "/usr/local/bin/mafft";
528         if ( MsaInferrer.isInstalled( path ) ) {
529             return path;
530         }
531         path = "/usr/bin/mafft";
532         if ( MsaInferrer.isInstalled( path ) ) {
533             return path;
534         }
535         path = "/bin/mafft";
536         if ( MsaInferrer.isInstalled( path ) ) {
537             return path;
538         }
539         path = "mafft";
540         if ( MsaInferrer.isInstalled( path ) ) {
541             return path;
542         }
543         return null;
544     }
545
546     private final void printTableHeader() {
547         if ( ( _step < 2 ) || ( _step_for_diagnostics < 2 ) ) {
548             System.out.print( ForesterUtil.pad( "Id", _longest_id_length, ' ', false ) );
549             System.out.print( "\t" );
550         }
551         System.out.print( "Seqs" );
552         System.out.print( "\t" );
553         System.out.print( "Length" );
554         System.out.print( "\t" );
555         System.out.print( "Gaps" );
556         System.out.print( "\t" );
557         if ( _report_aln_mean_identity ) {
558             System.out.print( "MSA qual" );
559             System.out.print( "\t" );
560         }
561         System.out.println();
562     }
563 }