libMems/MuscleInterface.cpp

Go to the documentation of this file.

/*******************************************************************************
 * $Id: MuscleInterface.cpp,v 1.27 2004/03/01 02:40:08 darling Exp $
 * This file is copyright 2002-2007 Aaron Darling and authors listed in the AUTHORS file.
 * Please see the file called COPYING for licensing, copying, and modification
 * Please see the file called COPYING for licensing details.
 * **************
 ******************************************************************************/

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "libMems/MuscleInterface.h"

#include "libGenome/gnFilter.h"
#include "libGenome/gnFASSource.h"
#include "libGenome/gnStringTools.h"
#include "libMUSCLE/muscle.h"
#include "libMUSCLE/params.h"
#include "libMUSCLE/msa.h"
#include "libMUSCLE/seq.h"
#include "libMUSCLE/seqvect.h"
#include "libMUSCLE/tree.h"
#include "libMUSCLE/clust.h"
#include "libMUSCLE/profile.h"
#include "libMUSCLE/distfunc.h"
#include "libMUSCLE/clustsetdf.h"
#include "libMUSCLE/textfile.h"
#include "libMUSCLE/types.h"
#include "boost/algorithm/string/erase.hpp"
#include "boost/algorithm/string/case_conv.hpp"

#include <sstream>
#include <fstream>

// this gets defined in muscle.cpp, but not declared in any headers
extern void MUSCLE(SeqVect &v, MSA &msaOut);
extern void RefineW(const MSA &msaIn, MSA &msaOut);

