inprogress
[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     private final short               _longest_id_length;
65     //
66     private String                    _maffts_opts              = "--auto";
67     private int                       _min_length               = -1;
68     //
69     private DeleteableMsa             _msa                      = null;
70     private boolean                   _norm                     = true;
71     private File                      _out_file_base            = null;
72     private MSA_FORMAT                _output_format            = MSA_FORMAT.FASTA;
73     private String                    _path_to_mafft            = null;
74     //
75     private boolean                   _realign                  = false;
76     private final SortedSet<String>   _removed_seq_ids;
77     private final ArrayList<Sequence> _removed_seqs;
78     private File                      _removed_seqs_out_base    = null;
79     private boolean                   _report_aln_mean_identity = false;
80     private int                       _step                     = -1;
81     private int                       _step_for_diagnostics     = -1;
82     static {
83         NF_4.setRoundingMode( RoundingMode.HALF_UP );
84         NF_3.setRoundingMode( RoundingMode.HALF_UP );
85     }
86
87     public MsaCompactor( final DeleteableMsa msa ) {
88         _msa = msa;
89         _removed_seq_ids = new TreeSet<String>();
90         _longest_id_length = _msa.determineMaxIdLength();
91         _removed_seqs = new ArrayList<Sequence>();
92     }
93
94     public final List<MsaProperties> chart( final int step, final boolean realign, final boolean norm )
95             throws IOException, InterruptedException {
96         final GapContribution stats[] = calcGapContribtionsStats( norm );
97         final List<String> to_remove_ids = new ArrayList<String>();
98         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
99         for( final GapContribution gap_gontribution : stats ) {
100             to_remove_ids.add( gap_gontribution.getId() );
101         }
102         printTableHeader();
103         final int x = ForesterUtil.roundToInt( _msa.getNumberOfSequences() / 20.0 );
104         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
105         msa_props.add( msa_prop );
106         printMsaProperties( "", msa_prop );
107         System.out.println();
108         int i = 0;
109         while ( _msa.getNumberOfSequences() > x ) {
110             final String id = to_remove_ids.get( i );
111             _msa.deleteRow( id, false );
112             if ( realign && isPrintMsaStatsWriteOutfileAndRealign( i ) ) {
113                 removeGapColumns();
114                 realignWithMafft();
115                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
116                 msa_props.add( msa_prop );
117                 printMsaProperties( id, msa_prop );
118                 System.out.print( "(realigned)" );
119                 System.out.println();
120             }
121             else if ( isPrintMsaStats( i ) ) {
122                 removeGapColumns();
123                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
124                 msa_props.add( msa_prop );
125                 printMsaProperties( id, msa_prop );
126                 System.out.println();
127             }
128             ++i;
129         }
130         return msa_props;
131     }
132
133     final public Msa getMsa() {
134         return _msa;
135     }
136
137     final public SortedSet<String> getRemovedSeqIds() {
138         return _removed_seq_ids;
139     }
140
141     public final List<MsaProperties> removeViaGapAverage( final double mean_gapiness ) throws IOException,
142             InterruptedException {
143         final GapContribution stats[] = calcGapContribtionsStats( _norm );
144         final List<String> to_remove_ids = new ArrayList<String>();
145         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
146         for( final GapContribution gap_gontribution : stats ) {
147             to_remove_ids.add( gap_gontribution.getId() );
148         }
149         printTableHeader();
150         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
151         msa_props.add( msa_prop );
152         printMsaProperties( "", msa_prop );
153         System.out.println();
154         int i = 0;
155         while ( MsaMethods.calcGapRatio( _msa ) > mean_gapiness ) {
156             final String id = to_remove_ids.get( i );
157             _removed_seq_ids.add( id );
158             final Sequence deleted = _msa.deleteRow( id, true );
159             _removed_seqs.add( deleted );
160             removeGapColumns();
161             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( MsaMethods.calcGapRatio( _msa ) <= mean_gapiness ) ) {
162                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
163                 msa_props.add( msa_prop );
164                 System.out.println();
165             }
166             else if ( isPrintMsaStats( i ) ) {
167                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
168                 msa_props.add( msa_prop );
169                 printMsaProperties( id, msa_prop );
170                 System.out.println();
171             }
172             ++i;
173         }
174         if ( _removed_seqs_out_base != null ) {
175             final String msg = writeAndAlignRemovedSeqs();
176             System.out.println();
177             System.out.println( msg );
178         }
179         return msa_props;
180     }
181
182     public List<MsaProperties> removeViaLength( final int length ) throws IOException, InterruptedException {
183         final GapContribution stats[] = calcGapContribtionsStats( _norm );
184         final List<String> to_remove_ids = new ArrayList<String>();
185         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
186         for( final GapContribution gap_gontribution : stats ) {
187             to_remove_ids.add( gap_gontribution.getId() );
188         }
189         printTableHeader();
190         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
191         msa_props.add( msa_prop );
192         printMsaProperties( "", msa_prop );
193         System.out.println();
194         int i = 0;
195         while ( _msa.getLength() > length ) {
196             final String id = to_remove_ids.get( i );
197             _removed_seq_ids.add( id );
198             final Sequence deleted = _msa.deleteRow( id, true );
199             _removed_seqs.add( deleted );
200             removeGapColumns();
201             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( _msa.getLength() <= length ) ) {
202                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
203                 msa_props.add( msa_prop );
204                 System.out.println();
205             }
206             else if ( isPrintMsaStats( i ) ) {
207                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
208                 printMsaProperties( id, msa_prop );
209                 msa_props.add( msa_prop );
210                 System.out.println();
211             }
212             ++i;
213         }
214         if ( _removed_seqs_out_base != null ) {
215             final String msg = writeAndAlignRemovedSeqs();
216             System.out.println();
217             System.out.println( msg );
218         }
219         return msa_props;
220     }
221
222     public final List<MsaProperties> removeWorstOffenders( final int to_remove ) throws IOException,
223             InterruptedException {
224         final GapContribution stats[] = calcGapContribtionsStats( _norm );
225         final List<String> to_remove_ids = new ArrayList<String>();
226         final List<MsaProperties> msa_props = new ArrayList<MsaProperties>();
227         for( int j = 0; j < to_remove; ++j ) {
228             to_remove_ids.add( stats[ j ].getId() );
229             _removed_seq_ids.add( stats[ j ].getId() );
230         }
231         printTableHeader();
232         MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
233         msa_props.add( msa_prop );
234         printMsaProperties( "", msa_prop );
235         System.out.println();
236         for( int i = 0; i < to_remove_ids.size(); ++i ) {
237             final String id = to_remove_ids.get( i );
238             _removed_seq_ids.add( id );
239             final Sequence deleted = _msa.deleteRow( id, true );
240             _removed_seqs.add( deleted );
241             removeGapColumns();
242             if ( isPrintMsaStatsWriteOutfileAndRealign( i ) || ( i == ( to_remove_ids.size() - 1 ) ) ) {
243                 msa_prop = printMsaStatsWriteOutfileAndRealign( _realign, id );
244                 msa_props.add( msa_prop );
245                 System.out.println();
246             }
247             else if ( isPrintMsaStats( i ) ) {
248                 msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
249                 msa_props.add( msa_prop );
250                 printMsaProperties( id, msa_prop );
251                 System.out.println();
252             }
253         }
254         if ( _removed_seqs_out_base != null ) {
255             final String msg = writeAndAlignRemovedSeqs();
256             System.out.println();
257             System.out.println( msg );
258         }
259         return msa_props;
260     }
261
262     public final void setGapRatio( final double gap_ratio ) {
263         _gap_ratio = gap_ratio;
264     }
265
266     public final void setMafftOptions( final String maffts_opts ) {
267         _maffts_opts = maffts_opts;
268     }
269
270     public final void setMinLength( final int min_length ) {
271         _min_length = min_length;
272     }
273
274     public final void setNorm( final boolean norm ) {
275         _norm = norm;
276     }
277
278     final public void setOutFileBase( final File out_file_base ) {
279         _out_file_base = out_file_base;
280     }
281
282     public final void setOutputFormat( final MSA_FORMAT output_format ) {
283         _output_format = output_format;
284     }
285
286     public void setPathToMafft( final String path_to_mafft ) {
287         _path_to_mafft = path_to_mafft;
288     }
289
290     public final void setRealign( final boolean realign ) {
291         _realign = realign;
292     }
293
294     public final void setRemovedSeqsOutBase( final File removed_seqs_out_base ) {
295         _removed_seqs_out_base = removed_seqs_out_base;
296     }
297
298     public final void setReportAlnMeanIdentity( final boolean report_aln_mean_identity ) {
299         _report_aln_mean_identity = report_aln_mean_identity;
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     final public String writeAndAlignRemovedSeqs() throws IOException, InterruptedException {
311         final StringBuilder msg = new StringBuilder();
312         final String n = _removed_seqs_out_base + "_" + _removed_seqs.size() + ".fasta";
313         SequenceWriter.writeSeqs( _removed_seqs, new File( n ), SEQ_FORMAT.FASTA, 100 );
314         msg.append( "wrote " + _removed_seqs.size() + " removed sequences to " + "\"" + n + "\"" );
315         if ( _realign ) {
316             final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
317             final List<String> opts = new ArrayList<String>();
318             for( final String o : _maffts_opts.split( "\\s" ) ) {
319                 opts.add( o );
320             }
321             final Msa removed_msa = mafft.infer( _removed_seqs, opts );
322             final Double gr = MsaMethods.calcGapRatio( removed_msa );
323             String s = _removed_seqs_out_base + "_" + removed_msa.getNumberOfSequences() + "_"
324                     + removed_msa.getLength() + "_" + ForesterUtil.roundToInt( gr * 100 );
325             final String suffix = obtainSuffix();
326             s += suffix;
327             writeMsa( removed_msa, s, _output_format );
328             msg.append( ", and as MSA of length " + removed_msa.getLength() + " to \"" + s + "\"" );
329         }
330         return msg.toString();
331     }
332
333     final public String writeMsa( final File outfile ) throws IOException {
334         final Double gr = MsaMethods.calcGapRatio( _msa );
335         final String s = outfile + "_" + _msa.getNumberOfSequences() + "_" + _msa.getLength() + "_"
336                 + ForesterUtil.roundToInt( gr * 100 );
337         writeMsa( _msa, s + obtainSuffix(), _output_format );
338         return s;
339     }
340
341     final int calcNonGapResidues( final Sequence seq ) {
342         int ng = 0;
343         for( int i = 0; i < seq.getLength(); ++i ) {
344             if ( !seq.isGapAt( i ) ) {
345                 ++ng;
346             }
347         }
348         return ng;
349     }
350
351     private final GapContribution[] calcGapContribtions( final boolean normalize_for_effective_seq_length ) {
352         final double gappiness[] = calcGappiness();
353         final GapContribution stats[] = new GapContribution[ _msa.getNumberOfSequences() ];
354         for( int row = 0; row < _msa.getNumberOfSequences(); ++row ) {
355             stats[ row ] = new GapContribution( _msa.getIdentifier( row ) );
356             for( int col = 0; col < _msa.getLength(); ++col ) {
357                 if ( !_msa.isGapAt( row, col ) ) {
358                     stats[ row ].addToValue( gappiness[ col ] );
359                 }
360             }
361             if ( normalize_for_effective_seq_length ) {
362                 stats[ row ].divideValue( calcNonGapResidues( _msa.getSequence( row ) ) );
363             }
364             else {
365                 stats[ row ].divideValue( _msa.getLength() );
366             }
367         }
368         return stats;
369     }
370
371     final private GapContribution[] calcGapContribtionsStats( final boolean norm ) {
372         final GapContribution stats[] = calcGapContribtions( norm );
373         Arrays.sort( stats );
374         return stats;
375     }
376
377     private final double[] calcGappiness() {
378         final int l = _msa.getLength();
379         final double gappiness[] = new double[ l ];
380         final int seqs = _msa.getNumberOfSequences();
381         for( int i = 0; i < l; ++i ) {
382             gappiness[ i ] = ( double ) MsaMethods.calcGapSumPerColumn( _msa, i ) / seqs;
383         }
384         return gappiness;
385     }
386
387     private final Phylogeny inferNJphylogeny( final PWD_DISTANCE_METHOD pwd_distance_method,
388                                               final Msa msa,
389                                               final boolean write_matrix,
390                                               final String matrix_name ) {
391         BasicSymmetricalDistanceMatrix m = null;
392         switch ( pwd_distance_method ) {
393             case KIMURA_DISTANCE:
394                 m = PairwiseDistanceCalculator.calcKimuraDistances( msa );
395                 break;
396             case POISSON_DISTANCE:
397                 m = PairwiseDistanceCalculator.calcPoissonDistances( msa );
398                 break;
399             case FRACTIONAL_DISSIMILARITY:
400                 m = PairwiseDistanceCalculator.calcFractionalDissimilarities( msa );
401                 break;
402             default:
403                 throw new IllegalArgumentException( "invalid pwd method" );
404         }
405         if ( write_matrix ) {
406             try {
407                 m.write( ForesterUtil.createBufferedWriter( matrix_name ) );
408             }
409             catch ( final IOException e ) {
410                 e.printStackTrace();
411             }
412         }
413         final NeighborJoiningF nj = NeighborJoiningF.createInstance( false, 5 );
414         final Phylogeny phy = nj.execute( m );
415         return phy;
416     }
417
418     private final boolean isPrintMsaStats( final int i ) {
419         return ( ( ( _step < 2 ) && ( _step_for_diagnostics < 2 ) ) || ( ( _step_for_diagnostics > 0 ) && ( ( ( i + 1 ) % _step_for_diagnostics ) == 0 ) ) );
420     }
421
422     private final boolean isPrintMsaStatsWriteOutfileAndRealign( final int i ) {
423         return ( ( ( _step < 2 ) && ( _step_for_diagnostics < 2 ) ) || ( ( _step > 0 ) && ( ( ( i + 1 ) % _step ) == 0 ) ) );
424     }
425
426     private final StringBuilder msaPropertiesAsSB( final MsaProperties msa_properties ) {
427         final StringBuilder sb = new StringBuilder();
428         sb.append( msa_properties.getNumberOfSequences() );
429         sb.append( "\t" );
430         sb.append( msa_properties.getLength() );
431         sb.append( "\t" );
432         sb.append( NF_4.format( msa_properties.getGapRatio() ) );
433         if ( _report_aln_mean_identity /*msa_properties.getAverageIdentityRatio() >= 0*/) {
434             sb.append( "\t" );
435             sb.append( NF_4.format( msa_properties.getAverageIdentityRatio() ) );
436         }
437         return sb;
438     }
439
440     private String obtainSuffix() {
441         if ( _output_format == MSA_FORMAT.FASTA ) {
442             return ".fasta";
443         }
444         else if ( _output_format == MSA_FORMAT.PHYLIP ) {
445             return ".aln";
446         }
447         return "";
448     }
449
450     private final Phylogeny pi( final String matrix ) {
451         final Phylogeny master_phy = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, _msa, true, matrix );
452         final int seed = 15;
453         final int n = 100;
454         final ResampleableMsa resampleable_msa = new ResampleableMsa( _msa );
455         final int[][] resampled_column_positions = BootstrapResampler.createResampledColumnPositions( _msa.getLength(),
456                                                                                                       n,
457                                                                                                       seed );
458         final Phylogeny[] eval_phys = new Phylogeny[ n ];
459         for( int i = 0; i < n; ++i ) {
460             resampleable_msa.resample( resampled_column_positions[ i ] );
461             eval_phys[ i ] = inferNJphylogeny( PWD_DISTANCE_METHOD.KIMURA_DISTANCE, resampleable_msa, false, null );
462         }
463         ConfidenceAssessor.evaluate( "bootstrap", eval_phys, master_phy, true, 1 );
464         PhylogenyMethods.extractFastaInformation( master_phy );
465         return master_phy;
466     }
467
468     private final void printMsaProperties( final String id, final MsaProperties msa_properties ) {
469         if ( ( _step < 2 ) || ( _step_for_diagnostics < 2 ) ) {
470             System.out.print( ForesterUtil.pad( id, _longest_id_length, ' ', false ) );
471             System.out.print( "\t" );
472         }
473         System.out.print( msaPropertiesAsSB( msa_properties ) );
474         System.out.print( "\t" );
475     }
476
477     final private MsaProperties printMsaStatsWriteOutfileAndRealign( final boolean realign, final String id )
478             throws IOException, InterruptedException {
479         if ( realign ) {
480             realignWithMafft();
481         }
482         final MsaProperties msa_prop = new MsaProperties( _msa, _report_aln_mean_identity );
483         printMsaProperties( id, msa_prop );
484         final String s = writeOutfile();
485         System.out.print( "-> " + s + ( realign ? "\t(realigned)" : "" ) );
486         return msa_prop;
487     }
488
489     private final void printTableHeader() {
490         if ( ( _step < 2 ) || ( _step_for_diagnostics < 2 ) ) {
491             System.out.print( ForesterUtil.pad( "Id", _longest_id_length, ' ', false ) );
492             System.out.print( "\t" );
493         }
494         System.out.print( "Seqs" );
495         System.out.print( "\t" );
496         System.out.print( "Length" );
497         System.out.print( "\t" );
498         System.out.print( "Gaps" );
499         System.out.print( "\t" );
500         if ( _report_aln_mean_identity ) {
501             System.out.print( "MSA qual" );
502             System.out.print( "\t" );
503         }
504         System.out.println();
505     }
506
507     final private void realignWithMafft() throws IOException, InterruptedException {
508         final MsaInferrer mafft = Mafft.createInstance( _path_to_mafft );
509         final List<String> opts = new ArrayList<String>();
510         for( final String o : _maffts_opts.split( "\\s" ) ) {
511             opts.add( o );
512         }
513         _msa = DeleteableMsa.createInstance( mafft.infer( _msa.asSequenceList(), opts ) );
514     }
515
516     final private void removeGapColumns() {
517         _msa.deleteGapOnlyColumns();
518     }
519
520     private final String writeOutfile() throws IOException {
521         final String s = writeMsa( _out_file_base );
522         return s;
523     }
524
525     // Returns null if not path found.
526     final public static String guessPathToMafft() {
527         String path;
528         if ( ForesterUtil.OS_NAME.toLowerCase().indexOf( "win" ) >= 0 ) {
529             path = "C:\\Program Files\\mafft-win\\mafft.bat";
530             if ( MsaInferrer.isInstalled( path ) ) {
531                 return path;
532             }
533         }
534         path = "/home/czmasek/SOFTWARE/MSA/MAFFT/mafft-7.130-without-extensions/scripts/mafft";
535         if ( MsaInferrer.isInstalled( path ) ) {
536             return path;
537         }
538         path = "/usr/local/bin/mafft";
539         if ( MsaInferrer.isInstalled( path ) ) {
540             return path;
541         }
542         path = "/usr/bin/mafft";
543         if ( MsaInferrer.isInstalled( path ) ) {
544             return path;
545         }
546         path = "/bin/mafft";
547         if ( MsaInferrer.isInstalled( path ) ) {
548             return path;
549         }
550         path = "mafft";
551         if ( MsaInferrer.isInstalled( path ) ) {
552             return path;
553         }
554         return null;
555     }
556
557     final private static void writeMsa( final Msa msa, final String outfile, final MSA_FORMAT format )
558             throws IOException {
559         final Writer w = ForesterUtil.createBufferedWriter( outfile );
560         msa.write( w, format );
561         w.close();
562     }
563 }