EnsemblServer/OrderOrient.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  : OrderOrient
// File Name    : OrderOrient.cpp
// Language     : C++
// Module Author: Anthony J. Cox (ac2@sanger.ac.uk)

// Description:
// OrderOrient
// This is an adapted version of SSAHAClient.cpp
// The contigs for a given clone are input in fasta format via standard 
// input. The client then links to a SSAHA server of reads and interprets 
// the matches returned to determine the order and orientation of the
// contigs. To do this it needs to know which read in the server hash
// table pairs to which. This information is held in an ASCII file
// with the following format
// numberOfPairForRead1  clipStartForRead1 clipEndForRead1 \n
// numberOfPairForRead2  clipStartForRead2 clipEndForRead2 \n
// ...
// numberOfPairForReadN  clipStartForReadN clipEndForReadN \n
// Reads which don't have a pair in the database have a pair number of zero.


// Includes:

#include "ClientServerUtils.h"
//#include "OrderOrient.h"
#include "SequenceReaderFasta.h"
#include "SequenceReaderString.h"
#include <string.h>
#include <iomanip>
#include <strstream>
#include <string>
#include <map>

// ### Function Definitions ###

static QueryHeader qinfo;

//map<SequenceNumber,SequenceNumber> pairs;

struct ReadInfo
{
  SequenceNumber pairNum;
  int            insertMean;
  int            insertSD;
  // int clipStart;
  // int clipEnd;
  ReadInfo( SequenceNumber p, int m, int s ) :
    pairNum(p), insertMean(m), insertSD(s) {}
  ReadInfo() { cout << "Default read info ctor call" << endl; }
};

typedef map<SequenceNumber,ReadInfo> ReadInfoType;
ReadInfoType readInfo;


map<SequenceNumber, int> clipStarts;
map<SequenceNumber, int> clipEnds;
vector<int> contigSizes;

// dodgyMatches stores the read and contig number of any
// matches that are anomalous for whatever reason.
// first is for read number, second is for contig number
typedef multimap<SequenceNumber,SequenceNumber> DodgyMatchesType;
DodgyMatchesType dodgyMatches;

// When both halves of a pair hit to the same contig it is possible 
// to compute their insert size. If so, it is placed in here, indexed
// by whichever read has the *lower* sequence number
typedef map<SequenceNumber,int> InsertSizesType;
InsertSizesType insertSizes;

// First thing received from the server is a list of sequence numbers
// and names for all reads involved in maches with the query contigs.
// They are placed in here

typedef map<SequenceNumber,std::string> NamesType;
NamesType names;

struct ReadHitOnContig
{
  ReadHitOnContig( int p, bool f ) : position(p), isForward(f) {}
  int position;
  bool isForward;
};

// This stores information for all reads hits on a given contig, indexed
// by the *read's* sequence number. Should only be one full length
// hit for each read on a given contig (if more than 1, the read is added
// to dodgyMatches), therefore only need a map not a multimap. 
typedef map<SequenceNumber,ReadHitOnContig> ContigHitsType;

struct ContigHitOnRead
{
  ContigHitOnRead ( SequenceNumber c, int p, bool f ) :
    contigNum(c), position(p), isForward(f) {}
  SequenceNumber contigNum;
  int position;
  bool isForward;
};

// Once a read hit on a contig is screened by checkContig an entry for
// it is placed in here. readHits is accumulated over all contigs in the query
// Indexing is by the sequence number of the *read*.
typedef multimap<SequenceNumber, ContigHitOnRead> ReadHitsType;
ReadHitsType readHits;