using namespace std;
using namespace genome;
namespace mems {

bool debug_muscle = false;

bool pipeExec( char** cmd_argv, const string& command, const string& input, string& output, string& error );

char** parseCommand( const string& cmd );
char** parseCommand( const string& cmd ){
        // first count tokens

        // tokenize on "
        stringstream qs( cmd );
        string cur_str;
        boolean in_quote = true;
        int token_count = 0;
        vector< string > cmd_tokens;
        while( getline( qs, cur_str, '"' ) ){
                // never start out in a quote
                in_quote = !in_quote;
                if( cur_str.length() == 0 )
                        continue;
                if( in_quote ){
                        cmd_tokens.push_back( cur_str );
                }else{
                        stringstream ss( cur_str );
                        string asdf;
                        while( ss >> asdf )
                                cmd_tokens.push_back( asdf );
                }
        }
        char ** cmd_array = new char*[ cmd_tokens.size() + 1 ];
        for( int tokI = 0; tokI < cmd_tokens.size(); tokI++ ){
                cmd_array[ tokI ] = new char[ cmd_tokens[ tokI ].length() + 1 ];
                strcpy( cmd_array[ tokI ], cmd_tokens[ tokI ].c_str() );
        }
        cmd_array[ cmd_tokens.size() ] = NULL;
        return cmd_array;
}

#if !defined(WIN32)
#include <unistd.h>
#include <fcntl.h>
#include <errno.h>
#include <sys/types.h>
#include <sys/wait.h>

// unix pipelined execution code
bool pipeExec( char** cmd_argv, const string& command, const string& input, string& output, string& error ){
        int stdin_pipe[2], stdout_pipe[2], stderr_pipe[2];
        boolean success = false;
        pid_t sid;
        pid_t pid1;
        const char* fail;
        char buf[1024];
        ssize_t bread = 0;
        int rval = 0;

        if((sid = setsid()) < 0)sid = getpgrp();
        
        if((sid < 0 && (fail = "sid"))
         || (pipe(stdin_pipe) < 0 && (fail = "stdin"))
         || (pipe(stdout_pipe) < 0 && (fail = "stdout"))
//       || (pipe(stderr_pipe) < 0 && (fail = "stderr"))
        )
    {
                fprintf(stderr, "Ouch, the world just collapsed (%s).\n", fail);
                perror("muscle:");
                goto cleanup;
        }
        
        fcntl(stdin_pipe[0], F_SETFL, fcntl(stdin_pipe[0], F_GETFL) & ~O_NONBLOCK);
        fcntl(stdin_pipe[1], F_SETFL, fcntl(stdin_pipe[1], F_GETFL) & ~O_NONBLOCK);
        fcntl(stdout_pipe[0], F_SETFL, fcntl(stdout_pipe[0], F_GETFL) & ~O_NONBLOCK);
        fcntl(stdout_pipe[1], F_SETFL, fcntl(stdout_pipe[1], F_GETFL) & ~O_NONBLOCK);
/*      fcntl(stderr_pipe[0], F_SETFL, fcntl(stderr_pipe[0], F_GETFL) & ~O_NONBLOCK);
        fcntl(stderr_pipe[1], F_SETFL, fcntl(stderr_pipe[1], F_GETFL) & ~O_NONBLOCK);
*/      
        if((pid1 = fork()) < 0)goto cleanup;    
        if(pid1)
                setpgid(pid1, sid);
        else
        {
                dup2(stdin_pipe[0], 0);
                dup2(stdout_pipe[1], 1);
//              dup2(stderr_pipe[1], 2);
                close( stdin_pipe[0] );
                close( stdin_pipe[1] );
                close( stdout_pipe[0] );
                close( stdout_pipe[1] );
//              close( stderr_pipe[0] );
//              close( stderr_pipe[1] );
                execvp(cmd_argv[0], cmd_argv);
                _exit(errno);
        }
        rval = write( stdin_pipe[1], input.c_str(), input.size() );
        if( rval == -1 )
                perror( "write: " );
        if( close( stdin_pipe[1] ) )
                perror( "close stdin_w: " );
        if( close( stdin_pipe[0] ) )
                perror( "close stdin_r: " );

        close( stdout_pipe[1] );
        // read the alignment
        while(true){
                bzero( buf, sizeof(buf) );
                bread = read( stdout_pipe[0], buf, 1023 );
                if( bread == 0 )
                        break;  // reached EOF
                if( bread == -1 ){
                        perror("muscle read: " );
                }
                output += buf;
        }
        wait( NULL );
        success = true;
        
cleanup:
        close(stdin_pipe[0]);
        close(stdin_pipe[1]);
        close(stdout_pipe[0]);
        close(stdout_pipe[1]);
//      close(stderr_pipe[0]);
//      close(stderr_pipe[1]);
        return success;
};


#else

//windows piping code
#include <windows.h>
#define bzero(a) memset(a,0,sizeof(a)) //easier -- shortcut

bool IsWinNT()  //check if we're running NT
{
  OSVERSIONINFO osv;
  osv.dwOSVersionInfoSize = sizeof(osv);
  GetVersionEx(&osv);
  return (osv.dwPlatformId == VER_PLATFORM_WIN32_NT);
}

void ErrorMessage(char *str)  //display detailed error info
{
  LPVOID msg;
  FormatMessage(
    FORMAT_MESSAGE_ALLOCATE_BUFFER | FORMAT_MESSAGE_FROM_SYSTEM,
    NULL,
    GetLastError(),
    MAKELANGID(LANG_NEUTRAL, SUBLANG_DEFAULT), // Default language
    (LPTSTR) &msg,
    0,
    NULL
  );
  printf("%s: %s\n",str,msg);
  LocalFree(msg);
}

bool pipeExec( char** cmd_argv, const string& command, const string& input, string& output, string& error ){

        char buf[1024];           //i/o buffer

        STARTUPINFO si;
        SECURITY_ATTRIBUTES sa;
        SECURITY_DESCRIPTOR sd;               //security information for pipes
        PROCESS_INFORMATION pi;
        HANDLE newstdin_w,newstdout_w,newstderr_w,newstdin_r,newstdout_r,newstderr_r;
        HANDLE read_stdout,read_stderr,write_stdin;  //pipe handles
        boolean success = false;

        if (IsWinNT())        //initialize security descriptor (Windows NT)
        {
                InitializeSecurityDescriptor(&sd,SECURITY_DESCRIPTOR_REVISION);
                SetSecurityDescriptorDacl(&sd, true, NULL, false);
                sa.lpSecurityDescriptor = &sd;
        }
        else sa.lpSecurityDescriptor = NULL;

        sa.nLength = sizeof(SECURITY_ATTRIBUTES);
        sa.bInheritHandle = true;         //allow inheritable handles

        if (!CreatePipe(&newstdin_r,&newstdin_w,&sa,0))   //create stdin pipe
        {
                ErrorMessage("CreatePipe");
                goto finito;
        }
        if (!CreatePipe(&newstdout_r,&newstdout_w,&sa,0))  //create stdout pipe
        {
                ErrorMessage("CreatePipe");
                goto finito;
        }
        if (!CreatePipe(&newstderr_r,&newstderr_w,&sa,0))  //create stdout pipe
        {
                ErrorMessage("CreatePipe");
                goto finito;
        }
        // Duplicate the write handle to the pipe so it is not inherited. 
        boolean fSuccess = DuplicateHandle(GetCurrentProcess(), newstdin_w, 
                GetCurrentProcess(), &write_stdin, 0, 
                FALSE,                  // not inherited 
                DUPLICATE_SAME_ACCESS); 
        if (! fSuccess){
                ErrorMessage("DuplicateHandle failed"); 
                goto finito;
        }
        CloseHandle(newstdin_w); 
        newstdin_w = INVALID_HANDLE_VALUE;

        // Duplicate the read handle to the pipe so it is not inherited. 
        fSuccess = DuplicateHandle(GetCurrentProcess(), newstdout_r, 
                GetCurrentProcess(), &read_stdout, 0, 
                FALSE,                  // not inherited 
                DUPLICATE_SAME_ACCESS); 
        if (! fSuccess){
                ErrorMessage("DuplicateHandle failed"); 
                goto finito;
        }
        CloseHandle(newstdout_r); 
        newstdout_r = INVALID_HANDLE_VALUE;

        // Duplicate the read handle to the pipe so it is not inherited. 
        fSuccess = DuplicateHandle(GetCurrentProcess(), newstderr_r, 
                GetCurrentProcess(), &read_stderr, 0, 
                FALSE,                  // not inherited 
                DUPLICATE_SAME_ACCESS); 
        if (! fSuccess){
                ErrorMessage("DuplicateHandle failed"); 
                goto finito;
        }
        CloseHandle(newstderr_r); 
        newstderr_r = INVALID_HANDLE_VALUE;

        GetStartupInfo(&si);      //set startupinfo for the spawned process
        /*
                The dwFlags member tells CreateProcess how to make the process.
                STARTF_USESTDHANDLES validates the hStd* members. STARTF_USESHOWWINDOW
                validates the wShowWindow member.
        */
        si.dwFlags = STARTF_USESTDHANDLES|STARTF_USESHOWWINDOW;
        si.wShowWindow = SW_HIDE;
        si.hStdOutput = newstdout_w;
        si.hStdError = newstderr_w;     //set the new handles for the child process
        si.hStdInput = newstdin_r;

        //spawn the child process
        char* cmd = new char[ command.length() + 1 ];
        strcpy( cmd, command.c_str() );
        if (!CreateProcess(NULL,cmd,NULL,NULL,TRUE,0,
                                                NULL,NULL,&si,&pi))
        {
                delete cmd;
                ErrorMessage("CreateProcess");
                goto finito;
        }
        delete cmd;

        unsigned long exit=0;  //process exit code
        unsigned long bread;   //bytes read
        unsigned long avail;   //bytes available

        WriteFile(write_stdin, input.c_str(), input.size(), &bread, NULL); //send data to stdin
        CloseHandle(write_stdin);
        write_stdin = INVALID_HANDLE_VALUE;

        // Wait until child process exits.
        while( true ){
                GetExitCodeProcess( pi.hProcess, &exit );
                if( exit != STILL_ACTIVE )
                        WaitForSingleObject( pi.hProcess, INFINITE );
                        
                // read anything that came to stdout
                PeekNamedPipe(read_stdout,buf,1023,&bread,&avail,NULL);
                if( avail == 0 )
                        Sleep(5);       // didn't get anything, so take a break to avoid hogging the CPU...
                while( avail > 0 ){
                        bzero(buf);
                        int read_size = 1023 < avail ? 1023 : avail;
                        ReadFile(read_stdout,buf,read_size,&bread,NULL);  //read the stdout pipe
                        avail -= bread;
                        output += buf;
                }

                // read anything that came to stderr
                PeekNamedPipe(read_stderr,buf,1023,&bread,&avail,NULL);
                while( avail > 0 ){
                        bzero(buf);
                        int read_size = 1023 < avail ? 1023 : avail;
                        ReadFile(read_stderr,buf,read_size,&bread,NULL);  //read the stdout pipe
                        avail -= bread;
                        error += buf;
                }

                if( exit != STILL_ACTIVE )
                        break;
        }
        // Wait until child process exits.
    WaitForSingleObject( pi.hProcess, INFINITE );
        success = true;

        //clean up and exit
finito:
    if( pi.hThread != INVALID_HANDLE_VALUE )
                CloseHandle(pi.hThread);
    if( pi.hProcess != INVALID_HANDLE_VALUE )
                CloseHandle(pi.hProcess);
    if( newstdin_r != INVALID_HANDLE_VALUE )
                CloseHandle(newstdin_r);
    if( newstdout_w != INVALID_HANDLE_VALUE )
                CloseHandle(newstdout_w);
    if( newstderr_w != INVALID_HANDLE_VALUE )
                CloseHandle(newstderr_w);
    if( read_stdout != INVALID_HANDLE_VALUE )
                CloseHandle(read_stdout);
    if( read_stderr != INVALID_HANDLE_VALUE )
                CloseHandle(read_stderr);
    if( write_stdin != INVALID_HANDLE_VALUE )
                CloseHandle(write_stdin);
        return success;
}

#endif



MuscleInterface& MuscleInterface::getMuscleInterface()
{
        static MuscleInterface m_ci;

        return m_ci;
}

MuscleInterface::MuscleInterface() : GappedAligner() {
        muscle_path = "muscle_aed";
        muscle_arguments = "-stable -quiet -seqtype DNA";
        muscle_cmdline = parseCommand( muscle_path + " " + muscle_arguments );
        max_alignment_length = 12500;
}

void MuscleInterface::ParseMusclePath( const char* argv0 ){
        // get the execution path
        string path_str = argv0;
        // trim quotes
        if( path_str[0] == '"' )
                path_str = path_str.substr( 1, path_str.size() - 2 );
        standardizePathString( path_str );
        string::size_type i = path_str.rfind('/');
        if( i != string::npos )
                path_str.erase(i+1, path_str.length() - (i+1));
        else
                path_str.clear();
        SetMusclePath( '"' + path_str + "muscle_aed\"");
}

void MuscleInterface::SetMusclePath( const string& path ){
        muscle_path = path;
        ClearCommandLine();
        muscle_cmdline = parseCommand( muscle_path + " " + muscle_arguments );
}

void MuscleInterface::SetExtraMuscleArguments( const string& args )
{
        extra_muscle_arguments = args;
}

void MuscleInterface::SetMuscleArguments( const string& args )
{
        ClearCommandLine();
        muscle_arguments = args + " " + extra_muscle_arguments;
        muscle_cmdline = parseCommand( muscle_path + " " + args + " " + extra_muscle_arguments );
}

MuscleInterface& MuscleInterface::operator=( const MuscleInterface& ci ){
        GappedAligner::operator =( ci );
        return *this;
}

//tjt: not the best way of doing this, should have just one Align function that takes an AbstractMatch*,
        //     not both Match* & AbstractMatch* in separate, nearly identical functions..
        //     Such a change would involve changes to GappedAligner, and would require some additional care taken
        //     with SeqCount & Multiplicity, as well as seq_table[ seqI ]->length()/seq_table[ 0 ]->length(i),
        //     for now, leave like this. hopefully sooner than later, make pretty!
boolean MuscleInterface::Align( GappedAlignment& cr, Match* r_begin, Match* r_end, vector< gnSequence* >& seq_table ){
        gnSeqI gap_size = 0;
        boolean create_ok = true;
        uint seq_count = seq_table.size();
        //seq_count = r_begin->Multiplicity();
        uint seqI;
        uint align_seqs = 0;
        vector< string > tmp_mat = vector< string >( seq_count );
try{

// 
//      Get the sequence in the intervening gaps between these two matches
//
        vector< string > seq_data;
        vector< int64 > starts;
        vector< uint > seqs;
        const gnFilter* rc_filter = gnFilter::DNAComplementFilter();
        
        for( seqI = 0; seqI < seq_count; seqI++ ){

                // skip this sequence if it's undefined
                if( (r_end != NULL && r_end->Start( seqI ) == NO_MATCH ) ||
                        (r_begin != NULL && r_begin->Start( seqI ) == NO_MATCH) ){
                        starts.push_back( NO_MATCH );
                        continue;
                }

                // determine the size of the gap
                int64 gap_start = 0;
                int64 gap_end = 0;
                getInterveningCoordinates( seq_table, r_begin, r_end, seqI, gap_start, gap_end );

                int64 diff = gap_end - gap_start;
                if( diff <= 0 || diff > max_alignment_length ){
                        starts.push_back( NO_MATCH );
                        continue;       // skip this sequence if it's either too big or too small
                }
                seqs.push_back( seqI );
// the gnSequence pointers are shared across threads and have a common ifstream
                // extract sequence data
                if( r_end == NULL || r_end->Start( seqI ) > 0 ){
                        starts.push_back( gap_start );
                        seq_data.push_back( seq_table[ seqI ]->ToString( diff , gap_start ) );
                }else{
                        // reverse complement the sequence data.
                        starts.push_back( -gap_start );
                        string cur_seq_data = seq_table[ seqI ]->ToString( diff , gap_start );
                        rc_filter->ReverseFilter( cur_seq_data );
                        seq_data.push_back( cur_seq_data );
                }
        }

        if( seqs.size() <= 1 )
                create_ok = false;

        if( create_ok ){
//              SetMuscleArguments( " -quiet -stable -seqtype DNA " );
                vector< string > aln_matrix;
                if( !CallMuscleFast( aln_matrix, seq_data ) ){
                        cout << "Muscle was unable to align:\n";
                        if( r_begin )
                                cout << "Left match: " << *r_begin << endl;
                        if( r_end )
                                cout << "Right match: " << *r_end << endl;
                        return false;
                }

                gnSeqI aln_length = aln_matrix.size() == 0 ? 0 : aln_matrix[0].length();
                cr = GappedAlignment( seq_count, aln_length );
                vector< string > aln_mat = vector< string >( seq_count );

                // set sequence starts
                for( uint seqI = 0; seqI < seqs.size(); seqI++ ){
                        cr.SetLength( seq_data[ seqI ].size(), seqs[ seqI ] );
                        aln_mat[ seqs[ seqI ] ] = aln_matrix[ seqI ];
                }
                for( uint seqI = 0; seqI < seq_count; seqI++ ){
                        cr.SetStart( seqI, starts[ seqI ] );
                        if( aln_mat[ seqI ].length() != aln_length )
                                aln_mat[ seqI ] = string( aln_length, '-' );
                }

                cr.SetAlignment( aln_mat );

                return true;
        }
}catch(exception& e){
        cerr << "At: " << __FILE__ << ":" << __LINE__ << endl;
        cerr << e.what();
}
        return false;
}

static int failure_count = 0;

boolean MuscleInterface::Align( GappedAlignment& cr, AbstractMatch* r_begin, AbstractMatch* r_end, vector< gnSequence* >& seq_table ){
        gnSeqI gap_size = 0;
        boolean create_ok = true;
        //tjt: set the seq_count to a match m's multiplicity
        //     even though all components n of match m could be 
        //     less than the k sequences
        //     if n == k, then perhaps there is 1 match component per sequence
        //     if k = 1, n == repeat match multiplicity, where n >= 2
        //     
        uint seq_count = r_begin->Multiplicity();
        uint seqI;
        uint align_seqs = 0;
        vector< string > tmp_mat = vector< string >( seq_count );
try{

// 
//      Get the sequence in the intervening gaps between these two matches
//
        vector< string > seq_data;
        vector< int64 > starts;
        vector< uint > seqs;
        const gnFilter* rc_filter = gnFilter::DNAComplementFilter();
        
        //std::cout << "getting regions between match components to align" << std::endl;
        for( seqI = 0; seqI < seq_count; seqI++ ){

                // skip this sequence if it's undefined
                if( (r_end != NULL && r_end->Start( seqI ) == NO_MATCH ) ||
                        (r_begin != NULL && r_begin->Start( seqI ) == NO_MATCH) ){
                        starts.push_back( NO_MATCH );
                        continue;
                }

                // determine the size of the gap
                int64 gap_start = 0;
                int64 gap_end = 0;
                // skip this sequence if it's undefined
                if( (r_end != NULL && r_end->Start( seqI ) == NO_MATCH) ||
                        (r_begin != NULL && r_begin->Start( seqI ) == NO_MATCH) ){
                        gap_start = 0;
                        gap_end = 0;
                        return true;
                }
                                
                // determine the size of the gap
                gap_end = r_end != NULL ? r_end->Start( seqI ) : seq_table[ seqI ]->length() + 1;
                gap_start = r_begin != NULL ? r_begin->End( seqI ) + 1 : 1;
                if( gap_end < 0 || gap_start < 0 ){
                        gap_end   = r_begin != NULL ? -r_begin->Start( seqI ) : seq_table[ 0 ]->length() + 1;
                        gap_start = r_end != NULL ? -r_end->Start( seqI ) + r_end->Length( seqI ) : 1;
                }
                if( gap_end <= 0 || gap_start <= 0 ){
                        // if either is still < 0 then there's a problem...
                        genome::ErrorMsg( "Error constructing intervening coordinates" );
                }

                int64 diff = gap_end - gap_start;
                
                //diff <= 0 ||
                if( diff <= 0 || diff > max_alignment_length ){
                        starts.push_back( NO_MATCH );
                        continue;       // skip this sequence if it's either too big or too small
                }
                seqs.push_back( seqI );

                // extract sequence data
                if (0 )
                {
                        starts.push_back( gap_start );
                        seq_data.push_back( "A" );
                        std::cout << "A" << std::endl;
                        diff = 1;
                }
// the gnSequence pointers are shared across threads and have a common ifstream
                if( r_end == NULL || r_end->Start( seqI ) > 0 ){
                        starts.push_back( gap_start );
                        //std::cout << seq_table[ 0 ]->ToString( diff , gap_start ) << std::endl;
                        //tjt: all sequences are concatenated together into 1 seq_table entry
                        //
                        seq_data.push_back( seq_table[ 0 ]->ToString( diff , gap_start ) );
                }else{
                        // reverse complement the sequence data.
                        starts.push_back( -gap_start );
                        //tjt: all sequences are concatenated together into 1 seq_table entry
                        //     
                        string cur_seq_data = seq_table[ 0 ]->ToString( diff , gap_start );
                        rc_filter->ReverseFilter( cur_seq_data );
                        seq_data.push_back( cur_seq_data );
                        //std::cout << cur_seq_data << std::endl;
                }
        }

        if( seqs.size() <= 1 )
                create_ok = false;

        if( create_ok ){
//              SetMuscleArguments( " -quiet -stable -seqtype DNA " );
                vector< string > aln_matrix;
                if( !CallMuscleFast( aln_matrix, seq_data ) ){
                        cout << "Muscle was unable to align:\n";
                        //if( r_begin )
                        //      cout << "Left match: " << *r_begin << endl;
                        //if( r_end )
                        //      cout << "Right match: " << *r_end << endl;
                        return false;
                }

                gnSeqI aln_length = aln_matrix.size() == 0 ? 0 : aln_matrix[0].length();
                cr = GappedAlignment( seq_count, aln_length );
                vector< string > aln_mat = vector< string >( seq_count );

                // set sequence starts
                for( uint seqI = 0; seqI < seqs.size(); seqI++ ){
                        cr.SetLength( seq_data[ seqI ].size(), seqs[ seqI ] );
                        aln_mat[ seqs[ seqI ] ] = aln_matrix[ seqI ];
                }
                for( uint seqI = 0; seqI < seq_count; seqI++ ){
                        cr.SetStart( seqI, starts[ seqI ] );
                        if( aln_mat[ seqI ].length() != aln_length )
                                aln_mat[ seqI ] = string( aln_length, '-' );
                }

                cr.SetAlignment( aln_mat );

                return true;
        }
}catch(exception& e){
        cerr << "At: " << __FILE__ << ":" << __LINE__ << endl;
        cerr << e.what();
}
        return false;
}

boolean MuscleInterface::CallMuscle( vector< string >& aln_matrix, const vector< string >& seq_table )
{
        gnSequence seq;

        try{
                ostringstream input_seq_stream;
                //istringstream muscle_input_seq_stream;
                for( uint seqI = 0; seqI < seq_table.size(); seqI++ ){
                        seq += seq_table[ seqI ];
                        seq.setContigName( seqI, "seq" );
                }
                gnFASSource::Write( seq, input_seq_stream, false, true );               
                // now open a pipe to Muscle
                string muscle_cmd = muscle_path + " " + muscle_arguments;
                string output;
                string error;
                boolean success = pipeExec( muscle_cmdline, muscle_cmd, input_seq_stream.str(), output, error );
                if( !success || output.size() == 0 )
                {
                        throw "b0rk3d";
                }

                istringstream output_aln_stream( output );
                string cur_line;

                // parse the fasta output
                while( getline( output_aln_stream, cur_line ) ){
                        if( cur_line[0] == '>' ){
                                aln_matrix.push_back( "" );
                                continue;
                        }
                        gnSeqI len = cur_line.size();
                        len = cur_line[ len - 1 ] == '\r' ? len - 1 : len;
                        uint seqI = aln_matrix.size() - 1;
                        aln_matrix[ seqI ] += cur_line.substr( 0, len );
                }

                return true;
        }catch( gnException& gne ){
        }catch( exception& e ){
        }catch(...){
        }
        cerr << "muscle failed!  saving failed input data to muscle_failure_" << failure_count << ".txt\n";
        cerr << "Please contact the Mauve developers about this problem\n";
        stringstream debug_fname;
        debug_fname << "muscle_failure_" << failure_count++ << ".txt";
        ofstream debug_file( debug_fname.str().c_str() );
        gnFASSource::Write(seq, debug_file, false);
        debug_file.close();
        return false;
}

// version 2 of this code: attempt to call muscle without performing costly disk I/O!!
boolean MuscleInterface::CallMuscleFast( vector< string >& aln_matrix, const vector< string >& seq_table )
{
        g_SeqType.get() = SEQTYPE_DNA;  // we're operating on DNA
        g_uMaxIters.get() = 1;                  // and we don't want to refine the alignment...yet
        g_bStable.get() = true;                 // we want output seqs in the same order as input
        g_bQuiet.get() = true;                  // and don't print anything to the console
        g_SeqWeight1.get() = SEQWEIGHT_ClustalW;        // not sure what weighting scheme works best for DNA

        SetMaxIters(g_uMaxIters.get());
        SetSeqWeightMethod(g_SeqWeight1.get());

        // now construct a SeqVect containing input sequences
        SeqVect sv;
        const char* seqname = "seq00000";
        for( size_t seqI = 0; seqI < seq_table.size(); seqI++ )
        {
                Seq curseq;
                curseq.SetId(seqI);
                curseq.SetName(seqname);
                curseq.resize(seq_table[seqI].size());
                std::copy(seq_table[seqI].begin(), seq_table[seqI].end(), curseq.begin());
                sv.AppendSeq(curseq);
        }

        MSA msaTmp;
        MUSCLE(sv,msaTmp);

        // now extract the alignment
        aln_matrix.clear();
        aln_matrix.resize(msaTmp.GetSeqCount());
        for( size_t seqI = 0; seqI < msaTmp.GetSeqCount(); seqI++ )
        {
                unsigned indie = msaTmp.GetSeqIndex(seqI);
                const char* buf = msaTmp.GetSeqBuffer(indie);
                string curseq(buf, msaTmp.GetColCount());
                swap(aln_matrix[seqI],curseq);
        }
        return true;    // how can it possibly fail? :)
}

bool MuscleInterface::Refine( GappedAlignment& ga, size_t windowsize )
{
        const vector< string >& seq_table = GetAlignment( ga, vector< gnSequence* >() );
        vector< string > aln_table;
        for( uint seqI = 0; seqI < ga.SeqCount(); seqI++ )
        {
                if( ga.LeftEnd(seqI) != NO_MATCH )
                {
                        aln_table.push_back( seq_table[seqI] );
                }
        }
        vector< string > aln_matrix;
        if( windowsize == 0 )
                SetMuscleArguments( " -quiet -refine -seqtype DNA " );
        else
        {
                stringstream sstr;
                sstr << " -quiet -seqtype DNA -refinew -refinewindow " << windowsize << " ";
                SetMuscleArguments( sstr.str() );
        }
        bool success = CallMuscle( aln_matrix, aln_table );
        if( success )
        {
                aln_table.clear();
                uint alnI = 0;
                for( uint seqI = 0; seqI < ga.SeqCount(); seqI++ )
                {
                        if( ga.LeftEnd(seqI) != NO_MATCH )
                                aln_table.push_back( aln_matrix[alnI++] );
                        else
                                aln_table.push_back( string( aln_matrix[0].size(), '-' ) );
                }
                ga.SetAlignment( aln_table );
        }
        return success;
}

void msaFromSeqTable(MSA& msa, const vector< string >& seq_table, unsigned id_base = 0)
{
        msa.SetSize(seq_table.size(), seq_table[0].size());
        for( uint seqI = 0; seqI < seq_table.size(); seqI++ )
        {
                stringstream ss;
                ss << "seq" << seqI;
                msa.SetSeqName(seqI, ss.str().c_str());
                msa.SetSeqId(seqI,seqI+id_base);
                for(size_t i = 0; i < seq_table[seqI].size(); i++)
                        msa.SetChar(seqI, i, seq_table[seqI][i]);
        }
}


bool MuscleInterface::RefineFast( GappedAlignment& ga, size_t windowsize )
{
        const vector< string >& seq_table = GetAlignment( ga, vector< gnSequence* >() );
        vector< string > aln_table;
        for( uint seqI = 0; seqI < ga.SeqCount(); seqI++ )
        {
                if( ga.LeftEnd(seqI) != NO_MATCH )
                {
                        aln_table.push_back( seq_table[seqI] );
                }
        }

        g_SeqType.get() = SEQTYPE_DNA;  // we're operating on DNA
        g_uMaxIters.get() = 1;                  // and we don't want to refine the alignment...yet
        g_bStable.get() = true;                 // we want output seqs in the same order as input
        g_bQuiet.get() = true;                  // and don't print anything to the console
        g_SeqWeight1.get() = SEQWEIGHT_ClustalW;        // not sure what weighting scheme works best for DNA

        g_uRefineWindow.get() = windowsize;
        g_uWindowTo.get() = 0;

        SetMaxIters(g_uMaxIters.get());
        SetSeqWeightMethod(g_SeqWeight1.get());

        MSA::SetIdCount(seq_table.size());

        // create an MSA
        MSA msa;
        msaFromSeqTable(msa, seq_table);

        SetAlpha(ALPHA_DNA);
        msa.FixAlpha();
        SetPPScore(PPSCORE_SPN);
        SetMuscleInputMSA(msa);

        Tree GuideTree;
        TreeFromMSA(msa, GuideTree, g_Cluster2.get(), g_Distance2.get(), g_Root2.get());
        SetMuscleTree(GuideTree);

        MSA msaOut;
        MSA* finalMsa;

        if(windowsize == 0)
        {
                if (g_bAnchors.get())
                        RefineVert(msa, GuideTree, g_uMaxIters.get());
                else
                        RefineHoriz(msa, GuideTree, g_uMaxIters.get(), false, false);
                finalMsa = &msa;
        }else{
                RefineW(msa, msaOut);
                finalMsa = &msaOut;
        }


        ValidateMuscleIds(*finalMsa);
        ValidateMuscleIds(GuideTree);

        // now extract the alignment
        vector< string > aln_matrix;
        aln_matrix.resize(finalMsa->GetSeqCount());
        for( size_t seqI = 0; seqI < finalMsa->GetSeqCount(); seqI++ )
        {
                unsigned indie = finalMsa->GetSeqIndex(seqI);
                const char* buf = finalMsa->GetSeqBuffer(indie);
                string curseq(buf, finalMsa->GetColCount());
                swap(aln_matrix[seqI],curseq);
        }

        ga.SetAlignment( aln_matrix );
        return true;
}


void stripGapColumns( std::vector< std::string >& aln )
{
        size_t cur_col = 0;
        size_t gap_seq = 0;
        for( size_t colI = 0; colI < aln[0].size(); colI++ )
        {
                gap_seq = 0;
                for( ; gap_seq < aln.size(); gap_seq++ )
                        if( aln[gap_seq][colI] != '-' )
                                break;
                if( gap_seq != aln.size() )
                {
                        for( gap_seq = 0; gap_seq < aln.size(); gap_seq++ )
                                aln[gap_seq][cur_col] = aln[gap_seq][colI];
                        cur_col++;
                }
        }
        for( gap_seq = 0; gap_seq < aln.size(); gap_seq++ )
                aln[gap_seq].resize(cur_col);
}

void stripGaps( std::string& str )
{
        std::string::iterator striter = std::remove(str.begin(), str.end(), '-');
        str.resize(striter - str.begin());
}

bool MuscleInterface::ProfileAlign( const GappedAlignment& ga1, const GappedAlignment& ga2, GappedAlignment& aln, bool anchored )
{
        try{
                const vector< string >& aln1 = GetAlignment( ga1, vector< gnSequence* >() );
                const vector< string >& aln2 = GetAlignment( ga2, vector< gnSequence* >() );
                vector< uint > order;
                ostringstream input_seq_stream;
                gnSequence seq;
                vector< string > aln11( ga1.Multiplicity() );
                vector< string > aln22( ga2.Multiplicity() );
                size_t curI = 0;
                for( uint seqI = 0; seqI < aln1.size(); seqI++ )
                {
                        if( ga1.LeftEnd(seqI) != NO_MATCH )
                        {
                                aln11[curI++] = aln1[seqI];
                                order.push_back(seqI);
                        }
                }
                curI = 0;
                for( uint seqI = 0; seqI < aln2.size(); seqI++ )
                {
                        if( ga2.LeftEnd(seqI) != NO_MATCH )
                        {
                                aln22[curI++] = aln2[seqI];
                                order.push_back(seqI);
                        }
                }
// strip the gap columns only if we're doing unanchored PP alignment
                if( !anchored )
                {
                        stripGapColumns(aln11);
                        stripGapColumns(aln22);
                }
                for( uint seqI = 0; seqI < aln11.size(); seqI++ )
                {
                        seq += aln11[ seqI ];
                        seq.setContigName( seq.contigListLength()-1, "seq" );
                }

                gnFASSource::Write( seq, input_seq_stream, false, true );
                input_seq_stream << "=\n";

                gnSequence seq2;
                for( uint seqI = 0; seqI < aln22.size(); seqI++ )
                {
                        seq2 += aln22[ seqI ];
                        seq2.setContigName( seq2.contigListLength()-1, "seq" );
                }

                gnFASSource::Write( seq2, input_seq_stream, false, true );
                input_seq_stream << "=\n";

                if( debug_muscle )
                {
                        // for debugging: write the anchored profiles to a file
                        stringstream debug_fname;
                        debug_fname << "muscle_debug_" << failure_count++ << ".txt";
                        ofstream debug_file( debug_fname.str().c_str() );
                        debug_file << input_seq_stream.str();
                        debug_file.close();
                }

                // now open a pipe to Muscle
                string musc_args = "-quiet -seqtype DNA -profile -ProfileOnStdIn ";
                if( anchored )
                        musc_args += "-AnchoredPP ";
                SetMuscleArguments( musc_args );
                string output;
                string error;
                string muscle_cmd = muscle_path + " " + muscle_arguments;
                if( debug_muscle )
                {
                        cerr << "Running " << muscle_cmd << endl;
                }
                boolean success = pipeExec( muscle_cmdline, muscle_cmd, input_seq_stream.str(), output, error );
                if( !success || output.size() == 0 )
                {
                        if( output.size() == 0 )
                                cerr << "\nmuscle nothing\n";
                        else
                                cerr << "\nunsuccessful muscle\n";
                        return false;
                }

                istringstream output_aln_stream( output );
                string cur_line;

                // parse the fasta output
                vector< string > aln_matrix( ga1.SeqCount() );
                int ordI = -1;
                while( getline( output_aln_stream, cur_line ) ){
                        if( cur_line[0] == '>' ){
                                ordI++;
                                continue;
                        }
                        gnSeqI len = cur_line.size();
                        len = cur_line[ len - 1 ] == '\r' ? len - 1 : len;
                        uint seqI = aln_matrix.size() - 1;
                        aln_matrix[ order[ordI] ] += cur_line.substr( 0, len );
                }
                for( size_t i = 0; i < aln_matrix.size(); i++ )
                {
                        if( aln_matrix[i].size() == 0 )
                                aln_matrix[i].resize( aln_matrix[order[0]].size(), '-' );
                }

                aln.SetAlignment( aln_matrix );
                for( uint seqI = 0; seqI < ga1.SeqCount(); seqI++ )
                        if( ga1.LeftEnd(seqI) != NO_MATCH )
                        {
                                aln.SetLeftEnd(seqI, ga1.LeftEnd(seqI));
                                aln.SetLength(ga1.Length(seqI), seqI);
                        }
                for( uint seqI = 0; seqI < ga2.SeqCount(); seqI++ )
                        if( ga2.LeftEnd(seqI) != NO_MATCH )
                        {
                                aln.SetLeftEnd(seqI, ga2.LeftEnd(seqI));
                                aln.SetLength(ga2.Length(seqI), seqI);
                        }
                return true;
        }catch( gnException& gne ){
        }catch( exception& e ){
        }catch(...){
        }
        return false;
}


bool MuscleInterface::ProfileAlignFast( const GappedAlignment& ga1, const GappedAlignment& ga2, GappedAlignment& aln, bool anchored )
{
        try{
                const vector< string >& aln1 = GetAlignment( ga1, vector< gnSequence* >() );
                const vector< string >& aln2 = GetAlignment( ga2, vector< gnSequence* >() );
                vector< uint > order;
                vector< string > aln11( ga1.Multiplicity() );
                vector< string > aln22( ga2.Multiplicity() );
                size_t curI = 0;
                for( uint seqI = 0; seqI < aln1.size(); seqI++ )
                {
                        if( ga1.LeftEnd(seqI) != NO_MATCH )
                        {
                                aln11[curI++] = aln1[seqI];
                                order.push_back(seqI);
                        }
                }
                curI = 0;
                for( uint seqI = 0; seqI < aln2.size(); seqI++ )
                {
                        if( ga2.LeftEnd(seqI) != NO_MATCH )
                        {
                                aln22[curI++] = aln2[seqI];
                                order.push_back(seqI);
                        }
                }
// strip the gap columns only if we're doing unanchored PP alignment
                if( !anchored )
                {
                        stripGapColumns(aln11);
                        stripGapColumns(aln22);
                }

                g_SeqType.get() = SEQTYPE_DNA;  // we're operating on DNA
                g_uMaxIters.get() = 1;                  // and we don't want to refine the alignment...yet
                g_bStable.get() = true;                 // we want output seqs in the same order as input
                g_bQuiet.get() = true;                  // and don't print anything to the console
                g_SeqWeight1.get() = SEQWEIGHT_ClustalW;        // not sure what weighting scheme works best for DNA

                SetMaxIters(g_uMaxIters.get());
                SetSeqWeightMethod(g_SeqWeight1.get());

                MSA::SetIdCount(order.size());

                MSA msa1;
                MSA msa2;
                MSA msaOut;
                msaFromSeqTable(msa1, aln11);
                msaFromSeqTable(msa2, aln22, msa1.GetSeqCount());

                SetAlpha(ALPHA_DNA);
                msa1.FixAlpha();
                msa2.FixAlpha();
                SetPPScore(PPSCORE_SPN);

                if(anchored)
                {
                        AnchoredProfileProfile(msa1, msa2, msaOut);
                }else{
                        ProfileProfile(msa1, msa2, msaOut);
                }

                // get the output
                vector< string > aln_matrix( aln1.size() );
                for( size_t seqI = 0; seqI < msaOut.GetSeqCount(); seqI++ )
                {
                        unsigned indie = msaOut.GetSeqIndex(seqI);
                        const char* buf = msaOut.GetSeqBuffer(indie);
                        string curseq(buf, msaOut.GetColCount());
                        swap(aln_matrix[order[indie]],curseq);

                        // debugging, check that sequences came out in the same order they went in!
/*                      string inseq = aln1[order[indie]];
                        string outseq = aln_matrix[order[indie]];
                        stripGaps(inseq);
                        stripGaps(outseq);
                        if(inseq != outseq)
                        {
                                unsigned indie = msaOut.GetSeqIndex(seqI);
                                cerr << "bad indie " << indie << endl;
                                genome::breakHere();
                        }
*/
                }
                // fill empty seqs with gaps
                for( size_t seqI = 0; seqI < aln_matrix.size(); seqI++ )
                        if(aln_matrix[seqI].size() == 0)
                                aln_matrix[seqI].resize(msaOut.GetColCount(), '-');

                aln.SetAlignment( aln_matrix );
                for( uint seqI = 0; seqI < ga1.SeqCount(); seqI++ )
                        if( ga1.LeftEnd(seqI) != NO_MATCH )
                        {
                                aln.SetLeftEnd(seqI, ga1.LeftEnd(seqI));
                                aln.SetLength(ga1.Length(seqI), seqI);
                        }
                for( uint seqI = 0; seqI < ga2.SeqCount(); seqI++ )
                        if( ga2.LeftEnd(seqI) != NO_MATCH )
                        {
                                aln.SetLeftEnd(seqI, ga2.LeftEnd(seqI));
                                aln.SetLength(ga2.Length(seqI), seqI);
                        }
                return true;

        }catch( gnException& gne ){
        }catch( exception& e ){
        }catch(...){
        }
        return false;
}


void MuscleInterface::CreateTree( const NumericMatrix<double>& distances, const std::string& tree_filename  )
{
        DistFunc df;
        df.SetCount( distances.rows() );
        for( size_t i = 0; i < distances.rows(); i++ )
                for( size_t j = 0; j < distances.rows(); j++ )
                        df.SetDist( i, j, distances(i,j) );

        for( size_t i = 0; i < distances.rows(); i++ )
        {
                stringstream ss;
                ss << "seq";
                ss << i + 1;
                df.SetName( i, ss.str().c_str() );
                df.SetId( i, i );
        }
        ClustSetDF csdf( df );
        Clust crusty;
        crusty.Create( csdf, CLUSTER_NeighborJoining );
        Tree tt;
        tt.FromClust( crusty );
        TextFile tf( tree_filename.c_str(), true );
        tt.ToFile( tf );
}


}

---