#endif
}
//***************************************************************************************************************
-Sequence* Ccode::print(ostream& out, ostream& outAcc) {
+Sequence Ccode::print(ostream& out, ostream& outAcc) {
try {
ofstream out2;
//free memory
for (int i = 0; i < closest.size(); i++) { delete closest[i].seq; }
- return NULL;
+ return *querySeq;
}
catch(exception& e) {
m->errorOut(e, "Ccode", "print");
}
#ifdef USE_MPI
//***************************************************************************************************************
-Sequence* Ccode::print(MPI_File& out, MPI_File& outAcc) {
+Sequence Ccode::print(MPI_File& out, MPI_File& outAcc) {
try {
string outMapString = "";
//free memory
for (int i = 0; i < closest.size(); i++) { delete closest[i].seq; }
- return NULL;
+ return *querySeq;
}
catch(exception& e) {
m->errorOut(e, "Ccode", "print");
~Ccode();
int getChimeras(Sequence* query);
- Sequence* print(ostream&, ostream&);
+ Sequence print(ostream&, ostream&);
#ifdef USE_MPI
- Sequence* print(MPI_File&, MPI_File&);
+ Sequence print(MPI_File&, MPI_File&);
#endif
private:
float dist;
int index;
};
+/***********************************************************************/
+struct SeqCompare {
+ Sequence seq;
+ float dist;
+ int index;
+};
//********************************************************************************************************************
//sorts lowest to highest
inline bool compareRegionStart(results left, results right){
return (left.dist < right.dist);
}
//********************************************************************************************************************
-
+//sorts lowest to highest
+inline bool compareSeqCompare(SeqCompare left, SeqCompare right){
+ return (left.dist < right.dist);
+}
+//********************************************************************************************************************
struct sim {
string leftParent;
string rightParent;
virtual void printHeader(ostream&){};
virtual int getChimeras(Sequence*){ return 0; }
virtual int getChimeras(){ return 0; }
- virtual Sequence* print(ostream&, ostream&){ return NULL; }
- virtual Sequence* print(ostream&, ostream&, data_results, data_results) { return NULL; }
+ virtual Sequence print(ostream&, ostream&){ Sequence temp; return temp; }
+ virtual Sequence print(ostream&, ostream&, data_results, data_results) { Sequence temp; return temp; }
virtual int print(ostream&, ostream&, string){ return 0; }
virtual data_results getResults() { data_results results; return results; }
#ifdef USE_MPI
- virtual Sequence* print(MPI_File&, MPI_File&){ return 0; }
- virtual Sequence* print(MPI_File&, MPI_File&, data_results, data_results){ return NULL; }
+ virtual Sequence print(MPI_File&, MPI_File&){ Sequence temp; return temp; }
+ virtual Sequence print(MPI_File&, MPI_File&, data_results, data_results){ Sequence temp; return temp; }
virtual int print(MPI_File&, MPI_File&, string){ return 0; }
#endif
}
}
//***************************************************************************************************************
-Sequence* ChimeraCheckRDP::print(ostream& out, ostream& outAcc) {
+Sequence ChimeraCheckRDP::print(ostream& out, ostream& outAcc) {
try {
m->mothurOut("Processing: " + querySeq->getName()); m->mothurOutEndLine();
}
}
- return NULL;
+ return *querySeq;
}
catch(exception& e) {
m->errorOut(e, "ChimeraCheckRDP", "print");
}
#ifdef USE_MPI
//***************************************************************************************************************
-Sequence* ChimeraCheckRDP::print(MPI_File& out, MPI_File& outAcc) {
+Sequence ChimeraCheckRDP::print(MPI_File& out, MPI_File& outAcc) {
try {
cout << "Processing: " << querySeq->getName() << endl;
}
}
- return NULL;
+ return *querySeq;
}
catch(exception& e) {
m->errorOut(e, "ChimeraCheckRDP", "print");
~ChimeraCheckRDP();
int getChimeras(Sequence*);
- Sequence* print(ostream&, ostream&);
+ Sequence print(ostream&, ostream&);
#ifdef USE_MPI
- Sequence* print(MPI_File&, MPI_File&);
+ Sequence print(MPI_File&, MPI_File&);
#endif
private:
realign = r;
trimChimera = trim;
- decalc = new DeCalculator();
-
doPrep();
}
catch(exception& e) {
trimChimera = trim;
priority = prior;
- decalc = new DeCalculator();
-
createFilter(templateSeqs, 0.0); //just removed columns where all seqs have a gap
if (searchMethod == "distance") {
}
}
//***************************************************************************************************************
-vector<Sequence*> ChimeraSlayer::getTemplate(Sequence* q, vector<Sequence*>& userTemplateFiltered) {
+vector<Sequence*> ChimeraSlayer::getTemplate(Sequence q, vector<Sequence*>& userTemplateFiltered) {
try {
//when template=self, the query file is sorted from most abundance to least abundant
//userTemplate grows as the query file is processed by adding sequences that are not chimeric and more abundant
vector<Sequence*> userTemplate;
- int myAbund = priority[q->getName()];
+ int myAbund = priority[q.getName()];
for (int i = 0; i < templateSeqs.size(); i++) {
//***************************************************************************************************************
ChimeraSlayer::~ChimeraSlayer() {
- delete decalc;
if (templateFileName != "self") {
if (searchMethod == "kmer") { delete databaseRight; delete databaseLeft; }
else if (searchMethod == "blast") { delete databaseLeft; }
out << "Name\tLeftParent\tRightParent\tDivQLAQRB\tPerIDQLAQRB\tBootStrapA\tDivQLBQRA\tPerIDQLBQRA\tBootStrapB\tFlag\tLeftWindow\tRightWindow\n";
}
//***************************************************************************************************************
-Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc) {
+Sequence ChimeraSlayer::print(ostream& out, ostream& outAcc) {
try {
- Sequence* trim = NULL;
- if (trimChimera) { trim = new Sequence(trimQuery.getName(), trimQuery.getAligned()); }
+ Sequence trim;
+ if (trimChimera) { trim.setName(trimQuery.getName()); trim.setAligned(trimQuery.getAligned()); }
if (chimeraFlags == "yes") {
string chimeraFlag = "no";
if (chimeraFlag == "yes") {
if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) {
- m->mothurOut(querySeq->getName() + "\tyes"); m->mothurOutEndLine();
- outAcc << querySeq->getName() << endl;
+ m->mothurOut(querySeq.getName() + "\tyes"); m->mothurOutEndLine();
+ outAcc << querySeq.getName() << endl;
- if (templateFileName == "self") { chimericSeqs.insert(querySeq->getName()); }
+ if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); }
if (trimChimera) {
int lengthLeft = chimeraResults[0].winLEnd - chimeraResults[0].winLStart;
int lengthRight = chimeraResults[0].winREnd - chimeraResults[0].winRStart;
- string newAligned = trim->getAligned();
+ string newAligned = trim.getAligned();
if (lengthLeft > lengthRight) { //trim right
for (int i = (chimeraResults[0].winRStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
}else { //trim left
for (int i = 0; i < chimeraResults[0].winLEnd; i++) { newAligned[i] = '.'; }
}
- trim->setAligned(newAligned);
+ trim.setAligned(newAligned);
}
}
}
printBlock(chimeraResults[0], chimeraFlag, out);
out << endl;
}else {
- out << querySeq->getName() << "\tno" << endl;
+ out << querySeq.getName() << "\tno" << endl;
}
return trim;
}
}
//***************************************************************************************************************
-Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftPiece, data_results rightPiece) {
+Sequence ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftPiece, data_results rightPiece) {
try {
- Sequence* trim = NULL;
+ Sequence trim;
if (trimChimera) {
string aligned = leftPiece.trimQuery.getAligned() + rightPiece.trimQuery.getAligned();
- trim = new Sequence(leftPiece.trimQuery.getName(), aligned);
+ trim.setName(leftPiece.trimQuery.getName()); trim.setAligned(aligned);
}
if ((leftPiece.flag == "yes") || (rightPiece.flag == "yes")) {
if (leftPiece.flag == "yes") { if ((leftPiece.results[0].bsa >= minBS) || (leftPiece.results[0].bsb >= minBS)) { leftChimeric = true; } }
if (rightChimeric || leftChimeric) {
- m->mothurOut(querySeq->getName() + "\tyes"); m->mothurOutEndLine();
- outAcc << querySeq->getName() << endl;
+ m->mothurOut(querySeq.getName() + "\tyes"); m->mothurOutEndLine();
+ outAcc << querySeq.getName() << endl;
- if (templateFileName == "self") { chimericSeqs.insert(querySeq->getName()); }
+ if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); }
if (trimChimera) {
- string newAligned = trim->getAligned();
+ string newAligned = trim.getAligned();
//right side is fine so keep that
if ((leftChimeric) && (!rightChimeric)) {
}
}
- trim->setAligned(newAligned);
+ trim.setAligned(newAligned);
}
}
printBlock(leftPiece, rightPiece, leftChimeric, rightChimeric, chimeraFlag, out);
out << endl;
}else {
- out << querySeq->getName() << "\tno" << endl;
+ out << querySeq.getName() << "\tno" << endl;
}
return trim;
#ifdef USE_MPI
//***************************************************************************************************************
-Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results leftPiece, data_results rightPiece) {
+Sequence ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results leftPiece, data_results rightPiece) {
try {
MPI_Status status;
bool results = false;
string outAccString = "";
string outputString = "";
- Sequence* trim = NULL;
+ Sequence trim;
if (trimChimera) {
string aligned = leftPiece.trimQuery.getAligned() + rightPiece.trimQuery.getAligned();
- trim = new Sequence(leftPiece.trimQuery.getName(), aligned);
+ trim.setName(leftPiece.trimQuery.getName()); trim.setAligned(aligned);
}
if (leftPiece.flag == "yes") { if ((leftPiece.results[0].bsa >= minBS) || (leftPiece.results[0].bsb >= minBS)) { leftChimeric = true; } }
if (rightChimeric || leftChimeric) {
- cout << querySeq->getName() << "\tyes" << endl;
- outAccString += querySeq->getName() + "\n";
+ cout << querySeq.getName() << "\tyes" << endl;
+ outAccString += querySeq.getName() + "\n";
results = true;
- if (templateFileName == "self") { chimericSeqs.insert(querySeq->getName()); }
+ if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); }
//write to accnos file
int length = outAccString.length();
delete buf2;
if (trimChimera) {
- string newAligned = trim->getAligned();
+ string newAligned = trim.getAligned();
//right side is fine so keep that
if ((leftChimeric) && (!rightChimeric)) {
}
}
- trim->setAligned(newAligned);
+ trim.setAligned(newAligned);
}
}
delete buf;
}else {
- outputString += querySeq->getName() + "\tno\n";
+ outputString += querySeq.getName() + "\tno\n";
//write to output file
int length = outputString.length();
}
}
//***************************************************************************************************************
-Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
+Sequence ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
try {
MPI_Status status;
bool results = false;
string outAccString = "";
string outputString = "";
- Sequence* trim = NULL;
- if (trimChimera) { trim = new Sequence(trimQuery.getName(), trimQuery.getAligned()); }
+ Sequence trim;
+ if (trimChimera) { trim.setName(trimQuery.getName()); trim.setAligned(trimQuery.getAligned()); }
if (chimeraFlags == "yes") {
string chimeraFlag = "no";
if (chimeraFlag == "yes") {
if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) {
- cout << querySeq->getName() << "\tyes" << endl;
- outAccString += querySeq->getName() + "\n";
+ cout << querySeq.getName() << "\tyes" << endl;
+ outAccString += querySeq.getName() + "\n";
results = true;
- if (templateFileName == "self") { chimericSeqs.insert(querySeq->getName()); }
+ if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); }
//write to accnos file
int length = outAccString.length();
int lengthLeft = chimeraResults[0].winLEnd - chimeraResults[0].winLStart;
int lengthRight = chimeraResults[0].winREnd - chimeraResults[0].winRStart;
- string newAligned = trim->getAligned();
+ string newAligned = trim.getAligned();
if (lengthLeft > lengthRight) { //trim right
for (int i = (chimeraResults[0].winRStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
}else { //trim left
for (int i = 0; i < (chimeraResults[0].winLEnd-1); i++) { newAligned[i] = '.'; }
}
- trim->setAligned(newAligned);
+ trim.setAligned(newAligned);
}
}
}
delete buf;
}else {
- outputString += querySeq->getName() + "\tno\n";
+ outputString += querySeq.getName() + "\tno\n";
//write to output file
int length = outputString.length();
chimeraFlags = "no";
printResults.flag = "no";
- querySeq = query;
+ querySeq = *query;
//you must create a template
vector<Sequence*> thisTemplate;
vector<Sequence*> thisFilteredTemplate;
if (templateFileName != "self") { thisTemplate = templateSeqs; thisFilteredTemplate = filteredTemplateSeqs; }
- else { thisTemplate = getTemplate(query, thisFilteredTemplate); } //fills this template and creates the databases
+ else { thisTemplate = getTemplate(*query, thisFilteredTemplate); } //fills this template and creates the databases
if (m->control_pressed) { return 0; }
if (thisTemplate.size() == 0) { return 0; } //not chimeric
//moved this out of maligner - 4/29/11
- vector<Sequence*> refSeqs = getRefSeqs(query, thisTemplate, thisFilteredTemplate);
+ vector<Sequence> refSeqs = getRefSeqs(*query, thisTemplate, thisFilteredTemplate);
Maligner maligner(refSeqs, match, misMatch, divR, minSim, minCov);
Slayer slayer(window, increment, minSim, divR, iters, minSNP, minBS);
if (m->control_pressed) { return 0; }
- string chimeraFlag = maligner.getResults(query, decalc);
+ string chimeraFlag = maligner.getResults(*query, decalc);
if (m->control_pressed) { return 0; }
vector<results> Results = maligner.getOutput();
- for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
+ //for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
if (chimeraFlag == "yes") {
// cout << query->getAligned() << endl;
//get sequence that were given from maligner results
- vector<SeqDist> seqs;
+ vector<SeqCompare> seqs;
map<string, float> removeDups;
map<string, float>::iterator itDup;
map<string, string> parentNameSeq;
for (itDup = removeDups.begin(); itDup != removeDups.end(); itDup++) {
itSeq = parentNameSeq.find(itDup->first);
- Sequence* seq = new Sequence(itDup->first, itSeq->second);
+ Sequence seq(itDup->first, itSeq->second);
- SeqDist member;
+ SeqCompare member;
member.seq = seq;
member.dist = itDup->second;
seqs.push_back(member);
//limit number of parents to explore - default 3
if (Results.size() > parents) {
//sort by distance
- sort(seqs.begin(), seqs.end(), compareSeqDist);
+ sort(seqs.begin(), seqs.end(), compareSeqCompare);
//prioritize larger more similiar sequence fragments
reverse(seqs.begin(), seqs.end());
- for (int k = seqs.size()-1; k > (parents-1); k--) {
- delete seqs[k].seq;
- seqs.pop_back();
- }
+ //for (int k = seqs.size()-1; k > (parents-1); k--) {
+ // delete seqs[k].seq;
+ //seqs.pop_back();
+ //}
}
//put seqs into vector to send to slayer
// cout << query->getAligned() << endl;
- vector<Sequence*> seqsForSlayer;
+ vector<Sequence> seqsForSlayer;
for (int k = 0; k < seqs.size(); k++) {
// cout << seqs[k].seq->getAligned() << endl;
seqsForSlayer.push_back(seqs[k].seq);
// cout << seqs[k].seq->getName() << endl;
}
- if (m->control_pressed) { for (int k = 0; k < seqs.size(); k++) { delete seqs[k].seq; } return 0; }
+ if (m->control_pressed) { return 0; }
//send to slayer
- chimeraFlags = slayer.getResults(query, seqsForSlayer);
+ chimeraFlags = slayer.getResults(*query, seqsForSlayer);
if (m->control_pressed) { return 0; }
chimeraResults = slayer.getOutput();
printResults.results = chimeraResults;
//free memory
- for (int k = 0; k < seqs.size(); k++) { delete seqs[k].seq; }
+ //for (int k = 0; k < seqs.size(); k++) { delete seqs[k].seq; }
}
//cout << endl << endl;
return 0;
//***************************************************************************************************************
void ChimeraSlayer::printBlock(data_struct data, string flag, ostream& out){
try {
- out << querySeq->getName() << '\t';
+ out << querySeq.getName() << '\t';
out << data.parentA.getName() << "\t" << data.parentB.getName() << '\t';
out << data.divr_qla_qrb << '\t' << data.qla_qrb << '\t' << data.bsa << '\t';
try {
if ((leftChimeric) && (!rightChimeric)) { //print left
- out << querySeq->getName() << '\t';
+ out << querySeq.getName() << '\t';
out << leftdata.results[0].parentA.getName() << "\t" << leftdata.results[0].parentB.getName() << '\t';
out << leftdata.results[0].divr_qla_qrb << '\t' << leftdata.results[0].qla_qrb << '\t' << leftdata.results[0].bsa << '\t';
out << flag << '\t' << leftdata.results[0].winLStart << "-" << leftdata.results[0].winLEnd << '\t' << leftdata.results[0].winRStart << "-" << leftdata.results[0].winREnd << '\t';
}else if ((!leftChimeric) && (rightChimeric)) { //print right
- out << querySeq->getName() << '\t';
+ out << querySeq.getName() << '\t';
out << rightdata.results[0].parentA.getName() << "\t" << rightdata.results[0].parentB.getName() << '\t';
out << rightdata.results[0].divr_qla_qrb << '\t' << rightdata.results[0].qla_qrb << '\t' << rightdata.results[0].bsa << '\t';
}else { //print both results
if (leftdata.flag == "yes") {
- out << querySeq->getName() + "_LEFT" << '\t';
+ out << querySeq.getName() + "_LEFT" << '\t';
out << leftdata.results[0].parentA.getName() << "\t" << leftdata.results[0].parentB.getName() << '\t';
out << leftdata.results[0].divr_qla_qrb << '\t' << leftdata.results[0].qla_qrb << '\t' << leftdata.results[0].bsa << '\t';
if (rightdata.flag == "yes") {
if (leftdata.flag == "yes") { out << endl; }
- out << querySeq->getName() + "_RIGHT"<< '\t';
+ out << querySeq.getName() + "_RIGHT"<< '\t';
out << rightdata.results[0].parentA.getName() << "\t" << rightdata.results[0].parentB.getName() << '\t';
out << rightdata.results[0].divr_qla_qrb << '\t' << rightdata.results[0].qla_qrb << '\t' << rightdata.results[0].bsa << '\t';
string out = "";
if ((leftChimeric) && (!rightChimeric)) { //get left
- out += querySeq->getName() + "\t";
+ out += querySeq.getName() + "\t";
out += leftdata.results[0].parentA.getName() + "\t" + leftdata.results[0].parentB.getName() + "\t";
out += toString(leftdata.results[0].divr_qla_qrb) + "\t" + toString(leftdata.results[0].qla_qrb) + "\t" + toString(leftdata.results[0].bsa) + "\t";
out += flag + "\t" + toString(leftdata.results[0].winLStart) + "-" + toString(leftdata.results[0].winLEnd) + "\t" + toString(leftdata.results[0].winRStart) + "-" + toString(leftdata.results[0].winREnd) + "\t";
}else if ((!leftChimeric) && (rightChimeric)) { //print right
- out += querySeq->getName() + "\t";
+ out += querySeq.getName() + "\t";
out += rightdata.results[0].parentA.getName() + "\t" + rightdata.results[0].parentB.getName() + "\t";
out += toString(rightdata.results[0].divr_qla_qrb) + "\t" + toString(rightdata.results[0].qla_qrb) + "\t" + toString(rightdata.results[0].bsa) + "\t";
}else { //print both results
if (leftdata.flag == "yes") {
- out += querySeq->getName() + "_LEFT\t";
+ out += querySeq.getName() + "_LEFT\t";
out += leftdata.results[0].parentA.getName() + "\t" + leftdata.results[0].parentB.getName() + "\t";
out += toString(leftdata.results[0].divr_qla_qrb) + "\t" + toString(leftdata.results[0].qla_qrb) + "\t" + toString(leftdata.results[0].bsa) + "\t";
if (rightdata.flag == "yes") {
if (leftdata.flag == "yes") { out += "\n"; }
- out += querySeq->getName() + "_RIGHT\t";
+ out += querySeq.getName() + "_RIGHT\t";
out += rightdata.results[0].parentA.getName() + "\t" + rightdata.results[0].parentB.getName() + "\t";
out += toString(rightdata.results[0].divr_qla_qrb) + "\t" + toString(rightdata.results[0].qla_qrb) + "\t" + toString(rightdata.results[0].bsa) + "\t";
string outputString = "";
- outputString += querySeq->getName() + "\t";
+ outputString += querySeq.getName() + "\t";
outputString += data.parentA.getName() + "\t" + data.parentB.getName() + "\t";
outputString += toString(data.divr_qla_qrb) + "\t" + toString(data.qla_qrb) + "\t" + toString(data.bsa) + "\t";
}
}
//***************************************************************************************************************
-vector<Sequence*> ChimeraSlayer::getRefSeqs(Sequence* q, vector<Sequence*>& thisTemplate, vector<Sequence*>& thisFilteredTemplate){
+vector<Sequence> ChimeraSlayer::getRefSeqs(Sequence q, vector<Sequence*>& thisTemplate, vector<Sequence*>& thisFilteredTemplate){
try {
- vector<Sequence*> refSeqs;
+ vector<Sequence> refSeqs;
if (searchMethod == "distance") {
//find closest seqs to query in template - returns copies of seqs so trim does not destroy - remember to deallocate
- Sequence* newSeq = new Sequence(q->getName(), q->getAligned());
+ Sequence* newSeq = new Sequence(q.getName(), q.getAligned());
runFilter(newSeq);
- refSeqs = decalc->findClosest(newSeq, thisTemplate, thisFilteredTemplate, numWanted, minSim);
+ refSeqs = decalc.findClosest(*newSeq, thisTemplate, thisFilteredTemplate, numWanted, minSim);
delete newSeq;
}else if (searchMethod == "blast") {
refSeqs = getBlastSeqs(q, thisTemplate, numWanted); //fills indexes
}
}
//***************************************************************************************************************/
-vector<Sequence*> ChimeraSlayer::getBlastSeqs(Sequence* q, vector<Sequence*>& db, int num) {
+vector<Sequence> ChimeraSlayer::getBlastSeqs(Sequence q, vector<Sequence*>& db, int num) {
try {
- vector<Sequence*> refResults;
+ vector<Sequence> refResults;
//get parts of query
- string queryUnAligned = q->getUnaligned();
+ 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
//cout << "whole length = " << queryUnAligned.length() << '\t' << "left length = " << leftQuery.length() << '\t' << "right length = "<< rightQuery.length() << endl;
- Sequence* queryLeft = new Sequence(q->getName(), leftQuery);
- Sequence* queryRight = new Sequence(q->getName(), rightQuery);
+ Sequence* queryLeft = new Sequence(q.getName(), leftQuery);
+ Sequence* queryRight = new Sequence(q.getName(), rightQuery);
vector<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1, minSim);
vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1, minSim);
for (int i = 0; i < mergedResults.size(); i++) {
//cout << q->getName() << mergedResults[i] << '\t' << db[mergedResults[i]]->getName() << endl;
- if (db[mergedResults[i]]->getName() != q->getName()) {
- Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
+ if (db[mergedResults[i]]->getName() != q.getName()) {
+ Sequence temp(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
refResults.push_back(temp);
-
}
}
//cout << endl << endl;
delete queryRight;
delete queryLeft;
- if (refResults.size() == 0) { m->mothurOut("[WARNING]: mothur found 0 potential parents, so we are not able to check " + q->getName() + ". This could be due to formatdb.exe not being setup properly, please check formatdb.log for errors."); m->mothurOutEndLine(); }
+ if (refResults.size() == 0) { m->mothurOut("[WARNING]: megablast found 0 potential parents, so we are not able to check " + q.getName() + ". This could be due to formatdb.exe not being setup properly, please check formatdb.log for errors."); m->mothurOutEndLine(); }
return refResults;
}
}
}
//***************************************************************************************************************
-vector<Sequence*> ChimeraSlayer::getKmerSeqs(Sequence* q, vector<Sequence*>& db, int num) {
+vector<Sequence> ChimeraSlayer::getKmerSeqs(Sequence q, vector<Sequence*>& db, int num) {
try {
- vector<Sequence*> refResults;
+ vector<Sequence> refResults;
//get parts of query
- string queryUnAligned = q->getUnaligned();
+ 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);
+ Sequence* queryLeft = new Sequence(q.getName(), leftQuery);
+ Sequence* queryRight = new Sequence(q.getName(), rightQuery);
vector<int> tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, num);
vector<int> tempIndexesRight = databaseRight->findClosestSequences(queryRight, num);
for (int i = 0; i < mergedResults.size(); i++) {
//cout << mergedResults[i] << '\t' << db[mergedResults[i]]->getName() << endl;
- if (db[mergedResults[i]]->getName() != q->getName()) {
- Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
+ if (db[mergedResults[i]]->getName() != q.getName()) {
+ Sequence temp(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
refResults.push_back(temp);
}
#include "maligner.h"
#include "slayer.h"
+
+
//***********************************************************************/
//This class was modeled after the chimeraSlayer written by the Broad Institute
/***********************************************************************/
~ChimeraSlayer();
int getChimeras(Sequence*);
- Sequence* print(ostream&, ostream&);
- Sequence* print(ostream&, ostream&, data_results, data_results);
+ Sequence print(ostream&, ostream&);
+ Sequence print(ostream&, ostream&, data_results, data_results);
void printHeader(ostream&);
int doPrep();
data_results getResults() { return printResults; }
#ifdef USE_MPI
- Sequence* print(MPI_File&, MPI_File&);
- Sequence* print(MPI_File&, MPI_File&, data_results, data_results);
+ Sequence print(MPI_File&, MPI_File&);
+ Sequence print(MPI_File&, MPI_File&, data_results, data_results);
#endif
private:
- Sequence* querySeq;
+ Sequence querySeq;
Sequence trimQuery;
- DeCalculator* decalc;
+ DeCalculator decalc;
Database* databaseRight;
Database* databaseLeft;
map<string, int> priority; //for template=self, seqname, seqAligned, abundance
string getBlock(data_struct, string);
string getBlock(data_results, data_results, bool, bool, string);
//int readNameFile(string);
- vector<Sequence*> getTemplate(Sequence*, vector<Sequence*>&);
- vector<Sequence*> getRefSeqs(Sequence*, vector<Sequence*>&, vector<Sequence*>&);
- vector<Sequence*> getBlastSeqs(Sequence*, vector<Sequence*>&, int);
- vector<Sequence*> getKmerSeqs(Sequence*, vector<Sequence*>&, int);
+ vector<Sequence*> getTemplate(Sequence, vector<Sequence*>&);
+ vector<Sequence> getRefSeqs(Sequence, vector<Sequence*>&, vector<Sequence*>&);
+ vector<Sequence> getBlastSeqs(Sequence, vector<Sequence*>&, int);
+ vector<Sequence> getKmerSeqs(Sequence, vector<Sequence*>&, int);
};
data_results rightResults = chimera->getResults();
//if either piece is chimeric then report
- Sequence* trimmed = chimera->print(out, out2, leftResults, rightResults);
- if (trim) { trimmed->printSequence(out3); delete trimmed; }
+ Sequence trimmed = chimera->print(out, out2, leftResults, rightResults);
+ if (trim) { trimmed.printSequence(out3); }
delete left; delete right;
}else { //already chimeric
//print results
- Sequence* trimmed = chimera->print(out, out2);
- if (trim) { trimmed->printSequence(out3); delete trimmed; }
+ Sequence trimmed = chimera->print(out, out2);
+ if (trim) { trimmed.printSequence(out3); }
}
data_results rightResults = chimera->getResults();
//if either piece is chimeric then report
- Sequence* trimmed = chimera->print(outMPI, outAccMPI, leftResults, rightResults);
+ Sequence trimmed = chimera->print(outMPI, outAccMPI, leftResults, rightResults);
if (trim) {
- string outputString = ">" + trimmed->getName() + "\n" + trimmed->getAligned() + "\n";
- delete trimmed;
+ string outputString = ">" + trimmed.getName() + "\n" + trimmed.getAligned() + "\n";
//write to accnos file
int length = outputString.length();
}else {
//print results
- Sequence* trimmed = chimera->print(outMPI, outAccMPI);
+ Sequence trimmed = chimera->print(outMPI, outAccMPI);
if (trim) {
- string outputString = ">" + trimmed->getName() + "\n" + trimmed->getAligned() + "\n";
- delete trimmed;
+ string outputString = ">" + trimmed.getName() + "\n" + trimmed.getAligned() + "\n";
//write to accnos file
int length = outputString.length();
}
//***************************************************************************************************************
//gets closest matches to each end, since chimeras will most likely have different parents on each end
-vector<Sequence*> DeCalculator::findClosest(Sequence* querySeq, vector<Sequence*>& thisTemplate, vector<Sequence*>& thisFilteredTemplate, int numWanted, int minSim) {
+vector<Sequence> DeCalculator::findClosest(Sequence querySeq, vector<Sequence*>& thisTemplate, vector<Sequence*>& thisFilteredTemplate, int numWanted, int minSim) {
try {
//indexes.clear();
- vector<Sequence*> seqsMatches;
+ vector<Sequence> seqsMatches;
vector<SeqDist> distsLeft;
vector<SeqDist> distsRight;
Dist* distcalculator = new eachGapDist();
- string queryUnAligned = querySeq->getUnaligned();
+ string queryUnAligned = querySeq.getUnaligned();
int numBases = int(queryUnAligned.length() * 0.33);
string leftQuery = ""; //first 1/3 of the sequence
string rightQuery = ""; //last 1/3 of the sequence
- string queryAligned = querySeq->getAligned();
+ string queryAligned = querySeq.getAligned();
//left side
bool foundFirstBase = false;
}
rightQuery = queryAligned.substr(rightSpot, (lastBaseSpot-rightSpot+1)); //sequence from pos spot to end
- Sequence queryLeft(querySeq->getName(), leftQuery);
- Sequence queryRight(querySeq->getName(), rightQuery);
+ Sequence queryLeft(querySeq.getName(), leftQuery);
+ Sequence queryRight(querySeq.getName(), rightQuery);
//cout << querySeq->getName() << '\t' << leftSpot << '\t' << rightSpot << '\t' << firstBaseSpot << '\t' << lastBaseSpot << endl;
//cout << queryUnAligned.length() << '\t' << queryLeft.getUnaligned().length() << '\t' << queryRight.getUnaligned().length() << endl;
for (int i = 0; i < dists.size(); i++) {
// cout << db[dists[i].index]->getName() << '\t' << dists[i].dist << endl;
- if ((thisTemplate[dists[i].index]->getName() != querySeq->getName()) && (((1.0-dists[i].dist)*100) >= minSim)) {
- Sequence* temp = new Sequence(thisTemplate[dists[i].index]->getName(), thisTemplate[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother.
+ if ((thisTemplate[dists[i].index]->getName() != querySeq.getName()) && (((1.0-dists[i].dist)*100) >= minSim)) {
+ Sequence temp(thisTemplate[dists[i].index]->getName(), thisTemplate[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother.
//cout << querySeq->getName() << '\t' << thisTemplate[dists[i].index]->getName() << '\t' << dists[i].dist << endl;
seqsMatches.push_back(temp);
}
}
}
/***************************************************************************************************************/
-map<int, int> DeCalculator::trimSeqs(Sequence* query, vector<Sequence*> topMatches) {
+map<int, int> DeCalculator::trimSeqs(Sequence& query, vector<Sequence>& topMatches) {
try {
int frontPos = 0; //should contain first position in all seqs that is not a gap character
- int rearPos = query->getAligned().length();
+ int rearPos = query.getAligned().length();
//********find first position in topMatches that is a non gap character***********//
//find first position all query seqs that is a non gap character
for (int i = 0; i < topMatches.size(); i++) {
- string aligned = topMatches[i]->getAligned();
+ string aligned = topMatches[i].getAligned();
int pos = 0;
//find first spot in this seq
}
- string aligned = query->getAligned();
+ string aligned = query.getAligned();
int pos = 0;
//find first position in query that is a non gap character
//********find last position in topMatches that is a non gap character***********//
for (int i = 0; i < topMatches.size(); i++) {
- string aligned = topMatches[i]->getAligned();
+ string aligned = topMatches[i].getAligned();
int pos = aligned.length();
//find first spot in this seq
}
- aligned = query->getAligned();
+ aligned = query.getAligned();
pos = aligned.length();
//find last position in query that is a non gap character
map<int, int> trimmedPos;
//check to make sure that is not whole seq
if ((rearPos - frontPos - 1) <= 0) {
- query->setAligned("");
+ query.setAligned("");
//trim topMatches
for (int i = 0; i < topMatches.size(); i++) {
- topMatches[i]->setAligned("");
+ topMatches[i].setAligned("");
}
}else {
//trim query
- string newAligned = query->getAligned();
+ string newAligned = query.getAligned();
newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
- query->setAligned(newAligned);
+ query.setAligned(newAligned);
//trim topMatches
for (int i = 0; i < topMatches.size(); i++) {
- newAligned = topMatches[i]->getAligned();
+ newAligned = topMatches[i].getAligned();
newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
- topMatches[i]->setAligned(newAligned);
+ topMatches[i].setAligned(newAligned);
}
for (int i = 0; i < newAligned.length(); i++) {
DeCalculator() { m = MothurOut::getInstance(); }
~DeCalculator() {};
- vector<Sequence*> findClosest(Sequence*, vector<Sequence*>&, vector<Sequence*>&, int, int); //takes querySeq, a reference db, filteredRefDB, numWanted, minSim
+ vector<Sequence> findClosest(Sequence, vector<Sequence*>&, vector<Sequence*>&, int, int); //takes querySeq, a reference db, filteredRefDB, numWanted, minSim
Sequence* findClosest(Sequence*, vector<Sequence*>);
set<int> getPos() { return h; }
void setMask(string);
void setAlignmentLength(int l) { alignLength = l; }
void runMask(Sequence*);
void trimSeqs(Sequence*, Sequence*, map<int, int>&);
- map<int, int> trimSeqs(Sequence*, vector<Sequence*>);
+ map<int, int> trimSeqs(Sequence&, vector<Sequence>&);
void removeObviousOutliers(vector< vector<float> >&, int);
vector<float> calcFreq(vector<Sequence*>, string);
vector<int> findWindows(Sequence*, int, int, int&, int);
#include "maligner.h"
/***********************************************************************/ //int num, int match, int misMatch, , string mode, Database* dataLeft, Database* dataRight
-Maligner::Maligner(vector<Sequence*> temp, int match, int misMatch, float div, int ms, int minCov) : db(temp), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov) {
+Maligner::Maligner(vector<Sequence> temp, int match, int misMatch, float div, int ms, int minCov) : db(temp), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov) {
//numWanted(num), , searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight)
m = MothurOut::getInstance();
}
/***********************************************************************/
-string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
+string Maligner::getResults(Sequence q, DeCalculator decalc) {
try {
outputResults.clear();
//make copy so trimming doesn't destroy query from calling class - remember to deallocate
- query = new Sequence(q->getName(), q->getAligned());
+ query.setName(q.getName()); query.setAligned(q.getAligned());
string chimera;
//copy refSeqs so that filter does not effect original
for(int i = 0; i < db.size(); i++) {
- Sequence* newSeq = new Sequence(db[i]->getName(), db[i]->getAligned());
+ Sequence newSeq(db[i].getName(), db[i].getAligned());
refSeqs.push_back(newSeq);
}
refSeqs = minCoverageFilter(refSeqs);
if (refSeqs.size() < 2) {
- for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
+ //for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
percentIdenticalQueryChimera = 0.0;
return "unknown";
}
if (m->control_pressed) { return chimera; }
//free memory
- delete query;
+ //delete query;
- for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
+ //for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; }
return chimera;
}
}
}
/***********************************************************************/
-string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
+string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator decalc) {
try {
string chimera;
-
//trims seqs to first non gap char in all seqs and last non gap char in all seqs
- spotMap = decalc->trimSeqs(query, refSeqs);
-
+ spotMap = decalc.trimSeqs(query, refSeqs);
+
//you trimmed the whole sequence, skip
- if (query->getAligned() == "") { return "no"; }
+ if (query.getAligned() == "") { return "no"; }
- vector<Sequence*> temp = refSeqs;
-
-// for(int i=0;i<refSeqs.size();i++){
-// cout << refSeqs[i]->getName() << endl;
-// }
-
+ vector<Sequence> temp = refSeqs;
temp.push_back(query);
- verticalFilter(temp);
+ temp = verticalFilter(temp);
+ query = temp[temp.size()-1];
+ for (int i = 0; i < temp.size()-1; i++) { refSeqs[i] = temp[i]; }
//for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl ; }//<< refSeqs[i]->getAligned() << endl
- vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
+ vector< vector<score_struct> > matrix = buildScoreMatrix(query.getAligned().length(), refSeqs.size()); //builds and initializes
if (m->control_pressed) { return chimera; }
if (m->control_pressed) { return chimera; }
- vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
+ vector<trace_struct> trace = mapTraceRegionsToAlignment(path);
if (trace.size() > 1) { chimera = "yes"; }
else { chimera = "no"; return chimera; }
int traceStart = path[0].col;
int traceEnd = path[path.size()-1].col;
- string queryInRange = query->getAligned();
+ string queryInRange = query.getAligned();
queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
// cout << queryInRange << endl;
string chimeraSeq = constructChimericSeq(trace, refSeqs);
// cout << regionStart << '\t' << regionEnd << '\t' << seqIndex << endl;
results temp;
- temp.parent = refSeqs[seqIndex]->getName();
- temp.parentAligned = db[seqIndex]->getAligned();
+ temp.parent = refSeqs[seqIndex].getName();
+ temp.parentAligned = db[seqIndex].getAligned();
temp.nastRegionStart = spotMap[regionStart];
temp.nastRegionEnd = spotMap[regionEnd];
temp.regionStart = unalignedMap[regionStart];
temp.regionEnd = unalignedMap[regionEnd];
- string parentInRange = refSeqs[seqIndex]->getAligned();
+ string parentInRange = refSeqs[seqIndex].getAligned();
parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
temp.queryToParent = computePercentID(queryInRange, parentInRange);
temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
- string queryInRegion = query->getAligned();
+ string queryInRegion = query.getAligned();
queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
- string parentInRegion = refSeqs[seqIndex]->getAligned();
+ string parentInRegion = refSeqs[seqIndex].getAligned();
parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
}
/***********************************************************************/
//removes top matches that do not have minimum coverage with query.
-vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
+vector<Sequence> Maligner::minCoverageFilter(vector<Sequence> ref){
try {
- vector<Sequence*> newRefs;
+ vector<Sequence> newRefs;
- string queryAligned = query->getAligned();
+ string queryAligned = query.getAligned();
for (int i = 0; i < ref.size(); i++) {
- string refAligned = ref[i]->getAligned();
+ string refAligned = ref[i].getAligned();
int numBases = 0;
int numCovered = 0;
//if coverage above minimum
if (coverage > minCoverage) {
newRefs.push_back(ref[i]);
- }else {
- delete ref[i];
- }
+ }//else {
+ //delete ref[i];
+ //}
}
return newRefs;
int Maligner::computeChimeraPenalty() {
try {
- int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
+ int numAllowable = ((1.0 - (1.0/minDivR)) * query.getNumBases());
// if(numAllowable < 1){ numAllowable = 1; }
}
/***********************************************************************/
//this is a vertical filter
-void Maligner::verticalFilter(vector<Sequence*> seqs) {
+vector<Sequence> Maligner::verticalFilter(vector<Sequence> seqs) {
try {
- vector<int> gaps; gaps.resize(query->getAligned().length(), 0);
+ vector<int> gaps; gaps.resize(query.getAligned().length(), 0);
- string filterString = (string(query->getAligned().length(), '1'));
+ string filterString = (string(query.getAligned().length(), '1'));
//for each sequence
for (int i = 0; i < seqs.size(); i++) {
- string seqAligned = seqs[i]->getAligned();
+ string seqAligned = seqs[i].getAligned();
for (int j = 0; j < seqAligned.length(); j++) {
//if this spot is a gap
//zero out spot where all sequences have blanks
int numColRemoved = 0;
- for(int i = 0; i < seqs[0]->getAligned().length(); i++){
+ for(int i = 0; i < seqs[0].getAligned().length(); i++){
if(gaps[i] == seqs.size()) { filterString[i] = '0'; numColRemoved++; }
}
//for each sequence
for (int i = 0; i < seqs.size(); i++) {
- string seqAligned = seqs[i]->getAligned();
+ string seqAligned = seqs[i].getAligned();
string newAligned = "";
int count = 0;
}
}
- seqs[i]->setAligned(newAligned);
+ seqs[i].setAligned(newAligned);
}
- string query = seqs[seqs.size()-1]->getAligned();
+ string query = seqs[seqs.size()-1].getAligned();
int queryLength = query.length();
unalignedMap.resize(queryLength, 0);
}
spotMap = newMap;
+
+ return seqs;
}
catch(exception& e) {
m->errorOut(e, "Maligner", "verticalFilter");
//***************************************************************************************************************
-void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
+void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence> seqs, int penalty) {
try{
//get matrix dimensions
- int numCols = query->getAligned().length();
+ int numCols = query.getAligned().length();
int numRows = seqs.size();
// cout << numRows << endl;
//initialize first col
- string queryAligned = query->getAligned();
+ string queryAligned = query.getAligned();
for (int i = 0; i < numRows; i++) {
- string subjectAligned = seqs[i]->getAligned();
+ string subjectAligned = seqs[i].getAligned();
//are you both gaps?
if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
// for (int i = 0; i < 1; i++) { //iterate through matrix rows
for (int i = 0; i < numRows; i++) { //iterate through matrix rows
- string subjectAligned = seqs[i]->getAligned();
+ string subjectAligned = seqs[i].getAligned();
int matchMisMatchScore = 0;
//are you both gaps?
try {
//get matrix dimensions
- int numCols = query->getAligned().length();
+ int numCols = query.getAligned().length();
int numRows = ms.size();
}
}
//***************************************************************************************************************
-vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
+vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path) {
try {
vector<trace_struct> trace;
*/
//***************************************************************************************************************
-string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
+string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence> seqs) {
try {
string chimera = "";
for (int i = 0; i < trace.size(); i++) {
// cout << i << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
- string seqAlign = seqs[trace[i].row]->getAligned();
+ string seqAlign = seqs[trace[i].row].getAligned();
seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
chimera += seqAlign;
}
//***************************************************************************************************************
-string Maligner::constructAntiChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
+string Maligner::constructAntiChimericSeq(vector<trace_struct> trace, vector<Sequence> seqs) {
try {
string antiChimera = "";
int oppositeIndex = trace.size() - i - 1;
- string seqAlign = seqs[trace[oppositeIndex].row]->getAligned();
+ string seqAlign = seqs[trace[oppositeIndex].row].getAligned();
seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
antiChimera += seqAlign;
}
public:
- Maligner(vector<Sequence*>, int, int, float, int, int); //int, int, int, , string, Database*, Database*
+ Maligner(vector<Sequence>, int, int, float, int, int); //int, int, int, , string, Database*, Database*
~Maligner() {};
- string getResults(Sequence*, DeCalculator*);
+ string getResults(Sequence, DeCalculator);
float getPercentID() { return percentIdenticalQueryChimera; }
vector<results> getOutput() { return outputResults; }
private:
- Sequence* query;
- vector<Sequence*> refSeqs;
- vector<Sequence*> db;
+ Sequence query;
+ vector<Sequence> refSeqs;
+ vector<Sequence> db;
int minCoverage, minSimilarity, matchScore, misMatchPenalty;
float minDivR, percentIdenticalQueryChimera;
vector<results> outputResults;
map<int, int> spotMap;
vector<int> unalignedMap;
- vector<Sequence*> minCoverageFilter(vector<Sequence*>); //removes top matches that do not have minimum coverage with query.
+ vector<Sequence> minCoverageFilter(vector<Sequence>); //removes top matches that do not have minimum coverage with query.
int computeChimeraPenalty();
- void verticalFilter(vector<Sequence*>);
+ vector<Sequence> verticalFilter(vector<Sequence>);
vector< vector<score_struct> > buildScoreMatrix(int, int);
- void fillScoreMatrix(vector<vector<score_struct> >&, vector<Sequence*>, int);
+ void fillScoreMatrix(vector<vector<score_struct> >&, vector<Sequence>, int);
vector<score_struct> extractHighestPath(vector<vector<score_struct> >);
- vector<trace_struct> mapTraceRegionsToAlignment(vector<score_struct>, vector<Sequence*>);
- string constructChimericSeq(vector<trace_struct>, vector<Sequence*>);
- string constructAntiChimericSeq(vector<trace_struct>, vector<Sequence*>);
+ vector<trace_struct> mapTraceRegionsToAlignment(vector<score_struct>);
+ string constructChimericSeq(vector<trace_struct>, vector<Sequence>);
+ string constructAntiChimericSeq(vector<trace_struct>, vector<Sequence>);
float computePercentID(string, string);
- string chimeraMaligner(int, DeCalculator*);
+ string chimeraMaligner(int, DeCalculator);
MothurOut* m;
};
}
}
//***************************************************************************************************************
-Sequence* Pintail::print(ostream& out, ostream& outAcc) {
+Sequence Pintail::print(ostream& out, ostream& outAcc) {
try {
int index = ceil(deviation);
for (int m = 0; m < expectedDistance.size(); m++) { out << expectedDistance[m] << '\t'; }
out << endl;
- return NULL;
+ return *querySeq;
}
catch(exception& e) {
~Pintail();
int getChimeras(Sequence*);
- Sequence* print(ostream&, ostream&);
+ Sequence print(ostream&, ostream&);
void setCons(string c) { consfile = c; }
void setQuantiles(string q) { quanfile = q; }
#ifdef USE_MPI
- Sequence* print(MPI_File&, MPI_File&);
+ Sequence print(MPI_File&, MPI_File&);
#endif
private:
Slayer::Slayer(int win, int increment, int parentThreshold, float div, int i, int snp, int mi) :
minBS(mi), windowSize(win), windowStep(increment), parentFragmentThreshold(parentThreshold), divRThreshold(div), iters(i), percentSNPSample(snp){ m = MothurOut::getInstance(); }
/***********************************************************************/
-string Slayer::getResults(Sequence* query, vector<Sequence*> refSeqs) {
+string Slayer::getResults(Sequence query, vector<Sequence> refSeqs) {
try {
vector<data_struct> all; all.clear();
- myQuery = *query;
+ myQuery = query;
for (int i = 0; i < refSeqs.size(); i++) {
if (m->control_pressed) { return "no"; }
//make copies of query and each parent because runBellerophon removes gaps and messes them up
- Sequence* q = new Sequence(query->getName(), query->getAligned());
- Sequence* leftParent = new Sequence(refSeqs[i]->getName(), refSeqs[i]->getAligned());
- Sequence* rightParent = new Sequence(refSeqs[j]->getName(), refSeqs[j]->getAligned());
+ Sequence q(query.getName(), query.getAligned());
+ Sequence leftParent(refSeqs[i].getName(), refSeqs[i].getAligned());
+ Sequence rightParent(refSeqs[j].getName(), refSeqs[j].getAligned());
//cout << q->getName() << endl << q->getAligned() << endl << endl;
- //cout << leftParent->getName() << endl << leftParent->getAligned() << endl << endl;
+ //cout << leftParent.getName() << '\t' << leftParent.getAligned().length() << endl << endl;
+ //cout << rightParent.getName() << '\t' << rightParent.getAligned().length() << endl << endl;
+ //cout << q.getName() << '\t' << q.getAligned().length() << endl << endl;
//cout << rightParent->getName() << endl << rightParent->getAligned() << endl << endl;
//cout << " length = " << rightParent->getAligned().length() << endl;
map<int, int> spots; //map from spot in original sequence to spot in filtered sequence for query and both parents
vector<data_struct> divs = runBellerophon(q, leftParent, rightParent, spots);
- if (m->control_pressed) { delete q; delete leftParent; delete rightParent; return "no"; }
+ if (m->control_pressed) { return "no"; }
// cout << "examining:\t" << refSeqs[i]->getName() << '\t' << refSeqs[j]->getName() << endl;
vector<data_struct> selectedDivs;
for (int k = 0; k < divs.size(); k++) {
vector<snps> snpsLeft = getSNPS(divs[k].parentA.getAligned(), divs[k].querySeq.getAligned(), divs[k].parentB.getAligned(), divs[k].winLStart, divs[k].winLEnd);
vector<snps> snpsRight = getSNPS(divs[k].parentA.getAligned(), divs[k].querySeq.getAligned(), divs[k].parentB.getAligned(), divs[k].winRStart, divs[k].winREnd);
- if (m->control_pressed) { delete q; delete leftParent; delete rightParent; return "no"; }
+ if (m->control_pressed) { return "no"; }
int numSNPSLeft = snpsLeft.size();
int numSNPSRight = snpsRight.size();
float BS_A, BS_B;
bootstrapSNPS(snpsLeft, snpsRight, BS_A, BS_B, iters);
- if (m->control_pressed) { delete q; delete leftParent; delete rightParent; return "no"; }
+ if (m->control_pressed) { return "no"; }
divs[k].bsa = BS_A;
divs[k].bsb = BS_B;
//save selected
for (int mi = 0; mi < selectedDivs.size(); mi++) { all.push_back(selectedDivs[mi]); }
-
- delete q;
- delete leftParent;
- delete rightParent;
}
}
}
}
/***********************************************************************/
-vector<data_struct> Slayer::runBellerophon(Sequence* q, Sequence* pA, Sequence* pB, map<int, int>& spots) {
+vector<data_struct> Slayer::runBellerophon(Sequence q, Sequence pA, Sequence pB, map<int, int>& spots) {
try{
vector<data_struct> data;
//vertical filter
- vector<Sequence*> temp;
- temp.push_back(q); temp.push_back(pA); temp.push_back(pB);
+ //cout << q.getName() << endl << q.getAligned() << endl << endl;
+ //cout << pA.getName() << endl << pA.getUnaligned() << endl << endl;
+ //cout << pB.getName() << endl << pB.getUnaligned() << endl << endl;
//maps spot in new alignment to spot in alignment before filter
- spots = verticalFilter(temp); //fills baseSpots
+ spots = verticalFilter(q, pA, pB); //fills baseSpots
//get these to avoid numerous function calls
- string query = q->getAligned();
- string parentA = pA->getAligned();
- string parentB = pB->getAligned();
+ string query = q.getAligned();
+ string parentA = pA.getAligned();
+ string parentB = pB.getAligned();
int length = query.length();
-//cout << q->getName() << endl << q->getAligned() << endl << endl;
-//cout << pA->getName() << endl << pA->getUnaligned() << endl << endl;
-//cout << pB->getName() << endl << pB->getUnaligned() << endl << endl;
+//cout << q.getName() << endl << q.getAligned() << endl << endl;
+//cout << pA.getName() << endl << pA.getUnaligned() << endl << endl;
+//cout << pB.getName() << endl << pB.getUnaligned() << endl << endl;
//cout << " length = " << length << endl;
//check window size
member.winLEnd = breakpoint;
member.winRStart = breakpoint+1;
member.winREnd = length-1;
- member.querySeq = *(q);
- member.parentA = *(pA);
- member.parentB = *(pB);
+ member.querySeq = q;
+ member.parentA = pA;
+ member.parentB = pB;
member.bsa = 0;
member.bsb = 0;
member.bsMax = 0;
}
/***********************************************************************/
//remove columns that contain any gaps
-map<int, int> Slayer::verticalFilter(vector<Sequence*> seqs) {
+map<int, int> Slayer::verticalFilter(Sequence& q, Sequence& pA, Sequence& pB) {
try {
//find baseSpots
baseSpots.clear();
baseSpots.resize(3); //query, parentA, parentB
- vector<int> gaps; gaps.resize(seqs[0]->getAligned().length(), 0);
+ vector<int> gaps; gaps.resize(q.getAligned().length(), 0);
- string filterString = (string(seqs[0]->getAligned().length(), '1'));
+ string filterString = (string(q.getAligned().length(), '1'));
- //for each sequence
- for (int i = 0; i < seqs.size(); i++) {
+ string seqAligned = q.getAligned();
+ for (int j = 0; j < seqAligned.length(); j++) {
+ //if this spot is a gap
+ if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N')) { gaps[j]++; }
+ }
- string seqAligned = seqs[i]->getAligned();
-
- for (int j = 0; j < seqAligned.length(); j++) {
- //if this spot is a gap
- if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N')) { gaps[j]++; }
- }
+ seqAligned = pA.getAligned();
+ for (int j = 0; j < seqAligned.length(); j++) {
+ //if this spot is a gap
+ if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N')) { gaps[j]++; }
}
+ seqAligned = pB.getAligned();
+ for (int j = 0; j < seqAligned.length(); j++) {
+ //if this spot is a gap
+ if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N')) { gaps[j]++; }
+ }
+
+
//zero out spot where any sequences have blanks
int numColRemoved = 0;
int count = 0;
map<int, int> maskMap; maskMap.clear();
- for(int i = 0; i < seqs[0]->getAligned().length(); i++){
+ for(int i = 0; i < q.getAligned().length(); i++){
if(gaps[i] != 0) { filterString[i] = '0'; numColRemoved++; }
else {
maskMap[count] = i;
}
}
- //for each sequence
- for (int i = 0; i < seqs.size(); i++) {
-
- string seqAligned = seqs[i]->getAligned();
- string newAligned = "";
+ seqAligned = q.getAligned();
+ string newAligned = "";
- int baseCount = 0;
- int count = 0;
- for (int j = 0; j < seqAligned.length(); j++) {
- //are you a base
- if ((seqAligned[j] != '-') && (seqAligned[j] != '.') && (toupper(seqAligned[j]) != 'N')) { baseCount++; }
+ int baseCount = 0;
+ count = 0;
+ for (int j = 0; j < seqAligned.length(); j++) {
+ //are you a base
+ if ((seqAligned[j] != '-') && (seqAligned[j] != '.') && (toupper(seqAligned[j]) != 'N')) { baseCount++; }
- //if this spot is not a gap
- if (filterString[j] == '1') {
- newAligned += seqAligned[j];
- baseSpots[i][count] = baseCount;
- count++;
- }
+ //if this spot is not a gap
+ if (filterString[j] == '1') {
+ newAligned += seqAligned[j];
+ baseSpots[0][count] = baseCount;
+ count++;
}
+ }
- seqs[i]->setAligned(newAligned);
+ q.setAligned(newAligned);
+
+ seqAligned = pA.getAligned();
+ newAligned = "";
+
+ baseCount = 0;
+ count = 0;
+ for (int j = 0; j < seqAligned.length(); j++) {
+ //are you a base
+ if ((seqAligned[j] != '-') && (seqAligned[j] != '.') && (toupper(seqAligned[j]) != 'N')) { baseCount++; }
+
+ //if this spot is not a gap
+ if (filterString[j] == '1') {
+ newAligned += seqAligned[j];
+ baseSpots[1][count] = baseCount;
+ count++;
+ }
+ }
+
+ pA.setAligned(newAligned);
+
+ seqAligned = pB.getAligned();
+ newAligned = "";
+
+ baseCount = 0;
+ count = 0;
+ for (int j = 0; j < seqAligned.length(); j++) {
+ //are you a base
+ if ((seqAligned[j] != '-') && (seqAligned[j] != '.') && (toupper(seqAligned[j]) != 'N')) { baseCount++; }
+
+ //if this spot is not a gap
+ if (filterString[j] == '1') {
+ newAligned += seqAligned[j];
+ baseSpots[2][count] = baseCount;
+ count++;
+ }
}
+ pB.setAligned(newAligned);
+
+
return maskMap;
}
catch(exception& e) {
Slayer(int, int, int, float, int, int, int);
~Slayer() {};
- string getResults(Sequence*, vector<Sequence*>);
+ string getResults(Sequence, vector<Sequence>);
vector<data_struct> getOutput() { return outputResults; }
vector< map<int, int> > baseSpots;
Sequence myQuery;
- map<int, int> verticalFilter(vector<Sequence*>);
+ map<int, int> verticalFilter(Sequence&, Sequence&, Sequence&);
float computePercentID(string, string, int, int);
- vector<data_struct> runBellerophon(Sequence*, Sequence*, Sequence*, map<int, int>&);
+ vector<data_struct> runBellerophon(Sequence, Sequence, Sequence, map<int, int>&);
vector<snps> getSNPS(string, string, string, int, int);
int bootstrapSNPS(vector<snps>, vector<snps>, float&, float&, int);
float snpQA(vector<snps>);
projects / mothur.git / commitdiff