EnsemblServer/SSAHAClient.cpp

Go to the documentation of this file.

// #######################################################################

// SSAHA : Sequence Search and Alignment by Hashing Algorithm
// Version 3.2, released 1st March 2004
// Copyright (c) Genome Research 2002

// SSAHA is free software; you can redistribute it and/or modify 
// it under the terms of version 2 of the GNU General Public Licence
// as published by the Free Software Foundation.
 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public Licence for more details.
 
// You should have received a copy of the GNU General Public Licence
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
// or see the on-line version at http://www.gnu.org/copyleft/gpl.txt

// #######################################################################

// Module Name  : SSAHAClient
// File Name    : SSAHAClient.cpp
// Language     : C++
// Module Author: Anthony J. Cox (ac2@sanger.ac.uk)

// Includes:

#include "ClientServerUtils.h"
#include "SSAHAClient.h"
#include "SequenceReaderFasta.h"
#include "SequenceReaderString.h"
#include "MatchStore.h"
#include "MatchAligner.h"
#include <string.h>
#include <iomanip>
#include <strstream>


// ### Function Definitions ###

static QueryHeader qinfo;
static  Handshake hello;

MatchRemote::MatchRemote( vector<QueryInfo>& query ) :
  query_(query) {}  

MatchRemote::~MatchRemote() {}


string MatchRemote::getQueryName( void ) const {
  return query_[ data_.queryNum-1 ].second;
} 

int MatchRemote::getQuerySize( void ) const {
  WordSequence& thisSeq(query_[ data_.queryNum-1 ].first);
  return (thisSeq.size()==0) 
    ? 0 
    : (thisSeq.size()-1)*hello.wordLength + thisSeq.getNumBasesInLast();
}


void SourceReaderDummy::extractSource ( char** pSource, //vector<char>& source, 
                                                          SequenceNumber seqNum,
                                                          SequenceOffset seqStart,
                                                          SequenceOffset seqEnd ) { 
  *pSource = (char*) source_.c_str(); 
}

