Hi,
I'm trying to compile a program (forgeG assembler) but I stuck with the following errors:
makeHash.cc:230 :21: error: conflicting declaration ‘std::vector<Cl ipInfo2>& clipInfo’
makeHash.cc:229 :15: error: ‘clipInfo’ has a previous declaration as ‘ClipInfoMap& clipInfo’
makeHash.cc:245 :21: error: conflicting declaration ‘std::vector<Cl ipInfo2>& clipInfo’
makeHash.cc:244 :15: error: ‘clipInfo’ has a previous declaration as ‘ClipInfoMap& clipInfo’
makeHash.cc: In function ‘int main(int, char**)’:
makeHash.cc:630 :75: error: too many arguments to function ‘void processFileStag e1(std::string& , FILE*, std::ofstream&, uint32, bool, uint32, ClipInfoMap&)’
makeHash.cc:219 :6: note: declared here
makeHash.cc:654 :61: error: too many arguments to function ‘void processFileStag e2(FILE*, uint32, bool, uint32, ClipInfoMap&)’
makeHash.cc:236 :6: note: declared here
makeHash.cc:663 :39: error: too many arguments to function ‘void processFileStag e2(FILE*, uint32, bool, uint32, ClipInfoMap&)’
makeHash.cc:236 :6: note: declared here
makeHash.cc:532 :10: warning: unused variable ‘infEst1’ [-Wunused-variable]
make: *** [makeHash.o] Error 1
can someone help me...I'm not c/c++ programmer.. :(
Thanx
P
I'm trying to compile a program (forgeG assembler) but I stuck with the following errors:
makeHash.cc:230 :21: error: conflicting declaration ‘std::vector<Cl ipInfo2>& clipInfo’
makeHash.cc:229 :15: error: ‘clipInfo’ has a previous declaration as ‘ClipInfoMap& clipInfo’
makeHash.cc:245 :21: error: conflicting declaration ‘std::vector<Cl ipInfo2>& clipInfo’
makeHash.cc:244 :15: error: ‘clipInfo’ has a previous declaration as ‘ClipInfoMap& clipInfo’
makeHash.cc: In function ‘int main(int, char**)’:
makeHash.cc:630 :75: error: too many arguments to function ‘void processFileStag e1(std::string& , FILE*, std::ofstream&, uint32, bool, uint32, ClipInfoMap&)’
makeHash.cc:219 :6: note: declared here
makeHash.cc:654 :61: error: too many arguments to function ‘void processFileStag e2(FILE*, uint32, bool, uint32, ClipInfoMap&)’
makeHash.cc:236 :6: note: declared here
makeHash.cc:663 :39: error: too many arguments to function ‘void processFileStag e2(FILE*, uint32, bool, uint32, ClipInfoMap&)’
makeHash.cc:236 :6: note: declared here
makeHash.cc:532 :10: warning: unused variable ‘infEst1’ [-Wunused-variable]
make: *** [makeHash.o] Error 1
Code:
#include <cstdlib>
using namespace std;
#include <math.h>
#include <string>
#include <iostream>
#include <fstream>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <map>
#include <signal.h>
#include <sys/mman.h>
#include <sys/types.h>
#include <fcntl.h>
#include <sys/stat.h>
#include "forgeG.h"
#include "md5.h"
//extern "C" {
#include "mpi.h"
//}
// #define DEBUG XXXX
// #define DEBUG_READ XXX
int32 debug1 = 128888; // 231714;
//
// Switch on to support contamination screening
//
#define USE_CONTAM_CODE XXX
bool debugContam = false;
//
// Get a record of k-mers that exceed expected depth
#define DUMP_FREQSEQS XXX
//
// Definte this to store the actual sequence in the bucket value, not the hash value
#define STORE_KMER XXX
// For readFasta
extern char lastLine[];
extern char currLine[];
//
// Processor physical/logical organization
extern uint32 G_master;
extern vector<uint32> G_PXL_P2J;
extern vector<uint32> G_PXL_J2P;
uint32 bCastBufferSz = 10000000;
uint32 hSizeG = 0;
uint32 hSizeOrig;
extern uint32 genomeSize;
extern float genomeCoverage;
extern long long baseTotal;
uint32 maxHisto = 200;
uint32 repMerMinDepth = 99999999;
uint32 depthExpectMin = 99999999;
uint32 depthExpectMax = 99999999;
uint32 kMerTile = 0; // reset later to word size if not overridden
uint32 kMerContamTile = 1;
// This is set to write all k-mers into the repmer file regardless of depth
bool dumpWholeHashTable = false;
bool haveTextClip = false;
bool haveBinClip = false;
int max(int a,int b) { return a>b?a:b;}
//
// Key data structure for holding hits. We need gads of these, so
// it needs to be as small as possible.
//
struct Bucket {
// uint32 v; // Primary hash value
// uint32 v2; // Secondary from other bit of base
unsigned long long v;
uint32 count; // # times this has been seen in input
// uint32 pos;
Flag_t flag; // Used, unused, skipped etc
#if defined(USE_CONTAM_CODE)
uchar isContam; // Non zero if a contam item - stores contam entry # (numbered from 1)
#endif
// ZZZ
};
#if defined(USE_CONTAM_CODE)
string contamFileSuffix = ".screen";
//
// Structures and prototypes for contamination screening code
//
struct ContamEntry {
string name;
uint32 off;
uint32 len;
ContamEntry(string n,uint32 o,uint32 l) {
name = n;
off = o;
len = l;
}
ContamEntry(const ContamEntry &x) { *this = x; }
void operator = (const ContamEntry &x) {name = x.name; off = x.off; len = x.len;}
};
bool injectContam(
int myId,
string &fileName,
Bucket *bucket,
uint32 hSize,
uint32 hSizeG,
uint32 hBot,
uint32 hTop
);
char *loadContamFasta(
string &contamFile,
vector<ContamEntry> &entry
);
#endif
//
// Protos for functions to get a prime sized hash table..
//
int selectPrime(bool *isPrime,int preferred,int upto);
bool verbose = false;
//
// Flags for tracking what we've told the user.
//
bool dashWarning = false;
bool NWarning = false;
bool XWarning = false;
bool calcMapping = true;
bool useColorSpace = false;
volatile void usage(void)
{
cerr << endl;
cerr << "makeHash:" << endl;
cerr << " Collect hash values, and store hits to hash values" << endl;
cerr << "Usage: makeHash projRoot" << endl << endl;
cerr << "-e hash table size [auto estimated default]" << endl;
cerr << "-m skip mapping file generation (must already exist and be accurate)" << endl;
cerr << "-d minDepth maxDepth provide nominal depth range expectation to help with stats calculation" << endl;
cerr << "-F dull dump of all k-mers to proj.kmer (debugging purposes)" << endl;
cerr << "-q path specify cache path for output files" << endl;
cerr << "-v verbose" << endl;
cerr << "-T n tile k-mers at intervals of n" << endl;
cerr << "-r threshold for dumping into repmer file (norminally 2*x expectedDepth)" << endl;
cerr << "-x dump whole hash table into repmer file (debugging only) == -r 1" << endl;
cerr << "-C assume input in color space" << endl;
cerr << "-r repmer lower threshold" << endl;
MPI_Abort(MPI_COMM_WORLD,123);
}
#if defined(NEEDS_GETOPT)
extern "C" {
int getopt(...);
}
#endif
extern uint32 readCount; // From support.cc
extern uint32 estReadCount; // From support.cc
typedef map<basic_string<char>,off_t> QualPosMap;
void processFileStage1(
// Bucket *bucket,
// uint32 hSize,
// QualPosMap &qPos,
string &qTmpFName,
FILE *inf,
ofstream &outf2,
uint32 readCount,
bool expectQual,
uint32 indexOffset,
ClipInfoMap &clipInfo,
vector<ClipInfo2> &clipInfo
);
//
// Process one file into the hash table
//
void processFileStage2(
// Bucket *bucket,
// uint32 hSize,
FILE *inf,
// ofstream &outf,
uint32 readCount,
bool expectQual,
uint32 indexOffset,
ClipInfoMap &clipInfo,
vector<ClipInfo2> &clipInfo
);
long long rehashTotal = 0;
long long totalClippedBases = 0;
int rehashSample = 0;
int usedCount = 0;
int hitTotal = 0;
int hitCum = 0;
int hitSeqSample = 0;
uint32 hashEventsG = 0;
uint32 coincidentHitsG = 0;
int hitSeqTotal = 0;
int hitSample = 0;
bool fullKMerDump = false;
void worker(string &projRoot,uint32 nProc,uint32 myId,string &outCachePath);
ofstream *openMapping(bool calcMapping,string file)
{
if (calcMapping) {
return new ofstream(file.c_str());
}
else {
return new ofstream("/dev/null");
}
}
int main(int argc,char **argv)
{
int nProc;
int myId;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&nProc);
MPI_Comm_rank(MPI_COMM_WORLD,&myId);
cout << "resetting signal handler" << endl;
signal(10,SIG_DFL);
signal(11,SIG_DFL);
#if defined(USE_SPLIT_HKEY)
// Get this out of the way otherwise it will all end in tears ;)
// we're going to assume that (sizeof(key) - 3) is a multiple of 4
//
// e.g 15 / 2 -> 7 15/4 -> 3 15-15/4 -1 = 11
//
// 0123456789ABCDE
// 00
// * 11
// * 10
// * 01
//
form2("%-5d: wordSize = %d\n",myId,wordSize);
cout.flush();
assert((wordSize - 3) % 4 == 0); // Make sure split key will actually work
#endif
extern char *optarg;
extern int optind;
string outCachePath=".";
testMD5Assumptions();
uint32 i,j;
time_t start;
time_t end;
time(&start);
char c;
if (argc == 1) usage();
while ((c = getopt(argc, argv, "t:he:mvFT:q:xCr:d:")) != EOF) {
switch (c) {
case 'e': // Set hash table size
hSizeG = atoi(optarg);
break;
case 'm': // No mapping
calcMapping = false;
break;
case 'F': //
fullKMerDump = true;
break;
case 'T': //
kMerTile = atoi(optarg);
break;
case 'd': //
depthExpectMin = atoi(optarg);
depthExpectMax = atoi(argv[optind++]);
if (myId==0) cout << "main : setting depthExpectMin=" <<
depthExpectMin << " depthExpectMax=" <<
depthExpectMax << endl;
break;
case 'r': //
repMerMinDepth = atoi(optarg);
if (myId==0) cout << "main : repMerMinDepth=" << repMerMinDepth << endl;
break;
case 'v': //
verbose = true;
break;
case 'x': // Dump whole hash table
dumpWholeHashTable = true;
break;
case 'h': // Show help
if (myId == 0) usage(); else return 1;
break;
case 'C': // ABI color space
useColorSpace = true;
if (myId == 0) cout << "main : Using color space " << endl;
break;
case 'q':
outCachePath = optarg;
if (myId == 0) cout << "main : Using output cachePath "
<< outCachePath << endl;
break;
case '?':
if (myId == 0) {
cerr << "ERROR: unknown option" << argv[optind] << endl;
usage();
}
return 0;
}
}
if (optind >=argc) {
cerr << "ERROR: expected projfile argument" << endl;
if (myId == 0) usage(); else return 1;
}
string projRoot = argv[optind];
// form2("nproc=%d myid=%d\n",nProc,myId);
if (depthExpectMax != 99999999) {
maxHisto = max(10,2 * depthExpectMax);
if (myId == 0) cout << "main : setting maxHisto to " << maxHisto << endl;
}
if (kMerTile ==0) {
#if defined(USE_SLXA)
kMerTile = 1;
#else
kMerTile = wordSize;
#endif
}
// ************************************************************************
//
// Setup worker types
//
vector<ResJob> job;
ResJob jb;
// Add master job
jb.jobType = RJ_MASTER; jb.id = 0; job.push_back(jb);
// Add mem jobs
for(i=1;i<(uint32)nProc;i++) {
jb.id = i;
jb.jobType = RJ_CPU;
jb.needsFile = "";
job.push_back(jb);
}
// Now do processor allocation
if (myId == 0) {
cout << "main : start processor allocation as temp master" << endl;
procDiscoverMaster(nProc,job);
} else {
form1("%-5d: start processor allocation slave\n",myId);
procDiscoverSlave(myId,nProc);
}
// ************************************************************************
if (myId == (int32)G_master) {
string report = projRoot + ".hashReport";
string error;
if (!unlinkAndTest(report,error)) {
cerr << "main : could not delete " << report <<
" check permissions: " << error << endl;
MPI_Abort(MPI_COMM_WORLD,393);
}
//
// Load in dynamically calculated statistics for this assembly.
//
readDyn(projRoot + ".dynstats");
form1("Using readCount/nodeCount = %d\n",readCount);
form1("Suggested hash size = %d\n",hSizeG);
form1("coverage = %f\n",genomeCoverage);
form1("baseTotal = %lld\n",baseTotal);
form1("kMerTile = %d\n",kMerTile);
if (dumpWholeHashTable) {
cout << "main : Dumping whole hash table to repmer file(s) " << endl;
}
#if defined(USE_CONTAM_CODE)
string contamFile = projRoot + contamFileSuffix;
#endif
if (hSizeG == 0) {
// If hash table suggested size is 0, we work it out ourselves
uint32 baseTotalTmp = baseTotal;
#if defined(USE_CONTAM_CODE)
struct stat st;
// Allow for possible addition of contamination bases to hash table
if (stat(contamFile.c_str(),&st) != -1) {
cout << "main : Expanding hash table required to allow for contaminants by "
<< st.st_size << " bases " << endl;
baseTotalTmp += st.st_size * kMerTile / kMerContamTile; // We are going to step these every 1 bp
} else {
cout << "main : no contam file present" << endl;
}
#endif
hSizeG = (int)(baseTotalTmp / (float)kMerTile);
#if defined(USE_SPLIT_HKEY)
hSizeG *= 4; // multiple split keys
#endif
cerr << "main : setting hSizeG automatically to " << hSizeG << endl;
}
hSizeOrig = hSizeG;
}
if (myId != (int32)G_master) {
worker(projRoot,nProc,myId,outCachePath);
} else {
#if defined(USE_CONTAM_CODE)
//
// Load the contam file, only to get the names of the contaminants for later
//
vector<ContamEntry> contamEntry;
string contamFile = projRoot + contamFileSuffix;
// Load entries
struct stat statFile;
if (stat(contamFile.c_str(),&statFile) != -1) {
char *genome = loadContamFasta(contamFile,contamEntry);
delete [] genome;
}
#endif
//
// Two handles for the main genomic fasta (one for each pass) +
// one handle for the quality data
//
FILE *inf = fopen((projRoot + ".fasta").c_str(),"r"); // ios::binary | ios::in);
FILE *inf2 = fopen((projRoot + ".fasta").c_str(),"r"); // ios::binary | ios::in);
FILE *inf3 = fopen((projRoot + ".qual").c_str(),"r"); // ios::binary | ios::in);
cout << "main: allocating three buffers for file reading\n" ;
char *buffer1 = new char[10000000];
char *buffer2 = new char[10000000];
char *buffer3 = new char[10000000];
setvbuf(inf,buffer1,_IOFBF,10000000);
setvbuf(inf2,buffer2,_IOFBF,10000000);
setvbuf(inf3,buffer3,_IOFBF,10000000);
//
// Open a fasta file for EST data - it may not exist
//
FILE * infEst1 = fopen((projRoot + ".est.fasta").c_str(),"r"); // ,ios::binary | ios::in);
FILE * infEst2 = fopen((projRoot + ".est.fasta").c_str(),"r"); // ios::binary | ios::in);
//
// Make hash table
//
int upto = (int)(hSizeG * 1.6);
//
// Make all primes up to our desired size
//
bool *isPrime = genPrimes(upto);
//
// Choose a prime just before
//
hSizeG = selectPrime(isPrime,(int)(1.59 * hSizeG),upto);
delete [] isPrime;
form1("main : Actual hash size = %d buckets\n",hSizeG);
cout << "main : sending hash size" << endl; cout.flush();
MPI_Bcast(&hSizeG,1,MPI_UNSIGNED,G_master,MPI_COMM_WORLD);
cout << "main : sent hash size" << endl; cout.flush();
MPI_Bcast(&genomeCoverage,1,MPI_FLOAT,G_master,MPI_COMM_WORLD);
uint32 readCountTotal = readCount + estReadCount;
cout << "sending rc " << endl; cout.flush();
MPI_Bcast(&readCountTotal,1,MPI_UNSIGNED,G_master,MPI_COMM_WORLD);
cout << "sent rc " << endl; cout.flush();
basic_string<char> line;
basic_string<char> name;
vector<uchar> quals;
off_t startPos;
//
// Load clip info if any
//
ClipInfoMap clipInfo;
vector<ClipInfo2> clipInfoB; // if/Binary clip info
string clipFileB = projRoot + ".clipb";
cout << "main : checking for " << clipFileB << endl;
if (stat(clipFileB.c_str(),&statFile) != -1) {
cout << "main : Loading " << clipFileB << endl;
if (!readBinaryClip(clipFileB,clipInfoB)) MPI_Abort(MPI_COMM_WORLD,102);
haveBinClip = true;
} else {
cout << "main : Loading clip file " << (projRoot + ".clip") << " if avail" << endl;
int ret;
ret = loadClip((projRoot + ".clip").c_str(),clipInfo);
if (ret ==2 ) {
MPI_Abort(MPI_COMM_WORLD,192);
};
if (ret == 0) cout << "done loading clip file " << endl;
haveTextClip = true;
}
//
// Locate quality scores in file
//
cout << "Map qual file" << endl;
// QualPosMap qPos;
string qNameTmp = projRoot + ".tmpqmap";
ofstream qTmpF(qNameTmp.c_str());
ofstream &outf2 = *openMapping(calcMapping,projRoot + ".mapping");
if (calcMapping) {
int inc = max(readCount/10,1);
for(i=0;i<readCount;i++) {
if (!readQual(inf3,name,quals,startPos)) {
cout << "EOF reached prematurely: " << i << " of " << readCount << " qual records " << endl;
MPI_Abort(MPI_COMM_WORLD,969);
}
if (i % inc == 0) {
form1("%d%%.." , (i / inc) * 10);
cout.flush();
}
qTmpF << name << " " << startPos << endl;
// qPos[name] = startPos;
}
} else {
cout << "Warning: assuming " << projRoot + ".mapping" << " is correct" << endl;
}
qTmpF.close();
cout << endl;
// ===========================================
// Pass one, insert novel hash values
// ===========================================
cout << "Pass one: Insert novel hash values and map fasta file" << endl;
processFileStage1(qNameTmp,inf,outf2,readCount,true,0,clipInfo,clipInfoB);
/*
//
// If we saw EST reads during the initial stats survey, include
// these in the hash table of hits
//
if (estReadCount != 0) {
cout << "Pass 1b: Insert novel hash values for EST" << endl;
processFileStage1(qPos,infEst1,outf2,estReadCount,
false,readCount,clipInfo,clipInfoB);
}
*/
unlink(qNameTmp.c_str());
//
// ===========================================
// Pass two, measure hits to hash values
// ===========================================
//
cout << endl;
cout << "main : Pass two: collect hits to hash values" << endl;
processFileStage2(inf2,readCount,true,0,clipInfo,clipInfoB);
//
// If we saw EST reads during the initial stats survey, include
// these in the hash table of hits
//
if (estReadCount != 0) {
cout << "main : Pass 2b: Insert novel hash values for EST" << endl;
processFileStage2(infEst2,estReadCount,
false,readCount,clipInfo,clipInfoB);
}
cout << "main : done" << endl;
delete [] buffer1;
delete [] buffer2;
delete [] buffer3;
MPI_Status st;
#if defined(USE_CONTAM_CODE)
//
// ***************************************************************************
//
//
// Receive the contamination details
//
vector<uchar> tmpContamCount(readCount);
vector<uchar> tmpContamType(readCount);
vector<uchar> tmpContamAttempts(readCount);
vector<uint32> finalContamCount(readCount);
vector<uchar> finalContamType(readCount);
vector<uint32> finalContamAttempts(readCount);
for(i=0;i<readCount;i++) {
finalContamCount[i] = 0;
finalContamType[i] = 0;
finalContamAttempts[i] = 0;
}
cout << "main : collecting contamination info " << endl;
for(i=1;i<(uint32)nProc;i++) {
// Receive count from worker i
cout << "main : waiting on worker " << i << " (a)" << endl;
MPI_Recv(&tmpContamCount[0],readCount,MPI_CHAR,G_PXL_J2P[i],623,MPI_COMM_WORLD,&st);
cout << "main : waiting on worker " << i << " (b)" << endl;
MPI_Recv(&tmpContamAttempts[0],readCount,MPI_CHAR,G_PXL_J2P[i],624,MPI_COMM_WORLD,&st);
cout << "main : waiting on worker " << i << " (d)" << endl;
MPI_Recv(&tmpContamType[0],readCount,MPI_CHAR,G_PXL_J2P[i],625,MPI_COMM_WORLD,&st);
cout << "main : collected worker " << i << endl;
for(j=0;j<readCount;j++) {
finalContamCount[j] += (uint32)tmpContamCount[j];
finalContamAttempts[j] += (uint32)tmpContamAttempts[j];
if (tmpContamType[j]) finalContamType[j] = tmpContamType[j];
}
}
string contamResult = projRoot + ".contam";
string contamResultB = projRoot + ".contamb";
ofstream contamF(contamResult.c_str());
ofstream contamBF(contamResultB.c_str());
cout << "main : writing " << contamResult << " and " << contamResultB << endl;
if (!contamF) {
cerr << "main: ERROR, can't open file '" << contamResult << "' for writing" << endl;
cerr.flush();
MPI_Abort(MPI_COMM_WORLD,938);
}
if (!contamBF) {
cerr << "main: ERROR, can't open file '" << contamResultB << "' for writing" << endl;
cerr.flush();
MPI_Abort(MPI_COMM_WORLD,338);
}
uint32 cEntries = 0;
uint32 countContam = 0;
vector<uchar> contamDecision(readCount);
for(i=0;i<readCount;i++) {
contamDecision[i] = (finalContamCount[i]/(float)finalContamAttempts[i] > 0.4)?1:0;
if (contamDecision[i]) countContam++;
if (finalContamCount[i]) {
string name;
cout.flush();
cEntries++;
assert(finalContamType[i]-1 >=0 && finalContamType[i]-1 <(int32)contamEntry.size());
name = (finalContamType[i]< 255)?contamEntry[finalContamType[i]-1].name:"unknown";
contamF << i << "\t" << finalContamCount[i] << "\t" << finalContamAttempts[i] << "\t";
contamF.precision(4);
contamF << (float)finalContamCount[i]/finalContamAttempts[i] << "\t"
<< name << endl;
}
}
cout <<"main : closing " << contamResult << " wrote " << cEntries << " entries " << endl;
cout.flush();
contamF.close();
contamBF.write("FRGECNTM",8);
uint32 version = 0;
contamBF.write((char *)&version,4);
contamBF.write((char *)(&readCount),4);
contamBF.write((char *)(&contamDecision[0]),readCount);
contamBF.close();
// END contam writing
// ***************************************************************************
//
#endif
#if defined(DUMP_FREQSEQS)
uint32 numKmerTotal = 0;
uint32 numKmerTmp = 0;
uint32 numRepKmerTotal = 0;
uint32 numRepKmerTmp = 0;
cout << "main : getting number of rep kmers" << endl; cout.flush();
for(i=1;i<(uint32)nProc;i++) {
MPI_Recv(&numKmerTmp,1,MPI_UNSIGNED,G_PXL_J2P[i],620,MPI_COMM_WORLD,&st);
MPI_Recv(&numRepKmerTmp,1,MPI_UNSIGNED,G_PXL_J2P[i],621,MPI_COMM_WORLD,&st);
numKmerTotal += numKmerTmp;
numRepKmerTotal += numRepKmerTmp;
}
#endif
//
// Collect histogram of coverage
//
vector<uint32> histo(maxHisto);
vector<uint32> histoTotal(maxHisto);
for(i=0;i<maxHisto;i++) {
histoTotal[i] = 0;
}
// Estimate edges required for next stage
uint32 maxEdgeRequired = 0;
uint32 tmp;
cout << "main : getting required edge counts" << endl; cout.flush();
for(i=1;i<(uint32)nProc;i++) {
MPI_Recv(&tmp,1, MPI_UNSIGNED,
G_PXL_J2P[i],626, MPI_COMM_WORLD,&st);
maxEdgeRequired += tmp;
}
// sum up histogram
cout << "main : getting histogram pieces" << endl; cout.flush();
for(i=1;i<(uint32)nProc;i++) {
MPI_Recv(&histo[0],maxHisto, MPI_UNSIGNED,
G_PXL_J2P[i],627, MPI_COMM_WORLD,&st);
for(j=0;j<maxHisto;j++) {
histoTotal[j] += histo[j];
}
}
cout << "main : writing hash report" << endl; cout.flush();
//
// Write out final report
// --------------------------------------------------------------
ofstream repF((projRoot + ".hashReport").c_str());
if (!repF) { cerr << "failed to open hashReport" << endl;
MPI_Abort(MPI_COMM_WORLD,284);
}
repF << "Hash Report" << endl;
repF << "------------------------" << endl;
int32 hTotal = 0;
int32 hTotalSkipLow = 0;
uint32 modePos = 0;
uint32 mode = 0;
uint32 modeLower = (depthExpectMin != 99999999)?depthExpectMin/2:1;
for(i=1;i<maxHisto;i++) {
// Note skip first entry which really counts unused
// hash table slots
if (i>2) hTotalSkipLow += histoTotal[i];
hTotal += histoTotal[i];
if (i > modeLower && histoTotal[i] > mode) {
mode = histoTotal[i];
modePos = i;
}
}
float hashUsed = (hSizeG-histoTotal[0])/(float)hSizeG;
time(&end);
repF << "hashTableSize = " << hSizeG <<
" (~1.6 * -e param)" << endl;
repF << "-e = " << hSizeOrig << endl;
repF << "maxEdgeRequired = " << maxEdgeRequired <<
" (assumes -t set to " << (uint32)genomeCoverage*3 << ")" << endl;
repF << "kMerTile = " << kMerTile << endl;
repF << "totalHashHits = " << hTotal << endl;
repF << "Hash Table Used = " << (long long)hSizeG-(long long)histoTotal[0] << endl;
repF << "Hash Table Used = " << hashUsed * 100 << "%" << endl;
repF << "totalClippedBases = " << totalClippedBases << endl;
repF << "Implied coverage = " << totalClippedBases /
(float)genomeSize << endl;
#if defined(DUMP_FREQSEQS)
// repF << " Num kmers = " << numKmerTotal << endl;
float perc = numRepKmerTotal / (numKmerTotal + 1.0) * 100;
repF << "Num freq kmers = " << numRepKmerTotal << " (" << perc << ")%" << endl;
#endif
repF << "depthExpectMin = " << depthExpectMin << endl;
repF << "depthExpectMax = " << depthExpectMax << endl;
repF << "Mode coverage = " << modePos << endl;
repF << "ContamRemoved = " << countContam << endl;
cout << "totalClippedBases = " << totalClippedBases << endl;
cout << "Implied coverage = " << totalClippedBases /
(float)genomeSize << endl;
cout << "Mode coverage = " << modePos << endl;
cout << "Elapsed time = " << end - start << " seconds " << endl;
repF << "Elapsed time = " << end - start << " seconds " << endl;
if (hashUsed < 0.3) {
repF << endl;
repF << "Hash underused -e could be lowered to " <<
0.9 * (hSizeG-histoTotal[0]) << endl;
}
if (hashUsed > 0.8) {
repF << " Hash close to full -e should be increased to " <<
(hSizeG-histoTotal[0]) << endl;
}
int32 last;
float tail = 0;
uint32 tailC = 0;
uint32 tailMin = depthExpectMax == 99999999?0:(uint32)(depthExpectMax * 1.1);
for(last = maxHisto-1;last>=(int32)tailMin && tail < 0.01 ;last--) {
tail += (float)histo[last]/max(1,hTotalSkipLow);
tailC += histo[last];
}
repF << endl;
repF << "Coverage Count percent" << endl;
repF << "--------------------------" << endl;
for(i=1;(int32)i<=last;i++) {
repF.width(6);
repF << i;
repF.width(10);
repF << histoTotal[i];
repF.width(10);
repF.precision(2);
repF << (100 * ((float)histoTotal[i])/hTotal);
repF << endl;
}
repF << ">=";
repF.width(4);
repF << i;
repF.width(10);
repF << tailC;
repF.width(10);
repF << 100 * tail << "%" << endl;
cout << endl;
cout << endl;
cout << endl;
cout << endl;
cout << endl;
cout << " ********************************************* " << endl;
cout << " totalClippedBases = " << totalClippedBases << endl;
cout << " Implied coverage = " <<
totalClippedBases / (float)genomeSize << endl;
cout << " You might want to update " << projRoot <<
".dynstats if this is different to the pre clipped estimate "
<< endl;
cout << " Mode coverage = " << modePos << " (more reliable) "
<< endl;
#if defined(DUMP_FREQSEQS)
cout << " Num kmers = " << numKmerTotal << endl;
cout << " Num freq kmers = " << numRepKmerTotal << " (" << perc << ")%" << endl;
#endif
cout << endl;
cout << endl;
cout << " You can see " << (projRoot+".hashReport").c_str() <<
" for more information " << endl;
cout << " ********************************************* " << endl;
cout << endl;
cout << "main : at finalize" << endl;
cout.flush();
} // end if master id
// Both workers and master finished now
form1("%-5d: at finalize\n",myId); cout.flush();
MPI_Finalize();
}
//
// Process one file into the hash table
//
void processFileStage1(
string &qTmpFName,
FILE *inf,
ofstream &outf2,
uint32 readCount,
bool expectQual,
uint32 indexOffset,
ClipInfoMap &clipInfo,
vector<ClipInfo2> &clipInfoB
)
{
uint32 i;
off_t startPos=0;
basic_string<char> line;
basic_string<char> name;
lastLine[0] = 0;
currLine[0] = 0;
form1("main : allocating broadcast buffer %d bytes\n",bCastBufferSz);
char * bCastBuffer= new char [bCastBufferSz];
uint32 buffUsed = 0;
ifstream qTmp(qTmpFName.c_str());
cout << "main : starting to send reads \n";
for(i=0;i<readCount;i++) {
if (!readFasta(inf,name,line,startPos)) {
cout << "EOF3 reached prematurely: " << i << " of " << readCount << " reads" << endl;
MPI_Abort(MPI_COMM_WORLD,963);
}
if (verbose) cout << "Read=" << name << endl;
// cout << i << " " << line.size() << " " << startPos << endl;
off_t qOffset = 0;
string qName;
if (calcMapping) {
if (expectQual) {
// QualPosMap::iterator p;
// p=qPos.find(name);
qTmp >> qName >> qOffset;
if (qName != name) {
cerr << "No matching qual record for seq '" << name << "' hit qName=" << qName << endl;
MPI_Abort(MPI_COMM_WORLD,263);
}
// qOffset = p->second;
}
} // End if doing offset calculation
if (haveTextClip) {
//
// Check clipping info on this one
//
ClipInfoMap::iterator p;
p=clipInfo.find(name);
if (p != clipInfo.end() ) {
ClipInfo &c = p->second;
if (!c.good) {
line = "XX";
} else if (c.right-(int32)wordSize <= c.left) { // was 40
line = "XX";
} else {
line = line.substr(c.left+1,c.right-c.left - 1);
}
} else {
if (verbose) cerr << "no clip info '" << name << "'" << endl;
}
} else if (haveBinClip) {
if (!clipInfoB[i].good) {
line = "XX";
} else if (clipInfoB[i].right-(int32)wordSize <= clipInfoB[i].left) { // was 40
line = "XX";
} else {
line = line.substr(
clipInfoB[i].left+1,
clipInfoB[i].right-clipInfoB[i].left-1);
}
}
uint32 len = line.length();
totalClippedBases += len;
// Emit offsets to a file
if (calcMapping) {
// Write out size and offsets to each read.
outf2 << i + indexOffset << " " << line.length() << " " << startPos << " " << qOffset <<
" " << name << endl;
}
// Send buffer if it's overly full
if (buffUsed + len >= bCastBufferSz - 3) {
// Time to send buffer
bCastBuffer[buffUsed++] = 0; // terminate buffer here
bCastBuffer[buffUsed] = 0; // more to come
MPI_Bcast((void *)bCastBuffer,bCastBufferSz,MPI_CHAR,G_master,MPI_COMM_WORLD);
buffUsed = 0;
}
// Add this seq entry to buffer
strcpy(bCastBuffer+buffUsed,line.c_str());
// for(int xx=0;xx<line.size();xx++) {
// assert(line[xx] != '\n');
// }
buffUsed += len + 1;
} // End foreach sequence
//
// Finalise buffer
bCastBuffer[buffUsed] = 0; // terminate buffer
cout << "main : flushing send buffer\n";
cout.flush();
MPI_Bcast((void *)bCastBuffer,bCastBufferSz,MPI_CHAR,G_master,MPI_COMM_WORLD);
delete [] bCastBuffer;
} // End pass 1 over file
int better(char a,char b)
{
char c='X';
switch(b) {
case 'G': c = 'C'; break;
case 'A': c = 'T'; break;
case 'T': c = 'A'; break;
case 'C': c = 'G'; break;
case 'N': c = 'N'; break;
case 'X': c = 'X'; break;
default:
cerr << "ERROR: unexpected base '" << (char)b << "' in better()" << endl;
MPI_Abort(MPI_COMM_WORLD,987);
}
if (a < c) return -1;
if (a > c) return 1;
return 0;
}
inline void makeCanonical(unsigned char *canonical,unsigned char *input)
{
// Is the fwd or reverse complement the canonical form?
uint32 i,j;
bool useFwd = true;
for(i=0,j=wordSize-1;i<wordSize;i++,j--) {
//
// See which base would be smaller 5' or 3' after reversing.
// Note that we complement or not depending on the space. For
// color space it's a head to head competition - which base is
// smaller, for seq space, it's 5' versus complemented 3' base
int x = (useColorSpace)?(input[i]-input[j]) : better(input[i],input[j]);
if (x < 0) break;
if (x > 0) {
useFwd = false;
break;
}
}
if (useFwd) {
for(i=0;i<wordSize;i++) {
canonical[i] = input[i];
}
} else {
canonical[wordSize] = 0;
if (useColorSpace) {
// in color space .. Just flip the sequence around, don't complement
for(i=0,j=wordSize-1;i<wordSize;i++,j--) {
canonical[j] = input[i];
}
} else {
// Normal reverse complement
for(i=0,j=wordSize-1;i<wordSize;i++,j--) {
switch(input[i]) {
case 'G': canonical[j] = 'C'; break;
case 'A': canonical[j] = 'T'; break;
case 'T': canonical[j] = 'A'; break;
case 'C': canonical[j] = 'G'; break;
case 'N': canonical[j] = 'N'; break;
case 'X': canonical[j] = 'X'; break;
default:
cerr << "ERROR: unexpected base '" << (char)input[i] << "' in makeCanonical" << endl;
MPI_Abort(MPI_COMM_WORLD,938);
}
}
}
}
}
unsigned long long seq2Bits(const unsigned char *p,uint32 wordSize)
{
uint32 i;
long long v = 0;
for(i=0;i<wordSize;i++) {
v <<= 2;
switch(p[i]) {
case 'G': v += 0; break;
case 'A': v += 1; break;
case 'T': v += 2; break;
case 'C': v += 3; break;
case 'N': v += i&3; break;
case 'X': v += i&3; break;
#if defined(USE_SPLIT_HKEY)
case '*': v >>= 2; // split keys have gaps that we don't encode,
// so thwart attempt to move key over 2 bits.
break;
#endif
default:
cerr << "ERROR: unexpected byte '" << p[i] << "' (" << (int)p[i] << ") position " << i << endl;
throw BadSeqException();
// assert(0);
}
}
return v;
}
string bits2Seq(unsigned long long v,uint32 wordSize)
{
#if defined(USE_SPLIT_HKEY)
// Split keys need reintroduction of a gap character in the right place
// 0123456789ABCDE
// 00
// * 11
// * 10
// * 01
//
// uint32 left = wordSize/4;
//uint32 right = wordSize - 1 - left;
//uint32 mid = wordSize / 2;
int32 gapPos = -1;
switch (v & 3) {
case 0: gapPos = -1; break; // no gap
case 1: gapPos = (wordSize-1 - wordSize/4); break; // rhs
case 2: gapPos = wordSize/4; break; // lhs
case 3: gapPos = wordSize/2; break; // middle
}
v >>= 2; // get rid of leading key
#endif
string res;
uint32 i;
for(i=0;i<wordSize;i++) {
#if defined(USE_SPLIT_HKEY)
// Must reintroduce gap in split key at some stage
if ((int32)i == gapPos) {
res = "*" + res;
continue; // Skip decoding for this cycle, don't shift v
}
#endif
switch(v & 3) {
case 0: res = "G" + res; break;
case 1: res = "A" + res; break;
case 2: res = "T" + res; break;
case 3: res = "C" + res; break;
}
v >>= 2;
}
return res;
}
//
// Store an individual word into the hash table
//
inline void storeWord(
Bucket *bucket,
uint32 hSize,
unsigned char canonical[],
uint32 hSizeG,
uint32 hBot,
uint32 hTop
#if defined(USE_CONTAM_CODE)
,
uchar isContam
#endif
)
{
#if defined(USE_SPLIT_HKEY)
// Positions for skip characters
uint32 left = wordSize/4;
uint32 right = wordSize - 1 - left;
uint32 mid = wordSize / 2;
unsigned char canonicalOrig[wordSize+1];
// Loop here over different split keys
//
// we're going to assume that (sizeof(key) - 3) is a multiple of 4
//
// e.g 15 / 2 -> 7 15/4 -> 3 15-15/4 -1 = 11
//
//
// 0123456789ABCDE
// 00
// * 11
// * 10
// * 01
//
// make a specific copy of the key
// modify the sequence accordingly
uint32 sk;
for(sk=0;sk<4;sk++) {
switch(sk) {
case 0:
// Make copy, leave unchanged
strncpy((char *)canonicalOrig,(char *)canonical,wordSize+1);
break;
case 1:
// Retrieve copy and modify right
strncpy((char *)canonical,(char *)canonicalOrig,wordSize+1);
canonical[right] = '*';
break;
case 2:
// Retrieve copy and modify left
strncpy((char *)canonical,(char *)canonicalOrig,wordSize+1);
canonical[left] = '*';
break;
case 3:
// Retrieve copy and modify center
strncpy((char *)canonical,(char *)canonicalOrig,wordSize+1);
canonical[mid] = '*';
break;
default:
assert(0); // only 4 options
} // End set up key
assert(canonical[0]!=0);
// WARNING..
//!!!!! The code below is conditionally included in the loop if we are doing split
// keys, hence break in indentation..
#endif
uint32 k;
#if defined(STORE_KMER)
assert(canonical[0] != 0);
unsigned long long kMerBits = seq2Bits(canonical,wordSize);
#if defined(USE_SPLIT_HKEY)
// For split keys, we tack on 2 bits to say which type of key we
// are dealing with so we can print them out sensibly later
kMerBits <<= 2; // move key left
kMerBits |= sk; // add on 2 bit key
#endif
#endif
//
// Create a hash value from wordSize letters of DNA
//
// TODO: test djb2Hash(unsigned char *str) (above)
MD5_CTX context;
long long tmp1[2];
MD5Init (&context);
MD5Update (&context, canonical, wordSize);
MD5Final ((unsigned char *)tmp1, &context);
unsigned long long v = tmp1[0] ^ tmp1[1];
#if defined(DEBUG_KMER)
if (strncmp((char *)canonical,"GACTTCCTATTAAA",14)==0) {
cout << "XXX0: in storeword inserting fwd" << endl;
}
if (strncmp((char *)canonical,"TTTAATAGGAAGTC",14)==0) {
cout << "XXX0: in storeword inserting rev" << endl;
}
#endif
hashEventsG++;
// k = (v^v2) % hSizeG;
k = v % hSizeG;
if (k<hBot || k > hTop) return;
k -= hBot;
bool done = false;
bool found = false;
do {
rehashTotal++;
k = (k + 1) % hSize;
switch(bucket[k].flag) {
case USED:
case SKIPPED:
#if defined(STORE_KMER)
if ( bucket[k].v == kMerBits )
#else
if ( bucket[k].v == v )
#endif
{
found = true;
done = true;
coincidentHitsG++;
} else if ( bucket[k].flag == USED) {
bucket[k].flag = SKIPPED;
}
break;
case UNUSED:
done = true;
found = false;
bucket[k].flag = USED;
#if defined(STORE_KMER)
bucket[k].v = kMerBits;
#else
bucket[k].v = v;
#endif
#if defined(USE_CONTAM_CODE)
bucket[k].isContam = isContam;
#endif
usedCount++;
break;
default:
cerr << "ERROR: illegal flag '" <<
(int)(bucket[k].flag) <<"'" << endl;
exit(1);
} // End switch on bucket type
} while (!done);
rehashSample++;
hitSample++;
#if defined(USE_SPLIT_HKEY)
// Extra loop must be closed
} // End loop over possible split keys
#endif
}
void processFileStage1Par(
string &projRoot,
uint32 myId,
Bucket *bucket,
uint32 hSizeG,
uint32 hBot,
uint32 hTop,
uint32 readCount
)
{
uint32 ii=0,j,m;
uint32 hSize = hTop - hBot + 1;
uint32 lineSize = 100000;
form2("%-5d: allocating %d bytes for line\n",myId,lineSize);
char *line = new char [lineSize];
form2("%-5d: allocation %d bytes for bcastBuffer\n",myId,bCastBufferSz);
char *bCastBuffer= new char [bCastBufferSz];
#if defined(USE_CONTAM_CODE)
string contamFile = projRoot + contamFileSuffix;
struct stat st;
form2("%-5d: checking for presence of contamination file %s\n",myId,contamFile.c_str());
if (stat(contamFile.c_str(),&st) != -1) {
form2("%-5d: loading contamination file %s\n",myId,contamFile.c_str());
//
// Load a file of possible contaminants and insert it into the hash table
//
bool outcome = injectContam( myId,contamFile, bucket, hSize, hSizeG, hBot, hTop);
if (!outcome) {
cerr << endl << endl << endl << "ERROR: exitting due to contamination file load failure" << endl;
cerr.flush();
MPI_Abort(MPI_COMM_WORLD,238);
}
} else form1("%-5d: no contam file present\n",myId);
#endif
bool seqFound = false;
// uint32 v = 0;
// unsigned long long v = 0;
// uint32 v2 = 0;
verbose = true;
for(ii=0;ii<readCount;) {
try {
cout << "Read " << ii << endl;
uint32 len;
if (verbose) {
form3("%d: waiting for seq stage 1 processed=%d of %d\n",myId,ii,readCount);
cout.flush();
}
cout.flush();
MPI_Bcast(bCastBuffer,bCastBufferSz,MPI_CHAR,G_master,MPI_COMM_WORLD);
cout.flush();
if (verbose) {
form1("%d: stage 1 got buffer\n",myId);
cout.flush();
}
char c = bCastBuffer[0];
if (!(c=='G' || c=='C' || c=='T' || c == 'A' || c == 'X' || c == 'N'))
throw BadSeqException();
// assert(c=='G' || c=='C' || c=='T' || c == 'A' || c == 'X' || c == 'N');
// if (verbose) form1("%d: stage 1 got seq\n",myId);
// line[len] = 0;
for(m=0;bCastBuffer[m];m+=len + 1) {
if (ii== readCount) {
if (verbose) {
cout << "Hit end of data by count" << endl;
cout.flush();
}
break;
}
if (bCastBuffer[m] == 1) {
if (verbose) {
cout << "Hit end of data marker" << endl;
cout.flush();
}
break;
}
strcpy(line,bCastBuffer+m);
ii++; // Count sequences served // WARNING - preincremented
len = strlen(line);
// cout << "seq=" << line << " len=" << len << " " << (len%wordSize)/2 << endl;
// cout << "read seq " << len << endl;
// cout << line << endl;
// cout.flush();
if (len < wordSize) continue; // CONTINUE if too short
//
// Check hash table isn't getting overly full
//
if ((float)usedCount/hSize >=0.95) {
form1("%-5d: hash table hit 95%% exitting\n",myId);
cout.flush();
MPI_Abort(MPI_COMM_WORLD,838);
}
//
// Decide where to step from and how much to step each time.
// This might be contextual - we can take different strategies depending
// on read length.
//
bool centiRead = (len < 50); // Short reads get special treatment
const uint32 leftStart = centiRead?0:(len%wordSize)/2;
const uint32 thisStep = centiRead?1:kMerTile;
for(j=leftStart;j<=len-wordSize;j+=thisStep) {
// if (strncmp(line+j,"ATAAGATTTATAACTAATTTCAAAT",25) == 0) {
// cout << myId << " YES i=" << i << endl;
// }
//
// We want to store only canonical hash values, so we store
// the reverse complement seqeuence if it is the canonical form?
//
unsigned char canonical[wordSize+1];
makeCanonical(canonical,(unsigned char *)(line + j));
#if defined(DEBUG_KMER)
// Debugging
string tmp;
for(int zzz=0;zzz<wordSize;zzz++) tmp += line[j+zzz];
if (tmp == "GACTTCCTATTAAA" || tmp == "TTTAATAGGAAGTC") {
cout << "XXX1: i=" << i-1 << " " << tmp << endl;
}
#endif
//
// *******************************************
// Key placement of word into hash table
// *******************************************
//
assert(canonical[0] != 0);
storeWord(bucket,hSize,canonical,hSizeG,hBot,hTop
#if defined(USE_CONTAM_CODE)
,false // These are real k-mers
#endif
);
} // End look through seq
hitSeqSample++;
if (hitSeqSample/(float)hSize > 0.95) {
form3("%-5d: ERROR: hash table 95%% full %d used out of %d try",myId,hitSeqSample,hSize);
form1(" increasing -e (currently %d)\n",hSizeG);
MPI_Abort(MPI_COMM_WORLD,90);
}
if (seqFound) {
hitSeqTotal++;
}
// Periodically emit statistics
if (ii % 10000 == 0|| ii==readCount-1) {
if (hitSample > 0 && rehashSample > 0) {
form1("%-5d: ",myId);
// cout.width(5);
// cout << myId << ":";
cout <<
"pass1: " << ii <<
" usedbuckets= " << usedCount;
cout.precision(3);
cout << " " << 100*(float)usedCount/hSize << "%";
cout <<
// " avgHits= " << hitTotal / (float)hitSample <<
// " readsWithHits= " << hitSeqTotal / (float)hitSeqSample <<
" rehashTotal= " << rehashTotal << endl;
// " coincidentHits=" << coincidentHits << endl;;
hitTotal = 0;
hitSample = 0;
hitSeqTotal = 0;
hitSeqSample = 0;
rehashTotal = 0;
rehashSample = 0;
cout.flush();
}
} // End report every N seqs.
} // End foreach sequence in buffer
if (verbose) {
form1("%-5d: processed seq buffer\n",myId);
cout.flush();
}
} catch (BadSeqException &b) {
cerr << "ERROR: error processing seq id=" << ii << endl;
cerr.flush();
MPI_Abort(MPI_COMM_WORLD,368);
}
} // while more seqs expected
form4("%-5d: done processStage1Par used %d buckets for %d seqs %d hasheventsG\n",
myId,usedCount,readCount,hashEventsG);
// form1("%-5d: deleting line\n",myId);
delete [] line;
// form1("%-5d: deleting buffer\n",myId);
delete [] bCastBuffer;
form1("%-5d: done processFile1par\n",myId);
cout.flush();
} // End process file
//
// Process one file against hash table
//
void processFileStage2(
// Bucket *bucket,
// uint32 hSize,
FILE *inf2,
// ofstream &outf,
uint32 readCount,
bool expectQual,
uint32 indexOffset,
ClipInfoMap &clipInfo,
vector<ClipInfo2> &clipInfoB
)
{
uint32 i;
off_t startPos=0;
basic_string<char> line;
basic_string<char> name;
lastLine[0] = 0;
currLine[0] = 0;
cout << "main : start processStage\n";
cout.flush();
form1("allocating %d bytes for bCastBuffer\n",bCastBufferSz);
char * bCastBuffer= new char [bCastBufferSz];
uint32 buffUsed = 0;
uint32 seqsSent = 0;
for(i=0;i<readCount;i++) {
if (!readFasta(inf2,name,line,startPos)) {
cout << "EOF2 reached prematurely: " << i << " of " << readCount << " reads " << endl;
MPI_Abort(MPI_COMM_WORLD,968);
}
//
// Check clipping info on this one
//
if (haveTextClip) {
//
// Check clipping info on this one
//
ClipInfoMap::iterator p;
p=clipInfo.find(name);
if (p != clipInfo.end() ) {
ClipInfo &c = p->second;
if (!c.good) {
line = "XX";
} else if (c.right-(int32)wordSize <= c.left) {
line = "XX";
} else {
line = line.substr(c.left+1,c.right-c.left - 1);
}
} else {
if (verbose) cerr << "no clip info '" << name << "'" << endl;
}
} else if (haveBinClip) {
if (!clipInfoB[i].good) {
line = "XX";
} else if (clipInfoB[i].right-(int32)wordSize <= clipInfoB[i].left) {
line = "XX";
} else {
line = line.substr(
clipInfoB[i].left+1,
clipInfoB[i].right-clipInfoB[i].left-1);
}
}
uint32 len = line.length();
// Send buffer if it's overly full
if (buffUsed + len >= bCastBufferSz - 3) {
// Time to send buffer
// cout << "sending buffer with " << seqsSent << " seqs" << endl; cout.flush();
bCastBuffer[buffUsed++] = 0; // terminate buffer here
bCastBuffer[buffUsed] = 0; // more to come
MPI_Bcast((void *)bCastBuffer,bCastBufferSz,MPI_CHAR,G_master,MPI_COMM_WORLD);
seqsSent = 0;
buffUsed = 0;
}
// cout << "DDD" << endl; cout.flush();
// Add this seq entry to buffer
strcpy(bCastBuffer+buffUsed,line.c_str());
seqsSent++;
// cout << "EEE" << endl; cout.flush();
buffUsed += len + 1;
}
//
// Finalise buffer
bCastBuffer[buffUsed] = 0; // terminate buffer
bCastBuffer[buffUsed+1] = 1; // no more to come
cout << "main : flushing final send buffer with " << seqsSent << " seqs" << endl;
cout.flush();
MPI_Bcast((void *)bCastBuffer,bCastBufferSz,MPI_CHAR,G_master,MPI_COMM_WORLD);
delete [] bCastBuffer;
cout << "main :finished sending seq stage 2\n";
cout.flush();
}
int nextOrient(int orient)
{
if (orient == SAME) return OPPOSITE;
if (orient == OPPOSITE) return DONE;
assert(0);
return DONE;
}
void processFileStage2Par(
uint32 myId,
Bucket *bucket,
uint32 hSizeG,
uint32 hBot,
uint32 hTop,
ofstream &outf,
uint32 readCount
)
{
uint32 i,j,m;
uint32 hSize = hTop - hBot + 1;
int orientation;
bool seqFound = false;
uint32 lineSize = 100000;
form1("%-5d: start processStage2Par\n",myId);
cout.flush();
char *line = new char [lineSize];
char *bCastBuffer= new char [bCastBufferSz];
#if defined(USE_SPLIT_HKEY)
//
// Part of split key implementation
//
// Positions for skip characters
uint32 left = wordSize/4;
uint32 right = wordSize - 1 - left;
uint32 mid = wordSize / 2;
#endif
#if defined(USE_CONTAM_CODE)
//
// In this pass we want to note all reads that contain contaminated k-mers
//
vector<uchar> contamCount(readCount);
vector<uchar> contamAttempts(readCount);
vector<uchar> contamType(readCount);
for(i=0;i<readCount;i++) {
contamCount[i] = 0;
contamType[i] = 0;
contamAttempts[i] = 0;
}
uint32 contamHits = 0;
bool thisReadContam;
uint32 contamReads = 0;
#endif
for(i=0;i<readCount;) {
MPI_Bcast(bCastBuffer,bCastBufferSz,MPI_CHAR,G_master,MPI_COMM_WORLD);
if (verbose) {
form1("%d: stage 2 got buffer\n",myId);
cout.flush();
}
// if (verbose) form1("%d: stage 2 got seq\n",myId);
// line[len] = 0;
uint32 len;
for(m=0;bCastBuffer[m];m+=len + 1) {
if (i== readCount) {
form2("%-5d: Hit end of data by count i=%d\n",myId,i);
cout.flush();
if (verbose) {
cout << "Hit end of data by count" << endl;
cout.flush();
}
break;
}
if (bCastBuffer[m] == 1) {
form2("%-5d: Hit end of data marker i=%d\n",myId,i);
cout.flush();
if (verbose) {
cout << "Hit end of data marker" << endl;
cout.flush();
}
break;
}
assert(strlen(bCastBuffer+m)<lineSize-1); // don't overflow line buffer
strcpy(line,bCastBuffer+m);
len = strlen(line);
i++; // WARNING - i already incremented here, so read is i-1
#if defined(USE_CONTAM_CODE)
thisReadContam = false;
#endif
for(orientation=SAME;orientation != DONE;orientation = nextOrient(orientation)) {
// uint32 v = 0;
// uint32 v2 = 0;
basic_string<char> seq;
if (orientation==SAME) {
seq = line;
} else {
string tmp = line;
if (useColorSpace) {
seq = revOnly(tmp); // ABI specific color space assem
} else {
seq = revcomp(tmp);
}
}
if (seq.size() < wordSize) continue;
for(j=0;j<=seq.size()-wordSize;j++) {
#if defined(USE_SPLIT_HKEY)
unsigned char canonicalOrig[wordSize+1];
unsigned char canonical[wordSize+1];
// Loop here over different split keys
//
// we're going to assume that (sizeof(key) - 3) is a multiple of 4
//
// e.g 15 / 2 -> 7 15/4 -> 3 15-15/4 -1 = 11
//
//
// 0123456789ABCDE
// 00
// * 11
// * 10
// * 01
//
// make a specific copy of the key
// modify the sequence accordingly
uint32 sk;
for(sk=0;sk<4;sk++) {
switch(sk) {
case 0:
// Make copy, leave unchanged
strncpy((char *)canonicalOrig,(char *)seq.c_str()+j,wordSize);
strncpy((char *)canonical,(char *)canonicalOrig,wordSize);
break;
case 1:
// Retrieve copy and modify right
strncpy((char *)canonical,(char *)canonicalOrig,wordSize);
canonical[right] = '*';
break;
case 2:
// Retrieve copy and modify left
strncpy((char *)canonical,(char *)canonicalOrig,wordSize);
canonical[left] = '*';
break;
case 3:
// Retrieve copy and modify center
strncpy((char *)canonical,(char *)canonicalOrig,wordSize);
canonical[mid] = '*';
break;
default:
assert(0); // only 4 options
} // End set up key
//
// WARNING .. outdent here since the following code is conditionally in the loop above.
//
#endif
//
// Create a hash value from wordSize letters of DNA
//
// if (strncmp(seq.c_str() + j,"ATAAGATTTATAACTAATTTCAAAT",25) == 0) {
// cout << myId << "YES2 n=" << i-1 << endl;
// }
MD5_CTX context;
unsigned char digest[16];
MD5Init (&context);
const char *p = seq.c_str() + j;
// ---------------------------------------------------------------------------
#if defined(USE_SPLIT_HKEY)
// Just use the constructed canonical mer
MD5Update (&context, canonical, wordSize);
#else
// Use a section of the input sequence directly
MD5Update (&context, (unsigned char *)p, wordSize);
#endif
// ---------------------------------------------------------------------------
// Finish up digest
MD5Final (digest, &context);
// vvvvvvvvvvvv---------------------------------------------------------------
#if defined(DEBUG_KMER)
// Debugging
string tmp;
for(int zzz=0;zzz<wordSize;zzz++) tmp += p[zzz];
if (tmp == "GACTTCCTATTAAA" || tmp == "TTTAATAGGAAGTC") {
cout << "XXX2: i=" << i-1 << " same=" << (orientation==SAME) << " " <<
tmp << endl;
}
#endif
// ^^^^^^^^^^^^^--------------------------------------------------------------
uint32 k;
// vvvvvvvvvvvv---------------------------------------------------------------
#if defined(STORE_KMER)
#if defined(USE_SPLIT_HKEY)
// For split keys, we tack on 2 bits to say which type of key we
// are dealing with so we can print them out sensibly later
unsigned long long kMerBits = seq2Bits((const unsigned char *)canonical,wordSize);
kMerBits <<= 2; // move key left
kMerBits |= sk; // add on 2 bit key
#else
// Just normal sequence
unsigned long long kMerBits = seq2Bits((const unsigned char *)p,wordSize);
#endif
#endif
// ^^^^^^^^^^^^^--------------------------------------------------------------
unsigned long long *tmp1 = (unsigned long long *)digest;
unsigned long long v = tmp1[0] ^ tmp1[1];
if (true) {
// k = (v^v2) % hSizeG;
k = v % hSizeG;
if (k<hBot || k > hTop) continue;
k -=hBot;
bool done = false;
bool found = false;
do {
rehashTotal++;
k = (k + 1) % hSize;
switch(bucket[k].flag) {
case USED:
case SKIPPED:
#if defined(STORE_KMER)
if ( bucket[k].v == kMerBits )
#else
if ( bucket[k].v == v)
#endif
{
found = true;
done = true;
} else {
if (bucket[k].flag == USED) {
done = true;
}
}
break;
case UNUSED:
done = true;
found = false;
break;
default:
cerr << "ERROR: illegal flag '" <<
(int)(bucket[k].flag) <<"'" << endl;
exit(1);
} // End switch on bucket type
} while (!done);
rehashSample++;
hitSample++;
#if defined(USE_CONTAM_CODE)
//
// Track how many k-mres could have been positive
if (contamAttempts[i-1] < 254) contamAttempts[i-1]++;
#endif
if (found) {
seqFound = true;
hitTotal++;
hitCum++;
#if defined(USE_CONTAM_CODE)
//
// Did we hit a contaminant?
//
if (bucket[k].isContam) {
// i is pre incremented above, so i-1 is read #
if (contamCount[i-1] < 254) contamCount[i-1]++;
contamType[i-1] = bucket[k].isContam;
if (debugContam) {
char tmp[wordSize + 1];
strncpy(tmp,p,wordSize);
tmp[wordSize] = 0;
if (debugContam) {
cout << "CONTAMHIT: " << i-1 << " " << (int) contamCount[i-1] << " " << tmp <<
" " << (int32)bucket[k].isContam << endl;
}
}
contamHits++;
if (!thisReadContam) {
thisReadContam = true;
contamReads++;
}
}
#endif
#if defined(STORE_KMER)
assert(bucket[k].v == kMerBits);
#else
assert(bucket[k].v == v);
#endif
bucket[k].count++;
// Write out the seq number
uint32 tmp = i-1; // -1 b/c of pre incrementing i above
assert(tmp < readCount);
outf.write((char *)(&tmp),4);
// Write value (hash) position that was hit
outf.write((char *)&k,4);
// Offset is the distance between the first
// base in a read.
int offset;
//
// SAME Orientation
// |------ j ----->
// *=============VV========> (query read)
// |----offset --->
// OPP Orientation
// <-- offset --|
// <=============WW============| (query read)
// ==============VV============> (query read RC)
// |---- j ------->
//
//
//
if (orientation == SAME) {
offset = j;
} else {
offset = seq.size()-j-1;
}
assert(offset >=0 && (unsigned)offset <seq.size());
outf.write((char *)&offset,4);
unsigned char orientByte = orientation & 255;
outf.write((char *)&orientByte,1);
// if (myId == 1 && k==75370) {
// form4("hit %d %d %d %d %d\n",
// i,bucket[k].count,offset,(int)orientByte);
//
// }
// cout << "+"; cout.flush();
// Debugging
#if defined(DEBUG_READ)
if ((int32)tmp == debug1) {
form6("%-5d: XXXX readId=%d k=%d kmer=%s off=%d orient=%d\n",
myId,tmp,k,
bits2Seq(v,wordSize).c_str(),
offset,
orientation==SAME
);
}
#endif
}
} // End if enough bits loaded
} // End look through seq
} // End forward and reverse search
hitSeqSample++;
if (seqFound) {
hitSeqTotal++;
}
#if defined(USE_SPLIT_HKEY)
} // End each split key variant
#endif
//
// Stats to keep the human amused..
//
if (i % 10000 == 0 || i==readCount-1) {
if (hitSample > 0 && rehashSample > 0) {
cout.width(5);
cout << myId << ":"; // daz removed left output
cout <<
" pass2: " << i <<
// " usedbuckets= " << usedCount <<
" avgHits= " << hitTotal / (float)hitSample <<
// " readsWithHits= " << hitSeqTotal / (float)hitSeqSample <<
// " rehashTotal= " << rehashTotal <<
" hitsTotal= " << hitCum << endl;
// " " << rehashTotal / (float)rehashSample << endl;
hitTotal = 0;
hitSample = 0;
hitSeqTotal = 0;
hitSeqSample = 0;
rehashTotal = 0;
rehashSample = 0;
}
}
} // End foreach buffer receive
} // End foreach sequence
form1("%-5d: exitting seq processing\n",myId);
cout.flush();
#if defined(USE_CONTAM_CODE)
if (debugContam) {
for(i=0;i<readCount;i++) {
if (contamCount[i]) {
cout << "CONTAMREAD: " << i << " " << (int)contamCount[i] << " " << (int)contamAttempts[i] << endl;
}
}
cout << "CONTAM: total contamHits " << contamHits << endl;
cout << "CONTAM: total contamReads " << contamReads << endl;
}
// Send in final contaminant count for each read
form1("%-5d: send contam count\n",myId);
MPI_Send(&contamCount[0],readCount,MPI_CHAR,G_master,623,MPI_COMM_WORLD);
form1("%-5d: send contam attempts\n",myId);
MPI_Send(&contamAttempts[0],readCount,MPI_CHAR,G_master,624,MPI_COMM_WORLD);
form1("%-5d: send contam types\n",myId);
MPI_Send(&contamType[0],readCount,MPI_CHAR,G_master,625,MPI_COMM_WORLD);
#endif
delete [] line;
delete [] bCastBuffer;
form1("%-5d: done processFile2par\n",myId);
cout.flush();
} // End pass2 Parallel
void worker(string &projRoot,uint32 nProc,uint32 myId,string &outCachePath)
{
//
// Emits a hash table, plus hits to that table, and also a lookup
// table of offsets for the fasta file
//
char tmp[1000];
sprintf(tmp,"%s/%s.%d",outCachePath.c_str(),projRoot.c_str(),G_PXL_P2J[myId]);
string projRootOrig = projRoot;
projRoot = tmp;
cout << "ProjRoot" << projRoot << endl;
ofstream outf((projRoot + ".hits").c_str(),ios::binary | ios::out);
ofstream outf3((projRoot + ".hash").c_str(),ios::binary | ios::out);
uint32 hSize = 0;
uint32 hTop = 0;
uint32 hBot = 0;
uint32 readCount;
uint32 i,k;
MPI_Bcast(&hSize,1,MPI_UNSIGNED,G_master,MPI_COMM_WORLD);
form2("%-5d: hSize=%d\n",myId,hSize);
MPI_Bcast(&genomeCoverage,1,MPI_FLOAT,G_master,MPI_COMM_WORLD);
form2("%-5d: genomeCoverage=%f\n",myId,genomeCoverage);
MPI_Bcast(&readCount,1,MPI_UNSIGNED,G_master,MPI_COMM_WORLD);
form2("%-5d: readCount=%d\n",myId,readCount);
k = 0;
for(i=1;i<=G_PXL_P2J[myId];i++) {
// form2("%d:A k=%d\n",myId,k);
if (i==G_PXL_P2J[myId]) hBot = k;
// form2("%d:B k=%d\n",myId,k);
k = k + (int)(hSize/(float)(nProc - 1));
// form2("%d:C k=%d\n",myId,k);
if (i==G_PXL_P2J[myId]) hTop = k;
k++;
}
// uint32 hSizeG = hSize;
// hSizeG = hSizeG;
hSizeG = hSize;
hSize = hTop - hBot + 1;
// form4("xxx %d: hSize=%d hBot=%d hTop=%d\n",myId,hSize,hBot,hTop);
//
// This is how big the hash table is
//
form3("%-5d: about to allocate %d buckets %lld bytes\n",myId,hSize,sizeof(Bucket) * (long long)hSize);
cout.flush();
if (sizeof(Bucket) * (long long)hSize > 2147483647 && sizeof (int *) == 4) {
cerr << "WARNING: memory allocation likely to fail" << endl;
cerr << "WARNING: try using more CPUs to spread the hash table across more machines" << endl;
}
Bucket *bucket = new Bucket[hSize];
form1("%-5d: done allocation\n",myId);
// Starts out empty
for(i=0;i<hSize;i++) {
bucket[i].flag = UNUSED;
bucket[i].count = 0;
bucket[i].isContam = false;
}
processFileStage1Par(projRootOrig,myId,bucket,hSizeG,hBot,hTop,readCount);
processFileStage2Par(myId,bucket,hSizeG,hBot,hTop,outf,readCount);
form3("%-5d: Dump hash table %d buckets %lld bytes\n",myId,hSize,sizeof(Bucket) * (long long)hSize);
cout.flush();
for(i=0;i<hSize;i++) {
// if (bucket[i].flag == USED || bucket[i].flag == SKIPPED) {
// cout << i << " " << bucket[i].flag << endl;
#if 0
// We aren't preserving actual hash values
outf3.write((char *)&(bucket[i].v),4);
outf3.write((char *)&(bucket[i].v2),4);
#endif
outf3.write((char *)&(bucket[i].count),4);
// }
}
form1("%-5d: hash table dumped\n",myId);
cout.flush();
#if defined(DEBUG)
// Write out all hash table
ofstream outfDump("hashdump");
cout << "Writing out full hash table to hashdump " << (sizeof(Bucket) * hSize) << " bytes" << endl;
cout.flush();
outfDump.write((char *)bucket,sizeof(Bucket) * hSize);
#endif
//
// Dump all kmers out - debugging purpose
//
if (fullKMerDump) {
sprintf(tmp,"%s.kmer.%d",projRootOrig.c_str(),myId);
form2("%-5d: Dump all k-mers to %s\n",myId,tmp);
ofstream outf5(tmp);
if (!outf5) {
cerr << "ERROR: unable to open '" << tmp << "' for writing " << endl;
MPI_Abort(MPI_COMM_WORLD,822);
}
for(i=0;i<hSize;i++) {
#if defined(USE_CONTAM_CODE)
if (bucket[i].isContam) continue;
#endif
if (bucket[i].flag == USED || bucket[i].flag == SKIPPED) {
outf5 << i << " " << bucket[i].count << " " << bits2Seq(bucket[i].v,wordSize) << endl;
} // End if occupied bucket
} // End for each hash table entry
}
#if defined(DUMP_FREQSEQS)
form1("%-5d: Dump frequent k-mers\n",myId);
cout.flush();
sprintf(tmp,"%s.repmer.%d",projRootOrig.c_str(),myId);
ofstream outf4(tmp);
if (!outf4) {
cerr << "ERROR: unable to open '" << tmp << "' for writing " << endl;
cerr.flush();
MPI_Abort(MPI_COMM_WORLD,822);
}
uint32 repThresh = (uint32)(genomeCoverage * 2.0);
if (repMerMinDepth != 99999999) {
repThresh = repMerMinDepth;
}
if (dumpWholeHashTable) {
repThresh = 1;
}
uint32 numRepKmer = 0;
uint32 numKmer = 0;
for(i=0;i<hSize;i++) {
#if defined(USE_CONTAM_CODE)
if (bucket[i].isContam) continue;
#endif
if (bucket[i].flag == USED || bucket[i].flag == SKIPPED) {
numKmer++;
if (bucket[i].count >= repThresh ) {
outf4 << i << " " << bucket[i].count << " " << bits2Seq(bucket[i].v,wordSize) << endl;
numRepKmer++;
// outf3.write((char *)&(bucket[i].count),4);
}
} // End if occupied bucket
} // End for each hash table entry
MPI_Send(&numKmer,1,MPI_UNSIGNED,G_master,620,MPI_COMM_WORLD);
MPI_Send(&numRepKmer,1,MPI_UNSIGNED,G_master,621,MPI_COMM_WORLD);
// End dump overly frequent k-mers
#endif
form1("%-5d: hash dump done\n",myId);
cout.flush();
//
// Calculate hash histogram and predict graph edge requirements
// for next stage
// -------------------------------------------------------------
vector<uint32> histo(maxHisto);
for(i=0;i<maxHisto;i++) {
histo[i] = 0;
}
uint32 expectedMaxEdges = 0;
for(i=0;i<hSize;i++) {
#if defined(USE_CONTAM_CODE)
if (bucket[i].isContam) continue;
#endif
histo[min(maxHisto-1,bucket[i].count)]++;
if (bucket[i].count && bucket[i].count < genomeCoverage * 3) {
expectedMaxEdges += (bucket[i].count - 1) * bucket[i].count;
}
}
// End in to master
form2("%-5d: sending in expectedMaxEdges=%d\n",myId,expectedMaxEdges); cout.flush();
MPI_Send(&expectedMaxEdges,1,MPI_INT,G_master,626,MPI_COMM_WORLD);
// Bring histograms together
cout.flush();
form2("%-5d: sending in %d histo entries\n",myId,maxHisto); cout.flush();
MPI_Send(&histo[0],maxHisto,MPI_INT,G_master,627,MPI_COMM_WORLD);
form1("%-5d: exitting worker()\n",myId);
}
int selectPrime(bool *isPrime,int preferred,int upto) {
int i;
for (i=preferred;i<upto;i++) {
if (isPrime[i]) break;
}
return i;
}
#if defined(USE_CONTAM_CODE)
//-----------------------------------------------------
//
// Contamination screening system
//
//-----------------------------------------------------
char *loadContamFasta(
string &contamFile,
vector<ContamEntry> &entry
)
{
uint32 i,j;
// Find out size of file
struct stat st;
stat(contamFile.c_str(),&st);
uint32 size2 = st.st_size;
// Open handle
int fildes2 = open(contamFile.c_str(),O_RDONLY);
if (fildes2 == -1) {
cerr << endl << endl << "ERROR: Failed to open fasta file '" << contamFile << "'" << endl;
return 0;
}
//
// Map file into memory
//
int mapSize = size2;
if (debugContam) cout << "CONTAM: Allocating " << (mapSize) << " bytes for genome seq storage" << endl;
char *genome1 = (char *)mmap((caddr_t) 0, mapSize, (PROT_READ),MAP_PRIVATE,
fildes2, 0);
//
// Make second copy of the sequence for the wrapped version
//
cout << "Mapping contam fasta " << mapSize << ": bytes " << endl;
char *genome = new char[mapSize];
for(i=0,j=0;i<(uint32)mapSize;i++) {
if (genome1[i] == '>') {
// Update last length entry if there is a last one
if (entry.size()) {
ContamEntry &last = entry[entry.size()-1];
last.len = j - last.off;
}
// Note header sequence
string name;
bool seenSpace = false;
for(i++;i<(uint32)mapSize && genome1[i] != '\n' && genome1[i] != '\r';i++) {
if (seenSpace) continue;
if (genome1[i] == ' ') {
seenSpace = true;
continue;
}
name += genome1[i];
}
ContamEntry ce(name,j,0);
entry.push_back(ce);
}
if (genome1[i] == 'G' || genome1[i] == 'A' ||
genome1[i] == 'T' || genome1[i] == 'C')
{
genome[j++] = genome1[i];
} else if (genome1[i] == 'g' || genome1[i] == 'a' ||
genome1[i] == 't' || genome1[i] == 'c')
{
genome[j++] = toupper(genome1[i]);
}
}
// Update last length entry if there is a last one
if (entry.size()) {
ContamEntry &last = entry[entry.size()-1];
last.len = j - last.off;
}
// Release memory mapped qual scores after we are done with them
munmap(genome1,mapSize);
genome[j] = 0;
return genome;
}
bool injectContam(
int myId,
string &fileName,
Bucket *bucket,
uint32 hSize,
uint32 hSizeG,
uint32 hBot,
uint32 hTop
)
{
vector<ContamEntry> entry;
// Load entries
char *genome = loadContamFasta(fileName,entry);
if (!genome) return false;
//
// Foreach contaminant entry, insert the k-mers into the
// hash table
//
uint32 i,j;
// uint32 l = strlen(genome);
uint32 added = 0;
for(i=0;i<entry.size();i++) {
uint32 end = entry[i].off + entry[i].len;
char x = genome[end];
genome[end] = 0;
string thisContam(genome+entry[i].off);
uchar whichContam = (i+1) < 254?(uchar)i+1: 255;
// Extract section for this contaminant entry
if (debugContam) {
cout << "CONTAM: " << entry[i].off << " " << entry[i].len << " " << entry[i].name << endl;
}
// cout << endl << thisContam << endl;
// Can't process really short entries, so skip
if (thisContam.size() < wordSize) continue; // CONTINUE if too short
uint32 len = thisContam.size();
for(j=(len%wordSize)/2;j<=len-wordSize;j+=kMerContamTile) {
//
// We want to store only canonical hash values, so we store
// the reverse complement seqeuence if it is the canonical form?
//
unsigned char canonical[wordSize+1]; // dp - +1 - don't know why this worked I think this gets null terminated?
makeCanonical(canonical,(unsigned char *)(thisContam.c_str() + j));
//
// *******************************************
// Key placement of word into hash table
// *******************************************
//
assert(canonical[0] != 0);
storeWord(bucket,hSize,canonical,hSizeG,hBot,hTop,whichContam);
added++;
} // End look through seq
// Undo temporary null termination of this record
genome[end] = x;
}
form3("%-5d: Loaded contaminants: %d words usedCount= %d\n",myId,added,usedCount);
delete [] genome;
return true;
}
#endif
Thanx
P
Comment