struct ContigJoin
{
  enum { notDetermined=2000000000 };
  // This constructor ensures that a.contigNum 
  // is always less then b.contigNum
  ContigJoin
  ( ContigHitOnRead* pa, 
    ContigHitOnRead* pb,
    ReadInfo* pReadInfo) 
  {
    if (pa->contigNum<pb->contigNum)
    { 
      a=pa; b=pb;
    }
    else
    { 
      a=pb; b=pa;
    }
    
    int minInsert 
      =  ((a->isForward) 
      ? contigSizes[a->contigNum]-a->position // TBD +1 ?????
      : a->position )
      + ((b->isForward) 
      ? contigSizes[b->contigNum]-b->position // TBD +1 ?????
      : b->position);

    if (minInsert>pReadInfo->insertMean+(3*pReadInfo->insertSD))
    { 
      cout << "WARNING: lower bound for insert=" << minInsert
           << " mean=" << pReadInfo->insertMean
           << " SD=" << pReadInfo->insertSD << endl;
      // TBD throw an exception, catch outside ctor
    }

    gapSize=pReadInfo->insertMean-minInsert;
    gapSD=pReadInfo->insertSD;

    if (gapSize<0)
    {
      cout << "WARNING: estimated gap size=" << gapSize
           << endl;
      //      gapSize=0;
    }
    else cout << "Computed gap size of " << gapSize 
              << ", SD=" << gapSD << endl;


  }

  ContigHitOnRead* a;
  ContigHitOnRead* b;
  int gapSize;
  int gapSD;

  bool operator<( const ContigJoin& rhs ) const
  {
    return (    ( a->contigNum < rhs.a->contigNum )
             || (    ( a->contigNum == rhs.a->contigNum )
                  && ( b->contigNum < rhs.b->contigNum )) ); 
  }

};

typedef vector<ContigJoin> ContigJoinType;

// readHits is used to generate a set of contig joins. e.g. if the forward
// half of a read hits contigs 4, 8, 11 and the reverse hits contigs 5 & 10,
// contig joins will be added for 
// (4,5), (4,10), (5,8), (5,11), (8,10) and (10,11)
ContigJoinType contigJoins;

// matches must extend to within maxIn bases of clip points
// to be accepted
int maxIn=50;

// if a computed insert size (or lower bound) exceeds this
// something is deemed to be afoot...
int maxInsertSize=8000;

// if two computations of the insert size for a read differ
// by more than this, something is deemed to be afoot...
int maxInsertDiff=50;

// gets rid of all the extra stuff in the read name
string truncateReadName( const string& name )
{
  string::size_type nameEnd = name.find_first_of('\t');
  return (( nameEnd == string::npos ) ? name : name.substr(0,nameEnd));

}