void sendQuery(FILE *fp, int sockfd, SequenceReader& seqReader) {

  char  recvline[MAXLINE];
  WordSequence seq;

  // get handshake from server
  // initial time out set to 12 seconds
  SocketInterface socket( sockfd, 12 );
  
  socket.receiveStruct(&hello);

  socket.checkSocketEmpty();

  cerr << "Server is using hash table with " 
       << ((hello.tableType==e2bitDNA)?gBaseBits:gResidueBits)
       << " bits per symbol, " << hello.wordLength 
       << " symbols per word.\n";

  if ((hello.tableType==e2bitDNA)&&(qinfo.bitsPerSymbol==gResidueBits)) {
    cerr << "Error: can't run a protein query against a DNA database.\n";
    throw SSAHAException("Can't run protein query against DNA database");
  }
  
  if (hello.ssahaversion != SERV_VERSION) {
    cerr << "Error: SSAHA client/Server version mismatch.\n";
    throw SSAHAException("Can't continue: client/server version mismatch");
  }
  
  cerr << "Server will reject queries of more than "  << hello.maxBufferSize << " words in total.\n";

  if ((hello.tableType!=e2bitDNA) && (qinfo.bitsPerSymbol==gBaseBits)) {
    cerr << "Sending DNA query to protein/translated database, changing word length to " << gMaxBasesPerWord << ".\n";
    hello.wordLength=gMaxBasesPerWord;
  } 

  if (qinfo.numRepeats>hello.wordLength) {
    cerr << "Warning: only repeats of " << hello.wordLength << " bases or less can be masked for this data, proceeding" << "using this value.\n";
    qinfo.numRepeats = hello.wordLength;
  }

  cerr << "Server will attempt to screen for tandem repeats of "  << qinfo.numRepeats << " bases or less.\n";

  if (qinfo.maxInsert>=hello.wordLength) {
    cerr << "Warning: indels of up to " << hello.wordLength-1  << " bases only can be handled for this data\n";
    qinfo.maxInsert = hello.wordLength-1;
  }

  if ((qinfo.sortMode == eSortAndReturnSequence) && (qinfo.bandExtension > 15)) {
    cerr << "Warning: restricting band size for banded dynamic programming to 15\n";
    qinfo.bandExtension=15;
  }

  cerr << "Matches can contain up to "  << qinfo.maxInsert << " indels between successive hits.\n";

  cerr << "Matches can contain gaps of up to "  << qinfo.maxGap << " bases.\n";

  // read queries to local memory
  const QueryInfo dummy;
  vector<QueryInfo> query; 
  query.push_back(dummy);
  qinfo.numQueryWords = 0;
  
  SequenceReaderModeReplace mode((qinfo.bitsPerSymbol==gBaseBits)?'A':'X');
  seqReader.changeMode( &mode );
  
  while ( seqReader.getNextSequence( query.back().first, hello.wordLength ) != -1 ) {
    seqReader.getLastSequenceName( query.back().second );
    qinfo.numQueryWords += query.back().first.size();
    query.push_back(dummy);
  }
  query.pop_back();
  qinfo.numQuerySeqs = query.size();

  socket.sendStruct(&qinfo);

  for( vector<QueryInfo>::iterator i(query.begin()) ; i != query.end() ; ++i )  {
    socket.sendSequence(i->first);
    // cerr << "|" << printWord(i->first,12) << "|\n";
    // WordSequence w;
    // w = i->first;
  }


  MatchHeader response;
  // wait for 60 seconds for initial response
  socket.setTimeOut(600);
  socket.receiveStruct(&response);

  if (response.wasSuccessful==false)  throw NetworkException("Query request failed!!");

  cerr << "Expecting to receive " << response.numMatches << " matches among " << response.numSubjectNames << " subject sequences.\n";


  MatchStoreRemote results(query);
  SourceReaderDummy subjectReader;

  int numColumns(80);

  MatchAligner* pAligner;

  if (qinfo.bitsPerSymbol == gBaseBits) {
    if ( hello.tableType == e2bitDNA ) {
        pAligner = new MatchAlignerDNA
          ( numColumns, qinfo.bandExtension ); 
      } else if ( hello.tableType == e5bitProtein ) {
        pAligner = new MatchAlignerTranslatedProtein ( false, numColumns, qinfo.bandExtension );   
      } else {
        // ( hello.TableType == e5bitTranslatedDNA)
        pAligner = new MatchAlignerTranslatedDNA
          ( numColumns, qinfo.bandExtension );   
      }
  }  else { // query is protein
      if (hello.tableType == e5bitProtein) {
        pAligner = new MatchAlignerProtein ( numColumns, qinfo.bandExtension ); 
      } else {
       // table is translated DNA
        pAligner = new MatchAlignerTranslatedProtein ( true, numColumns, qinfo.bandExtension );   
      } 
  }

  assert( pAligner!=false);

  // ownership of *pAligner passed to aligner
  MatchTaskAlign aligner ( seqReader, subjectReader, pAligner, false, qinfo.sortMode==eSortAndReturnSequence );

  pair<SequenceNumber,std::string> p;
  for ( int i(0) ;  i < response.numSubjectNames ; i++ ) {
       socket.receiveStruct(&p.first);
       socket.receiveString(p.second);
       results.match_.names_.insert(p);
  } 
  
  cout << setprecision(2) << setiosflags(ios::fixed);

  cout << "OK: " << response.numMatches << " " << response.numSubjectNames << endl;

  for ( int i(0); i < response.numMatches ; i++ ) {
    socket.receiveStruct(&results.match_.data_);
    if (qinfo.sortMode == eSortAndReturnSequence ) {
      socket.receiveString( subjectReader.source_ );
    }
    aligner( results );
  } 
  socket.checkSocketEmpty();
} 



