/***********************************************************************/
Maligner::Maligner(vector<Sequence*> temp, int num, int match, int misMatch, float div, int ms, int minCov, string mode, Database* dataLeft, Database* dataRight) :
- db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov), searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight) {}
+ db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov), searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight) { m = MothurOut::getInstance(); }
/***********************************************************************/
string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
try {
refSeqs = getBlastSeqs(query, numWanted); //fills indexes
}else if (searchMethod == "kmer") {
refSeqs = getKmerSeqs(query, numWanted); //fills indexes
- }else { mothurOut("not valid search."); exit(1); } //should never get here
+ }else { m->mothurOut("not valid search."); exit(1); } //should never get here
+
+ if (m->control_pressed) { return chimera; }
refSeqs = minCoverageFilter(refSeqs);
//fills outputResults
chimera = chimeraMaligner(chimeraPenalty, decalc);
-
+
+ if (m->control_pressed) { return chimera; }
//free memory
delete query;
return chimera;
}
catch(exception& e) {
- errorOut(e, "Maligner", "getResults");
+ m->errorOut(e, "Maligner", "getResults");
exit(1);
}
}
//trims seqs to first non gap char in all seqs and last non gap char in all seqs
spotMap = decalc->trimSeqs(query, refSeqs);
+
+ //you trimmed the whole sequence, skip
+ if (query->getAligned() == "") { return "no"; }
vector<Sequence*> temp = refSeqs;
temp.push_back(query);
vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
+ if (m->control_pressed) { return chimera; }
+
fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
vector<score_struct> path = extractHighestPath(matrix);
+ if (m->control_pressed) { return chimera; }
+
vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
if (trace.size() > 1) { chimera = "yes"; }
percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
+ if (m->control_pressed) { return chimera; }
+
//save output results
for (int i = 0; i < trace.size(); i++) {
int regionStart = trace[i].col;
return chimera;
}
catch(exception& e) {
- errorOut(e, "Maligner", "chimeraMaligner");
+ m->errorOut(e, "Maligner", "chimeraMaligner");
exit(1);
}
}
//if coverage above minimum
if (coverage > minCoverage) {
newRefs.push_back(ref[i]);
+ }else {
+ delete ref[i];
}
}
return newRefs;
}
catch(exception& e) {
- errorOut(e, "Maligner", "minCoverageFilter");
+ m->errorOut(e, "Maligner", "minCoverageFilter");
exit(1);
}
}
}
catch(exception& e) {
- errorOut(e, "Maligner", "computeChimeraPenalty");
+ m->errorOut(e, "Maligner", "computeChimeraPenalty");
exit(1);
}
}
spotMap = newMap;
}
catch(exception& e) {
- errorOut(e, "Maligner", "verticalFilter");
+ m->errorOut(e, "Maligner", "verticalFilter");
exit(1);
}
}
return m;
}
catch(exception& e) {
- errorOut(e, "Maligner", "buildScoreMatrix");
+ m->errorOut(e, "Maligner", "buildScoreMatrix");
exit(1);
}
}
//***************************************************************************************************************
-void Maligner::fillScoreMatrix(vector<vector<score_struct> >& m, vector<Sequence*> seqs, int penalty) {
+void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
try{
//get matrix dimensions
//are you both gaps?
if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
- m[i][0].score = 0;
+ ms[i][0].score = 0;
}else if (queryAligned[0] == subjectAligned[0]) {
- m[i][0].score = matchScore;
+ ms[i][0].score = matchScore;
}else{
- m[i][0].score = 0;
+ ms[i][0].score = 0;
}
}
//compute score based on previous columns scores
for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
- int sumScore = matchMisMatchScore + m[prevIndex][j-1].score;
+ int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
//you are not at yourself
if (prevIndex != i) { sumScore += penalty; }
if (sumScore < 0) { sumScore = 0; }
- if (sumScore > m[i][j].score) {
- m[i][j].score = sumScore;
- m[i][j].prev = prevIndex;
+ if (sumScore > ms[i][j].score) {
+ ms[i][j].score = sumScore;
+ ms[i][j].prev = prevIndex;
}
}
}
}
catch(exception& e) {
- errorOut(e, "Maligner", "fillScoreMatrix");
+ m->errorOut(e, "Maligner", "fillScoreMatrix");
exit(1);
}
}
//***************************************************************************************************************
-vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > m) {
+vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
try {
//get matrix dimensions
int numCols = query->getAligned().length();
- int numRows = m.size();
+ int numRows = ms.size();
//find highest score scoring matrix
for (int i = 0; i < numRows; i++) {
for (int j = 0; j < numCols; j++) {
- if (m[i][j].score > highestScore) {
- highestScore = m[i][j].score;
- highestStruct = m[i][j];
+ if (ms[i][j].score > highestScore) {
+ highestScore = ms[i][j].score;
+ highestStruct = ms[i][j];
}
}
}
int score = highestStruct.score;
while (pos >= 0 && score > 0) {
- score_struct temp = m[rowIndex][pos];
+ score_struct temp = ms[rowIndex][pos];
score = temp.score;
if (score > 0) { path.push_back(temp); }
}
catch(exception& e) {
- errorOut(e, "Maligner", "extractHighestPath");
+ m->errorOut(e, "Maligner", "extractHighestPath");
exit(1);
}
}
}
catch(exception& e) {
- errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
+ m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
exit(1);
}
}
return chimera;
}
catch(exception& e) {
- errorOut(e, "Maligner", "constructChimericSeq");
+ m->errorOut(e, "Maligner", "constructChimericSeq");
exit(1);
}
}
try {
if (queryAlign.length() != chimera.length()) {
- mothurOut("Error, alignment strings are of different lengths: "); mothurOutEndLine();
- mothurOut(toString(queryAlign.length())); mothurOutEndLine(); mothurOutEndLine(); mothurOutEndLine(); mothurOutEndLine();
- mothurOut(toString(chimera.length())); mothurOutEndLine();
+ m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine();
+ m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine();
+ m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
return -1.0;
}
}
catch(exception& e) {
- errorOut(e, "Maligner", "computePercentID");
+ m->errorOut(e, "Maligner", "computePercentID");
exit(1);
}
}
try {
indexes.clear();
vector<Sequence*> refResults;
- //generate blastdb
- Database* database = new BlastDB(-2.0, -1.0, matchScore, misMatchPenalty);
- for (int i = 0; i < db.size(); i++) { database->addSequence(*db[i]); }
- database->generateDB();
- database->setNumSeqs(db.size());
-
+
//get parts of query
string queryUnAligned = q->getUnaligned();
string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence
string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence
- Sequence* queryLeft = new Sequence(q->getName(), leftQuery);
- Sequence* queryRight = new Sequence(q->getName(), rightQuery);
-
- vector<int> tempIndexesRight = database->findClosestMegaBlast(queryRight, num+1);
- vector<int> tempIndexesLeft = database->findClosestMegaBlast(queryLeft, num+1);
+ Sequence* queryLeft = new Sequence(q->getName()+"left", leftQuery);
+ Sequence* queryRight = new Sequence(q->getName()+"right", rightQuery);
- //if ((tempIndexesRight.size() != (num+1)) || (tempIndexesLeft.size() != (num+1))) { mothurOut("megablast returned " + toString(tempIndexesRight.size()) + " results for the right end, and " + toString(tempIndexesLeft.size()) + " for the left end. Needed " + toString(num+1) + ". Unable to porcess sequence " + q->getName()); mothurOutEndLine(); return refResults; }
+ vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1);
+ vector<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1);
+
+ //if ((tempIndexesRight.size() == 0) && (tempIndexesLeft.size() == 0)) { m->mothurOut("megablast returned " + toString(tempIndexesRight.size()) + " results for the right end, and " + toString(tempIndexesLeft.size()) + " for the left end. Needed " + toString(num+1) + ". Unable to process sequence " + q->getName()); m->mothurOutEndLine(); return refResults; }
vector<int> smaller;
vector<int> larger;
//cout << endl;
delete queryRight;
delete queryLeft;
- delete database;
-
+
return refResults;
}
catch(exception& e) {
- errorOut(e, "Maligner", "getBlastSeqs");
+ m->errorOut(e, "Maligner", "getBlastSeqs");
exit(1);
}
}
return refResults;
}
catch(exception& e) {
- errorOut(e, "Maligner", "getBlastSeqs");
+ m->errorOut(e, "Maligner", "getBlastSeqs");
exit(1);
}
}
projects / mothur.git / blobdiff