// checkContig: contigHits contains all reads that hit on a contig
// for which a) the match lasts the full length of the read.
// and b) there is a mate pair for the read. The hits are checked, any
// good ones are added to readHits, then contigHits is cleared ready for
// the next contig
void checkContig( ContigHitsType& contigHits, SequenceNumber lastQueryNum )
{
            cout << "\n[ All hits received for contig, checking them...\n";

            for( ContigHitsType::iterator i(contigHits.begin());
                 i!=contigHits.end();
                 ++i )
            {
              cout << i->first << " " << truncateReadName(names[i->first])
                   << " " << i->second.position
                   << " " << (i->second.isForward?'F':'R');

              ContigHitsType::iterator myPair 
                = contigHits.find(readInfo[i->first].pairNum);

              if (myPair==contigHits.end())
              {
                int minInsert = (i->second.isForward) 
                ? contigSizes[lastQueryNum] - i->second.position
                : i->second.position;
                // lower bound on insert size
                cout << " no pair, lower bd on insert=" << minInsert;
                if (minInsert>maxInsertSize)
                {
                    cout << " REJECTED! - lower bound > max insert size\n";
                    dodgyMatches.insert
                      (make_pair(i->first,lastQueryNum));
                    continue; // next entry, don't add to readHits
                }
              }
              else
              {
                cout << " pair found"; 
                if (i->first<myPair->first)
                  // only process for the half of the pair that has
                  // lowest subject num 
                {
                  int insertSize
                    = abs(i->second.position-myPair->second.position);
                  cout << ", insert size=" << insertSize;
                  if (insertSize>maxInsertSize)
                  {
                    cout << " REJECTED! - insert size too large\n";
                    dodgyMatches.insert
                      (make_pair(i->first,lastQueryNum));
                    dodgyMatches.insert
                      (make_pair(myPair->first,lastQueryNum));
                    continue; // next entry, don't add to readHits

                  }
                  else if (i->second.isForward==myPair->second.isForward)
                  {
                    cout << " REJECTED! - both hits in same direction\n";
                    dodgyMatches.insert
                      (make_pair(i->first,lastQueryNum));
                    dodgyMatches.insert
                      (make_pair(myPair->first,lastQueryNum));
                    continue; // next entry, don't add to readHits
                  }
                  else
                  {
                    // check insert size does not contradict a previous one
                    if (insertSizes.find(i->first)!=insertSizes.end())
                    {
                      if (   abs(insertSize-insertSizes[i->first])
                           > maxInsertDiff)
                      {
                        cout << " REJECTED! - inconsistent insert sizes ("
                             << insertSize << " and " 
                             << insertSizes[i->first]
                             << ")\n";
                        dodgyMatches.insert
                          (make_pair(i->first,lastQueryNum));
                        dodgyMatches.insert
                          (make_pair(myPair->first,lastQueryNum));
                        continue; // next entry, don't add to readHits
                      }
                      cout << " confirmed existing insert size value\n";
                    }
                    else insertSizes[i->first]=insertSize;


                  } // ~else
                } // ~if (i->first...

              } // ~else

              // if the read passes all this the hit gets added to readHits.
              // It can still be discarded if the read fails a check elsewhere
              // and gets added to dodgyMatches
              readHits.insert
              ( make_pair
                ( i->first,ContigHitOnRead
                  ( lastQueryNum, i->second.position, i->second.isForward )));

              cout << endl;
            } // ~for
            cout << "...finished checking hits for contig ]" << endl;

            contigHits.clear();


} // ~checkContig




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

  typedef pair< WordSequence, std::string> QueryInfo;

  char  recvline[MAXLINE];

  WordSequence seq;

  // get handshake from server

  SocketInterface socket( sockfd);
  
  Handshake hello;

  socket.receiveStruct(&hello);

  socket.checkSocketEmpty();

  cerr << "Server is using hash table with " << hello.bitsPerSymbol 
       << " bits per symbol, " << hello.wordLength 
       << " symbols per word.\n";

  if ((hello.bitsPerSymbol==gBaseBits)&&(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
  cerr << "Server will reject queries of more than " 
       << hello.maxBufferSize << " words in total.\n";

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

  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;
  }

  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

  //  int totalNumWords(0);
  const QueryInfo dummy;
  vector<QueryInfo> query; 
  query.push_back(dummy);
  qinfo.numQueryWords = 0;
  
  contigSizes.push_back();
  
  while ( seqReader.getNextSequence
          ( query.back().first, hello.wordLength ) != -1 )
  {
    seqReader.getLastSequenceName( query.back().second );
    qinfo.numQueryWords += query.back().first.size();

    contigSizes.push_back
    ( 
        ( query.back().first.size()-1) * hello.wordLength 
        + query.back().first.getNumBasesInLast()
    );
    query.push_back(dummy);
  }
  query.pop_back();

  //  QueryHeader qinfo;

  qinfo.numQuerySeqs = query.size();

  socket.sendStruct(&qinfo);

  for( vector<QueryInfo>::iterator i(query.begin()) ; 
       i != query.end() ; ++i )
  {
    socket.sendSequence(i->first);
  }


  MatchHeader response;
  socket.receiveStruct(&response);

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

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

  pair<SequenceNumber,std::string> p;

  for ( int i(0) ;  i < response.numSubjectNames ; i++ )
  {
       socket.receiveStruct(&p.first);
       socket.receiveString(p.second);
       names.insert(p);
       //       cout << i << ": " << p.first << " " << p.second << endl;
  } // ~for

  // Here we parse the pair file to get pair info
  // and parse the read names to get clipInfo

  int  pairBufSize(200);
  char pairBuf[pairBufSize];
  int numPairs(0);

  // lastWanted is the largest sequence number which appears in the match,
  // ie for which we are interested in pairing information

  SequenceNumber p1(0), p2(0), firstWanted(0), lastWanted(0);

  if (!names.empty()) 
  {
    NamesType::iterator i(names.begin());
    firstWanted=i->first;
    i=names.end();
    lastWanted=(--i)->first;
    cout << "Range of interest: " << firstWanted 
         << " to " << lastWanted << endl;
  }

  int insertMean, insertSD;

  cout << "Reading the rest of the pair file ..." << endl;

  while (pairFile.getline(pairBuf, pairBufSize, '\n'))
  {
    sscanf(pairBuf,"%d %d",&p1, &p2);
    // p1 always less than p2, so scan until p2>=firstWanted
    if (p2 >= firstWanted) break;
  }

  do
  {
    sscanf(pairBuf,"%d %d %d %d",&p1, &p2, &insertMean, &insertSD);
    // pairs are in order of p1, so finish once p1 gets past lastWanted
    if (p1 > lastWanted) break;
    if (names.find(p1)!=names.end()) 
    { // pairs[p1]=p2; 
      readInfo.insert(make_pair(p1,ReadInfo(p2,insertMean,insertSD)));
      numPairs++; 
    }
    if (names.find(p2)!=names.end()) 
    { // pairs[p2]=p1; 
      readInfo.insert(make_pair(p2,ReadInfo(p1,insertMean,insertSD)));
      numPairs++; 
    }
  }
  while (pairFile.getline(pairBuf, pairBufSize, '\n'));

  pairFile.close();

  cout << "Found pairing information for " << numPairs << " reads\n";

  string::size_type tabPos;

  // Now grab the clip info from the names

  for ( NamesType::iterator i(names.begin());
        i!=names.end(); ++i )
  {
    tabPos=i->second.find('\t');
    ++tabPos=i->second.find('\t',tabPos);
    sscanf
      (&(i->second.c_str()[tabPos]),
       "\tcl|%d\tcr|%d",
       &clipStarts[i->first],
       &clipEnds[i->first] );

  } // ~FOR


  // Now read in the match info from the socket.

  MatchInfo thisMatch;

  //  cout << setprecision(2) << setiosflags(ios::fixed);

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


  int numAccepted(0);
  SequenceNumber lastQueryNum(0);

  ContigHitsType contigHits;

  for ( int i(0); i < response.numMatches ; i++ )
  {

    socket.receiveStruct(&thisMatch);
    cout 
      << query[thisMatch.queryNum-1].second << "\t"
      //query[thisMatch.queryNum-1].second << "\t"
      << thisMatch.queryStart << "\t"
      << thisMatch.queryEnd << "\t"
      << thisMatch.subjectNum << "\t"
      << truncateReadName(names[thisMatch.subjectNum]) << "\t"
      << thisMatch.subjectStart << "\t"
      << thisMatch.subjectEnd << "\t"
      << ((thisMatch.isForward)?"F":"R") << "\t"
      << thisMatch.numBases << " "; 
    //     << "\t"
    //     << 100.0*thisMatch.  numBases/
    //       (thisMatch.queryEnd-thisMatch.queryStart+1);

    // Only use hits on reads for which there is a read pair.
    // Hits must extend to within maxIn pairs of clip points.

    //    if (    (pairs[thisMatch.subjectNum]!=0) // half of a pair
    if (    (readInfo.find(thisMatch.subjectNum)!=readInfo.end()) 
         && ( thisMatch.subjectStart 
              <   clipStarts[ thisMatch.subjectNum ] + maxIn )  
         && ( thisMatch.subjectEnd   
              >   clipEnds[   thisMatch.subjectNum ] - maxIn )  )
    {

      if ((thisMatch.queryNum!=lastQueryNum)&&(lastQueryNum!=0)) 
      {
        // Then we have moved onto the hits for the next contig.
        // Check the hits for the previous one
        checkContig(contigHits, lastQueryNum);
      } // ~if (thisMatch ..
      lastQueryNum=thisMatch.queryNum;

      // place details of this hit into contigHits
      if (!contigHits.insert
          ( make_pair
            ( thisMatch.subjectNum,ReadHitOnContig
              ( thisMatch.isForward 
                ?   thisMatch.queryStart
                - thisMatch.subjectStart
                + 1
                :   contigSizes[thisMatch.queryNum]
                - thisMatch.queryStart
                + thisMatch.subjectStart, 
                thisMatch.isForward
              )
              )
            ).second) 
      {
        // If this insert fails then the same read has a full length hit
        // twice in the same contig. It is therefore dodgy.
        cout << "REJECTED! - same read hits twice in same contig! "
             << query[thisMatch.queryNum-1].second << "0"
             << thisMatch.subjectNum << " "
             << truncateReadName(names[thisMatch.subjectNum]) << "\t";
        dodgyMatches.insert
          (make_pair(thisMatch.subjectNum,thisMatch.queryNum));


      } // ~if (!contigHits ....
      else
      {
        // Accepted means `accepted for further checking.' Hit can
        // still be discarded at the checkContigs stage

        cout << "accepted for checking, pair=" 
             << readInfo[thisMatch.subjectNum].pairNum;
        numAccepted++;
      }

    } // ~if

    cout << endl;
    

  } // ~for

  // Check the hits for the last contig
  if (contigHits.size()!=0) checkContig(contigHits, lastQueryNum);

  cout << numAccepted << " reads passed the initial checking stage" << endl;

  socket.checkSocketEmpty();

  if (dodgyMatches.size()!=0)
  {
    cout << "Found " << dodgyMatches.size() << " potential inconsistencies:\n";
    for (DodgyMatchesType::iterator i(dodgyMatches.begin());
         i!=dodgyMatches.end();++i)
    {
      cout << "read: " << i->first << " " << truncateReadName(names[i->first])
           << " contig:" << i->second << " " << query[i->second-1].second 
           << endl;
      // don't attempt to deduce any O&O from these
      readHits.erase(i->first);
    }
    cout 
    << "These reads will not be used for the order/orientation calculation\n";
  }

  cout << "Found " << readHits.size() << " potential join points:\n";

  ReadHitsType::iterator ub,i(readHits.begin());
  //  SequenceNumber pairNum;
  ReadInfo* pReadInfo;

  
  //  for (ReadHitsType::iterator i(readHits.begin());i!=readHits.end();++i)
  while(i!=readHits.end())
  {
    cout << " --- \n";
    ub = readHits.upper_bound(i->first);
    pReadInfo = &(readInfo[i->first]);
    //    pairNum = pairs[i->first];
    //    InsertSizesType::iterator myInsert;

    if ((pReadInfo->pairNum!=0)&&(i->first<pReadInfo->pairNum))
    {
      cout << "contigs hit by forward read: " <<  i->first << " " 
           << truncateReadName(names[i->first]) << endl; 

      //      myInsert = insertSizes.find(i->first);
      //      if (myInsert!=insertSizes.end())
      //      {
      //        cout << "Got an insert size of " << myInsert->second 
      //             << " for this read pair" << endl;
      //      }
      //      else cout << "No insert size for this read pair" << endl;

      for (ReadHitsType::iterator j(i);j!=ub;j++)
      {      
        cout << j->second.contigNum << " " 
             << query[j->second.contigNum-1].second << " "
             << j->second.position << " " 
             << (j->second.isForward?'F':'R') << endl;
        for (ReadHitsType::iterator 
               k(readHits.lower_bound(pReadInfo->pairNum));
             k!=readHits.upper_bound(pReadInfo->pairNum);k++)
        {      
          if (j->second.contigNum!=k->second.contigNum)
          {
            contigJoins.push_back(ContigJoin(&j->second,&k->second,pReadInfo));
            //      if (myInsert==insertSizes.end())
            //        contigJoins.push_back(ContigJoin(j->second,k->second));
            //      else
            //     contigJoins.push_back(ContigJoin(j->second,k->second, 
            //    myInsert->second));

          }
        } // ~for k

      } // ~for j
      cout << "contigs hit by reverse read: " 
           << pReadInfo->pairNum << " "
           << truncateReadName(names[pReadInfo->pairNum]) << endl; 
      for (ReadHitsType::iterator 
             j(readHits.lower_bound(pReadInfo->pairNum));
           j!=readHits.upper_bound(pReadInfo->pairNum);j++)
      {      
        cout << j->second.contigNum << " " 
             << query[j->second.contigNum-1].second << " "
             << j->second.position << " " 
             << (j->second.isForward?'F':'R') << endl;

      } // ~for j
    } // ~if

    i=readHits.upper_bound(i->first);



  } // ~while


    sort(contigJoins.begin(),contigJoins.end());

    cout << "Found " << contigJoins.size() 
         << " potential joins between contigs\n";

    for (ContigJoinType::iterator i(contigJoins.begin());
         i!=contigJoins.end();++i)
    { 
      cout << "Contig " << i->a->contigNum << " " 
           << query[i->a->contigNum-1].second 
           << i->a->position << "/"
           << contigSizes[i->a->contigNum] << " "
           << (i->a->isForward?'F':'R')
           << " joins contig "
           << i->b->contigNum << " " << query[i->b->contigNum-1].second 
           << i->b->position << "/" 
           << contigSizes[i->b->contigNum] << " "
           << (i->b->isForward?'F':'R') << " " 
           << i->gapSize << " "
           << i->gapSD << endl;

      //      if (i->distance!=ContigJoin::notDetermined)
      //      {
      //        cout << " distance=" << i->distance;
      //`      } // ~if
      cout << endl;

    } // ~for





} // ~sendQuery








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;
    //    buf << ">unnamedQuery\n" << cin.rdbuf();
    //  cout << a.rdbuf();
    istream istr(buf.rdbuf());



    if( numArgs!=4)
    {
      cerr 
<< "syntax: " << args[0] 
<< " serverMachine serverPort pairFilePrefixName\n"
<<"serverMachine: name of machine on which server is running\n"
<<"serverPort   : port number on that machine\n"
<<"pairFileName : names of all valid read pairs\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=200;//atoi(args[3]);

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


    // read in ASCII file containing pairing and clip information


    string pairFileName((string)args[3]+(string)".pair");
    ifstream pairFile(pairFileName.c_str());

    if (pairFile.fail())
    {
      cout << "Could not open pairs file " << args[3] << ".pair\n";
      throw SSAHAException("Could not open pair file");
    }  
  
    cout << "Reading from " << args[3] << ".pair ..." << endl;

    string buff;
    
    // print out any comments, leave the rest til later
    while (pairFile.peek()=='#')
    {
      getline(pairFile,buff,'\n');
      cout << buff << endl;
    }

//      int estimatedNumReads(21000000);

//      pairs.reserve(estimatedNumReads);
//      clipStarts.reserve(estimatedNumReads);
//      clipEnds.reserve(estimatedNumReads);

//      pairs.push_back();
//      clipStarts.push_back();
//      clipEnds.push_back();

//      while (pairFile.peek()!=EOF)
//      {
//        pairs.push_back();
//        pairFile >> pairs.back();
//        clipStarts.push_back();
//        pairFile >> clipStarts.back();
//        clipEnds.push_back();
//        pairFile >> clipEnds.back();
//      }
   
//      pairs.pop_back();
//      clipStarts.pop_back();
//      clipEnds.pop_back();

//      assert(pairs.size()==clipStarts.size());
//      assert(pairs.size()==clipEnds.size());

//      cout << " ... done\n";

//      cout << "According to the pairs file there are "
//       << pairs.size()-1 << " reads in the server hash table" << endl;

    // set up rest of query parameters (same as for SSAHAClient)

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

    string queryType("DNA");
    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");
    } // ~if

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

    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";
      pReader = new SequenceReaderFile(cin, '>', '>', pEncoder, cerr );
    } // ~if   
    else
    {
      cerr << "Assuming input text is a plain string of " 
           << queryType << " data.\n";

      //      ostrstream buf;
      buf << ">unnamedQuery\n" << cin.rdbuf();
      //  cout << a.rdbuf();
      istream istr(buf.rdbuf());
      pReader = new SequenceReaderFile(istr, '>', '>', pEncoder, cerr );

      //      cin >> source;
      //if(queryType=="DNA")pReader = new SequenceReaderStringDNA(source,cerr);
      //     else pReader = new SequenceReaderStringProtein(source,cerr);
                                                                    
    } // ~else

    //  bzero(&servaddr, sizeof(servaddr));
    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, pairFile);

    exit(0);

  }    
  catch (const NetworkException& err )
  {
    cerr << "Caught NetworkException: " << err.what() << "\n";
    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);
  }



}

// Name:
// Arguments:
// TYPE  NAME  IN/OUT COMMENT
// Returns: TYPE COMMENT

// End of file OrderOrient.cpp

---