int main(int numArgs, char* args[] ) {

try {
    int sockfd;
    struct sockaddr_in  servaddr;
    int fred;

    // only needed if we need to make our own fasta file
    ostrstream buf;
    istream istr(buf.rdbuf());

    if(( numArgs<11)||( numArgs>13))
    {
        cerr 
        <<"syntax: " << args[0] 
        <<" serverMachine serverPort minMatchSize maxGap maxInsert numRepeats\n"
        <<"queryType clipThreshold maxMatches sortMode "
        << " [substituteThreshold] [bandExtension]\n"
        <<"serverMachine: name of machine on which server is running\n"
        <<"serverPort   : port number on that machine\n"
        <<"minMatchSize : matches must contain at least this many matching symbols\n"
        <<"maxGap       : matches may contain gaps of up to this many symbols\n"
        <<"maxInsert    : max number of indels between successive hits in a match\n"
        <<"numRepeats   : screen out repeating motifs of up to this many symbols\n"
        <<"queryType    : DNA or protein\n"
        <<"clipThreshold: ignore words occurring this many more times than expected\n"
        <<"maxNumMatches: at most this many matches for each query sequence\n"
        <<"               (set to 0 to obtain all matches unsorted)\n"
        <<"sortMode     : none, size, percent or align\n"
        <<"               (ignored if maxNumMatches set to 0)\n"
        <<"substituteThreshold : allow 1 base/amino mismatch in words that occur\n"
        <<"               up to this many more times than expected (default 0)\n"
        <<"bandExtension: band size for banded dynamic programming extends this\n"
        <<"               far from the diagonal (ignored unless sortMode set to\n"
        <<"               align, defaults to 0)\n"; 
                                         
        throw SSAHAException("Invalid command line input to client");
    }
    memset((void*)&qinfo, 0, sizeof(qinfo));
    
    sockfd = Socket(AF_INET, SOCK_STREAM, 0);

    int portNumber(atoi(args[2]));

    cerr << "Queries will be sent via port number " << portNumber << ".\n";

    qinfo.minPrint=atoi(args[3]);

    cerr << "Only matches greater than " << qinfo.minPrint << " bases will be reported.\n";

    qinfo.maxGap=atoi(args[4]);
    qinfo.maxInsert=atoi(args[5]);
    qinfo.numRepeats=atoi(args[6]);
    qinfo.clipThreshold=atoi(args[8]);
    qinfo.maxMatches=atoi(args[9]);

    string queryType(args[7]);
    SequenceEncoder* pEncoder;

    if (queryType=="DNA")  {
      qinfo.bitsPerSymbol = gBaseBits;
      pEncoder = new SequenceEncoderDNA(12);
      cerr << "Query sequence is DNA.\n";
    } else if (queryType=="protein") { 
      qinfo.bitsPerSymbol = gResidueBits;
      pEncoder = new SequenceEncoderProtein(5);
      cerr << "Query sequence is protein.\n";
    } else {
      cerr  << "Unknown value for query type (must be \"DNA\" or \"protein\")\n";
      throw SSAHAException("Invalid value for queryType");
    }

    string sortMode(args[10]);
    if (sortMode=="none")  {
      qinfo.maxMatches=0;
    }
    else if (sortMode=="size")  {
      qinfo.sortMode = eSortByMatchLength;
    } else if (sortMode=="percent") {
      qinfo.sortMode = eSortByPercentMatch;
    } else if (sortMode=="align") {
      qinfo.sortMode = eSortAndReturnSequence;
    } else {
      cerr  << "Unknown value for sort mode " << "(must be \"none\", \"size\", \"percent\" or \"align\")\n";
      throw SSAHAException("Invalid value for sortMode");
    }

    qinfo.substituteThreshold = ((numArgs>11)?atoi(args[11]):0);
    qinfo.bandExtension = ((numArgs>12)?atoi(args[12]):0);

    while (cin.peek()==' ') cin.ignore(); // zap any leading spaces

    SequenceReader* pReader;
    string source;

    if (cin.peek()=='>')  {
      cerr  << "First nonspace character is a \">\", assuming input text " << "is in fasta format.\n";
      buf << cin.rdbuf();
      pReader = new SequenceReaderFile(istr, '>', '>', pEncoder, cerr );
    }  else {
      cerr << "Assuming input text is a plain string of " << queryType << " data.\n";
      buf << ">unnamedQuery\n" << cin.rdbuf();
      pReader = new SequenceReaderFile(istr, '>', '>', pEncoder, cerr );
    }

    memset((void*)&servaddr, 0, sizeof(servaddr));
    servaddr.sin_family = AF_INET;
    servaddr.sin_port = htons(portNumber);

    cerr << "Server is assumed to be running on machine " << args[1] << ".\n";

    struct hostent *hp;
    hp =  gethostbyname(args[1]);
    if (hp==NULL) throw NetworkException("Invalid host name");

      memcpy(&(servaddr.sin_addr.s_addr), *(hp->h_addr_list), sizeof(struct
      in_addr));

    Connect(sockfd, (SA *) &servaddr, sizeof(servaddr));
    
    sendQuery(stdin, sockfd, *pReader);

    exit(0);

} 
catch (const NetworkException& err ) {
    cerr << "Caught NetworkException: " << err.what() << "\n";
    if (((string)err.what()).substr(0,16)=="getAtLeast error") {
      cout << "ERROR: lost connection, please resubmit your query"  << endl;
    }  else {
      cout << "ERROR: " << err.what() << endl;
    } 
    exit(1);
  }
  catch (const SSAHAException& err ) {
    cerr << "Caught SSAHA exception: " << err.what() << "\n";
    cout << "ERROR: " << err.what() << endl;
    exit(1);
  }
  catch (const std::exception& err ) {
    cerr << "Caught exception: " << err.what() << "\n";
    cout << "ERROR: " << err.what() << endl;
    exit(1);
  }

}

// End of file SSAHAClient.cpp

---