forester/java/src/org/forester/archaeopteryx/tools/PhylogeneticInferrer.java

   1 // $Id:
   2 // forester -- software libraries and applications
   3 // for genomics and evolutionary biology research.
   4 //
   5 // Copyright (C) 2010 Christian M Zmasek
   6 // Copyright (C) 2010 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 // Contact: phylosoft @ gmail . com
  24 // WWW: www.phylosoft.org/forester
  25
  26 package org.forester.archaeopteryx.tools;
  27
  28 import java.io.BufferedWriter;
  29 import java.io.File;
  30 import java.io.FileWriter;
  31 import java.io.IOException;
  32 import java.util.ArrayList;
  33 import java.util.List;
  34 import java.util.regex.Matcher;
  35
  36 import javax.swing.JOptionPane;
  37
  38 import org.forester.archaeopteryx.AptxUtil;
  39 import org.forester.archaeopteryx.MainFrameApplication;
  40 import org.forester.evoinference.distance.NeighborJoining;
  41 import org.forester.evoinference.distance.PairwiseDistanceCalculator;
  42 import org.forester.evoinference.matrix.distance.BasicSymmetricalDistanceMatrix;
  43 import org.forester.evoinference.tools.BootstrapResampler;
  44 import org.forester.io.parsers.FastaParser;
  45 import org.forester.io.writers.SequenceWriter;
  46 import org.forester.io.writers.SequenceWriter.SEQ_FORMAT;
  47 import org.forester.msa.BasicMsa;
  48 import org.forester.msa.Mafft;
  49 import org.forester.msa.Msa;
  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.PhylogenyNode;
  55 import org.forester.phylogeny.data.Accession;
  56 import org.forester.phylogeny.iterators.PhylogenyNodeIterator;
  57 import org.forester.sequence.Sequence;
  58 import org.forester.tools.ConfidenceAssessor;
  59 import org.forester.util.ForesterUtil;
  60
  61 public class PhylogeneticInferrer implements Runnable {
  62
  63     private Msa                                _msa;
  64     private final MainFrameApplication         _mf;
  65     private final PhylogeneticInferenceOptions _options;
  66     private final List<Sequence>               _seqs;
  67     private final boolean DEBUG = true;
  68     public final static String                 MSA_FILE_SUFFIX = ".aln";
  69     public final static String                 PWD_FILE_SUFFIX = ".pwd";
  70
  71     public PhylogeneticInferrer( final List<Sequence> seqs,
  72                                  final PhylogeneticInferenceOptions options,
  73                                  final MainFrameApplication mf ) {
  74         _msa = null;
  75         _seqs = seqs;
  76         _mf = mf;
  77         _options = options;
  78     }
  79
  80     public PhylogeneticInferrer( final Msa msa,
  81                                  final PhylogeneticInferenceOptions options,
  82                                  final MainFrameApplication mf ) {
  83         _msa = msa;
  84         _seqs = null;
  85         _mf = mf;
  86         _options = options;
  87     }
  88
  89     private Msa inferMsa() throws IOException, InterruptedException {
  90         final File temp_seqs_file = File.createTempFile( "__msa__temp__", ".fasta" );
  91         if ( DEBUG ) {
  92         System.out.println();
  93         System.out.println( "temp file: " + temp_seqs_file );
  94         System.out.println();
  95         }
  96         //final File temp_seqs_file = new File( _options.getTempDir() + ForesterUtil.FILE_SEPARATOR + "s.fasta" );
  97         final BufferedWriter writer = new BufferedWriter( new FileWriter( temp_seqs_file ) );
  98         SequenceWriter.writeSeqs( _seqs, writer, SEQ_FORMAT.FASTA, 100 );
  99         writer.close();
 100         final List<String> opts = processMafftOptions();
 101         return runMAFFT( temp_seqs_file, opts );
 102     }
 103
 104     private List<String> processMafftOptions() {
 105         final String opts_str = _options.getMsaPrgParameters().trim().toLowerCase();
 106         final String[] opts_ary = opts_str.split( " " );
 107         final List<String> opts = new ArrayList<String>();
 108         boolean saw_quiet = false;
 109         for( final String opt : opts_ary ) {
 110             opts.add( opt );
 111             if ( opt.equals( "--quiet" ) ) {
 112                 saw_quiet = true;
 113             }
 114         }
 115         if ( !saw_quiet ) {
 116             opts.add( "--quiet" );
 117         }
 118         return opts;
 119     }
 120
 121     private Phylogeny inferPhylogeny( final Msa msa ) {
 122         BasicSymmetricalDistanceMatrix m = null;
 123         switch ( _options.getPwdDistanceMethod() ) {
 124             case KIMURA_DISTANCE:
 125                 m = PairwiseDistanceCalculator.calcKimuraDistances( msa );
 126                 break;
 127             case POISSON_DISTANCE:
 128                 m = PairwiseDistanceCalculator.calcPoissonDistances( msa );
 129                 break;
 130             case FRACTIONAL_DISSIMILARITY:
 131                 m = PairwiseDistanceCalculator.calcFractionalDissimilarities( msa );
 132                 break;
 133             default:
 134                 throw new RuntimeException( "invalid pwd method" );
 135         }
 136         if ( !ForesterUtil.isEmpty( _options.getIntermediateFilesBase() ) ) {
 137             BufferedWriter pwd_writer;
 138             try {
 139                 pwd_writer = new BufferedWriter( new FileWriter( _options.getIntermediateFilesBase() + PWD_FILE_SUFFIX ) );
 140                 m.write( pwd_writer );
 141                 pwd_writer.close();
 142             }
 143             catch ( final IOException e ) {
 144                 // TODO Auto-generated catch block
 145                 e.printStackTrace();
 146             }
 147         }
 148         final NeighborJoining nj = new NeighborJoining();
 149         final Phylogeny phy = nj.execute( m );
 150         PhylogeneticInferrer.extractFastaInformation( phy );
 151         return phy;
 152     }
 153
 154     private void infer() throws InterruptedException {
 155         //_mf.getMainPanel().getCurrentTreePanel().setWaitCursor();
 156         if ( ( _msa == null ) && ( _seqs == null ) ) {
 157             throw new IllegalArgumentException( "cannot run phylogenetic analysis with null msa and seq array" );
 158         }
 159         if ( _msa == null ) {
 160             Msa msa = null;
 161             try {
 162                 msa = inferMsa();
 163             }
 164             catch ( final IOException e ) {
 165                 JOptionPane.showMessageDialog( _mf,
 166                                                "Could not create multiple sequence alignment with \""
 167                                                        + _options.getMsaPrg() + "\" and the following parameters:\n\""
 168                                                        + _options.getMsaPrgParameters() + "\"\nError: "
 169                                                        + e.getLocalizedMessage(),
 170                                                "Failed to Calculate MSA",
 171                                                JOptionPane.ERROR_MESSAGE );
 172                 if ( DEBUG ) {
 173                 e.printStackTrace();
 174                 }
 175                 return;
 176             }
 177             catch ( final Exception e ) {
 178                 JOptionPane.showMessageDialog( _mf,
 179                                                "Could not create multiple sequence alignment with \""
 180                                                        + _options.getMsaPrg() + "\" and the following parameters:\n\""
 181                                                        + _options.getMsaPrgParameters() + "\"\nError: "
 182                                                        + e.getLocalizedMessage(),
 183                                                "Unexpected Exception During MSA Calculation",
 184                                                JOptionPane.ERROR_MESSAGE );
 185                 if ( DEBUG ) {
 186                 e.printStackTrace();
 187                 }
 188                 return;
 189             }
 190             if ( msa == null ) {
 191                 JOptionPane.showMessageDialog( _mf,
 192                                                "Could not create multiple sequence alignment with "
 193                                                        + _options.getMsaPrg() + "\nand the following parameters:\n\""
 194                                                        + _options.getMsaPrgParameters() + "\"",
 195                                                "Failed to Calculate MSA",
 196                                                JOptionPane.ERROR_MESSAGE );
 197                 return;
 198             }
 199             if ( DEBUG ) {
 200             System.out.println( msa.toString() );
 201             System.out.println( MsaMethods.calcBasicGapinessStatistics( msa ).toString() );
 202             }
 203             final MsaMethods msa_tools = MsaMethods.createInstance();
 204             if ( _options.isExecuteMsaProcessing() ) {
 205                 msa = msa_tools.removeGapColumns( _options.getMsaProcessingMaxAllowedGapRatio(),
 206                                                   _options.getMsaProcessingMinAllowedLength(),
 207                                                   msa );
 208                 if ( msa == null ) {
 209                     JOptionPane.showMessageDialog( _mf,
 210                                                    "Less than two sequences longer than "
 211                                                            + _options.getMsaProcessingMinAllowedLength()
 212                                                            + " residues left after MSA processing",
 213                                                    "MSA Processing Settings Too Stringent",
 214                                                    JOptionPane.ERROR_MESSAGE );
 215                     return;
 216                 }
 217             }
 218             if ( DEBUG ) {
 219             System.out.println( msa_tools.getIgnoredSequenceIds() );
 220             System.out.println( msa.toString() );
 221             System.out.println( MsaMethods.calcBasicGapinessStatistics( msa ).toString() );
 222             }
 223             _msa = msa;
 224         }
 225         final int n = _options.getBootstrapSamples();
 226         final long seed = _options.getRandomNumberGeneratorSeed();
 227         final Phylogeny master_phy = inferPhylogeny( _msa );
 228         if ( _options.isPerformBootstrapResampling() && ( n > 0 ) ) {
 229             final ResampleableMsa resampleable_msa = new ResampleableMsa( ( BasicMsa ) _msa );
 230             final int[][] resampled_column_positions = BootstrapResampler.createResampledColumnPositions( _msa
 231                     .getLength(), n, seed );
 232             final Phylogeny[] eval_phys = new Phylogeny[ n ];
 233             for( int i = 0; i < n; ++i ) {
 234                 resampleable_msa.resample( resampled_column_positions[ i ] );
 235                 eval_phys[ i ] = inferPhylogeny( resampleable_msa );
 236             }
 237             ConfidenceAssessor.evaluate( "bootstrap", eval_phys, master_phy, true, 1 );
 238         }
 239         _mf.getMainPanel().addPhylogenyInNewTab( master_phy, _mf.getConfiguration(), "nj", "njpath" );
 240         _mf.getMainPanel().getCurrentTreePanel().setArrowCursor();
 241         JOptionPane.showMessageDialog( _mf,
 242                                        "Inference successfully completed",
 243                                        "Inference Completed",
 244                                        JOptionPane.INFORMATION_MESSAGE );
 245     }
 246
 247     @Override
 248     public void run() {
 249         try {
 250             infer();
 251         }
 252         catch ( final InterruptedException e ) {
 253             // TODO need to handle this exception SOMEHOW!
 254             // TODO Auto-generated catch block
 255             e.printStackTrace();
 256         }
 257     }
 258
 259     private Msa runMAFFT( final File input_seqs, final List<String> opts ) throws IOException, InterruptedException {
 260         Msa msa = null;
 261         final MsaInferrer mafft = Mafft.createInstance();
 262         try {
 263             msa = mafft.infer( input_seqs, opts );
 264         }
 265         catch ( final IOException e ) {
 266             System.out.println( mafft.getErrorDescription() );
 267         }
 268         return msa;
 269     }
 270
 271     private void writeToFiles( final BasicSymmetricalDistanceMatrix m ) {
 272         if ( !ForesterUtil.isEmpty( _options.getIntermediateFilesBase() ) ) {
 273             try {
 274                 final BufferedWriter msa_writer = new BufferedWriter( new FileWriter( _options.getIntermediateFilesBase()
 275                         + MSA_FILE_SUFFIX ) );
 276                 _msa.write( msa_writer );
 277                 msa_writer.close();
 278                 final BufferedWriter pwd_writer = new BufferedWriter( new FileWriter( _options.getIntermediateFilesBase()
 279                         + PWD_FILE_SUFFIX ) );
 280                 m.write( pwd_writer );
 281                 pwd_writer.close();
 282             }
 283             catch ( final Exception e ) {
 284                 System.out.println( "Error: " + e.getMessage() );
 285             }
 286         }
 287     }
 288
 289     public static void extractFastaInformation( final Phylogeny phy ) {
 290         for( final PhylogenyNodeIterator iter = phy.iteratorExternalForward(); iter.hasNext(); ) {
 291             final PhylogenyNode node = iter.next();
 292             if ( !ForesterUtil.isEmpty( node.getName() ) ) {
 293                 final Matcher name_m = FastaParser.FASTA_DESC_LINE.matcher( node.getName() );
 294                 if ( name_m.lookingAt() ) {
 295                     System.out.println();
 296                     // System.out.println( name_m.group( 1 ) );
 297                     // System.out.println( name_m.group( 2 ) );
 298                     // System.out.println( name_m.group( 3 ) );
 299                     // System.out.println( name_m.group( 4 ) );
 300                     final String acc_source = name_m.group( 1 );
 301                     final String acc = name_m.group( 2 );
 302                     final String seq_name = name_m.group( 3 );
 303                     final String tax_sn = name_m.group( 4 );
 304                     if ( !ForesterUtil.isEmpty( acc_source ) && !ForesterUtil.isEmpty( acc ) ) {
 305                         AptxUtil.ensurePresenceOfSequence( node );
 306                         node.getNodeData().getSequence( 0 ).setAccession( new Accession( acc, acc_source ) );
 307                     }
 308                     if ( !ForesterUtil.isEmpty( seq_name ) ) {
 309                         AptxUtil.ensurePresenceOfSequence( node );
 310                         node.getNodeData().getSequence( 0 ).setName( seq_name );
 311                     }
 312                     if ( !ForesterUtil.isEmpty( tax_sn ) ) {
 313                         AptxUtil.ensurePresenceOfTaxonomy( node );
 314                         node.getNodeData().getTaxonomy( 0 ).setScientificName( tax_sn );
 315                     }
 316                 }
 317             }
 318         }
 319     }
 320 }