for(int i=alignmentLength-1;i>=0;i--){
flipSeq += queryAlignment[i];
}
-
+
return flipSeq;
}
catch(exception& e) {
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;
+// cout << querySeq->getName() << "\tyes" << endl;
outAccString += querySeq->getName() + "\n";
results = true;
}
//put seqs into vector to send to slayer
+
+// cout << query->getAligned() << endl;
vector<Sequence*> seqsForSlayer;
- for (int k = 0; k < seqs.size(); k++) { seqsForSlayer.push_back(seqs[k].seq); }
+ for (int k = 0; k < seqs.size(); k++) {
+// cout << seqs[k].seq->getAligned() << endl;
+ seqsForSlayer.push_back(seqs[k].seq);
+
+ }
if (m->control_pressed) { for (int k = 0; k < seqs.size(); k++) { delete seqs[k].seq; } return 0; }
}
}
+
+
+// for(int i=0;i<refResults.size();i++){
+// cout << refResults[i]->getName() << endl;
+// }
delete queryRight;
delete queryLeft;
CommandParameter pmincov("mincov", "Number", "", "70", "", "", "",false,false); parameters.push_back(pmincov);
CommandParameter pminsnp("minsnp", "Number", "", "10", "", "", "",false,false); parameters.push_back(pminsnp);
CommandParameter pminbs("minbs", "Number", "", "90", "", "", "",false,false); parameters.push_back(pminbs);
- CommandParameter psearch("search", "Multiple", "kmer-blast-distance", "distance", "", "", "",false,false); parameters.push_back(psearch);
+ CommandParameter psearch("search", "Multiple", "kmer-blast-distance", "blast", "", "", "",false,false); parameters.push_back(psearch);
CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
CommandParameter prealign("realign", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(prealign);
CommandParameter ptrim("trim", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(ptrim);
int traceStart = path[0].col;
int traceEnd = path[path.size()-1].col;
string queryInRange = query->getAligned();
- queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart));
+ queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
string chimeraSeq = constructChimericSeq(trace, refSeqs);
+// cout << queryInRange.length() << endl;
+// cout << chimeraSeq.length() << endl;
percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
temp.regionEnd = regionEnd;
string parentInRange = refSeqs[seqIndex]->getAligned();
- parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart));
+ parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
temp.queryToParent = computePercentID(queryInRange, parentInRange);
temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
string queryInRegion = query->getAligned();
- queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart));
+ queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
string parentInRegion = refSeqs[seqIndex]->getAligned();
- parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart));
+ parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
- //cout << query->getName() << '\t' << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << ", " << temp.divR << endl;
+// cout << query->getName() << '\t' << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << ", " << temp.divR << endl;
outputResults.push_back(temp);
}
seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
chimera += seqAlign;
}
-
- if (chimera != "") { chimera = chimera.substr(0, (chimera.length()-1)); }
+// cout << chimera << endl;
+// if (chimera != "") { chimera = chimera.substr(0, (chimera.length()-1)); } //this was introducing a fence post error
+// cout << chimera << endl;
return chimera;
}
catch(exception& e) {
return -1.0;
}
+// cout << queryAlign.length() << endl;
int numIdentical = 0;
int countA = 0;
int countB = 0;
if ((queryAlign[i] == 'G') || (queryAlign[i] == 'T') || (queryAlign[i] == 'A') || (queryAlign[i] == 'C')) { charA = true; }
if ((chimera[i] == 'G') || (chimera[i] == 'T') || (chimera[i] == 'A') || (chimera[i] == 'C')) { charB = true; }
+
if (charA || charB) {
if (charA) { countA++; }
numIdentical++;
}
}
+// cout << queryAlign[i] << '\t' << chimera[i] << '\t' << countA << '\t' << countB << endl;
+
}
}
+// cout << "pat\t" << countA << '\t' << countB << '\t' << numIdentical << endl;
+
+
float numBases = (countA + countB) /(float) 2;
if (numBases == 0) { return 0; }
float percentIdentical = (numIdentical/(float)numBases) * 100;
-
+
+// cout << percentIdentical << endl;
+
return percentIdentical;
}
selectedDivs.push_back(divs[k]);
}
}
-
+
//save selected
for (int mi = 0; mi < selectedDivs.size(); mi++) { all.push_back(selectedDivs[mi]); }
}
}
+
// compute bootstrap support
if (all.size() > 0) {
//sort them
int breakpoint = i;
int leftLength = breakpoint + 1;
int rightLength = length - leftLength;
-
+
float QLA = computePercentID(query, parentA, 0, breakpoint);
- float QRB = computePercentID(query, parentB, breakpoint+1, length - 1);
+ float QRB = computePercentID(query, parentB, breakpoint+1, length);
float QLB = computePercentID(query, parentB, 0, breakpoint);
- float QRA = computePercentID(query, parentA, breakpoint+1, length - 1);
+ float QRA = computePercentID(query, parentA, breakpoint+1, length);
float LAB = computePercentID(parentA, parentB, 0, breakpoint);
- float RAB = computePercentID(parentA, parentB, breakpoint+1, length - 1);
-
+ float RAB = computePercentID(parentA, parentB, breakpoint+1, length);
+
float AB = ((LAB*leftLength) + (RAB*rightLength)) / (float) length;
float QA = ((QLA*leftLength) + (QRA*rightLength)) / (float) length;
float QB = ((QLB*leftLength) + (QRB*rightLength)) / (float) length;
float divR_QLA_QRB = min((QLA_QRB/QA), (QLA_QRB/QB));
float divR_QLB_QRA = min((QLB_QRA/QA), (QLB_QRA/QB));
+
//cout << q->getName() << '\t';
//cout << pA->getName() << '\t';
//cout << pB->getName() << '\t';
member.winLEnd = breakpoint;
member.winRStart = breakpoint+1;
member.winREnd = length-1;
- member.querySeq = *(q);
+ member.querySeq = *(q);
member.parentA = *(pA);
member.parentB = *(pB);
member.bsa = 0;
try {
vector<snps> data;
-//cout << left << '\t' << right << endl;
+
for (int i = left; i <= right; i++) {
char A = parentA[i];
char B = parentB[i];
if ((A != Q) || (B != Q)) {
-//cout << "not equal " << Q << '\t' << A << '\t' << B << endl;
-
+
//ensure not neighboring a gap. change to 12/09 release of chimeraSlayer - not sure what this adds, but it eliminates alot of SNPS
+
+
if (
//did query loose a base here during filter??
( i == 0 || abs (baseSpots[0][i] - baseSpots[0][i-1]) == 1) &&
- ( i == query.length() || abs (baseSpots[0][i] - baseSpots[0][i+1]) == 1)
+ ( i == query.length()-1 || abs (baseSpots[0][i] - baseSpots[0][i+1]) == 1)
&&
//did parentA loose a base here during filter??
( i == 0 || abs (baseSpots[1][i] - baseSpots[1][i-1]) == 1) &&
- ( i == parentA.length() || abs (baseSpots[1][i] - baseSpots[1][i+1]) == 1)
+ ( i == parentA.length()-1 || abs (baseSpots[1][i] - baseSpots[1][i+1]) == 1)
&&
//did parentB loose a base here during filter??
( i == 0 || abs (baseSpots[2][i] - baseSpots[2][i-1]) == 1) &&
- ( i == parentB.length() || abs (baseSpots[2][i] - baseSpots[2][i+1]) == 1)
+ ( i == parentB.length()-1 || abs (baseSpots[2][i] - baseSpots[2][i+1]) == 1)
)
{
-
snps member;
member.queryChar = Q;
member.parentAChar = A;
member.parentBChar = B;
-//cout << "not neighboring a gap " << Q << '\t' << A << '\t' << B << '\t' << baseSpots[0][i] << '\t' << baseSpots[0][i+1] << '\t' << baseSpots[0][i-1] << '\t' << baseSpots[1][i] << '\t' << baseSpots[1][i+1] << '\t' << baseSpots[1][i-1] << '\t' << baseSpots[2][i] << '\t' << baseSpots[2][i+1] << '\t' << baseSpots[2][i-1] << endl;
data.push_back(member);
}
}
-// cout << i << '\t' << data.size() << endl;
}
return data;
int numLeft = max(1, int(left.size() * percentSNPSample/(float)100 + 0.5));
int numRight = max(1, int(right.size() * percentSNPSample/(float)100 + 0.5));
- //cout << numLeft << '\t' << numRight << endl;
+
for (int i = 0; i < numIters; i++) {
//random sampling with replacement.
}
}
}
+
}
float numBases = (countA + countB) /(float) 2;
if (numBases == 0) { return 0; }
float percentIdentical = (numIdentical/(float)numBases) * 100;
-
+
return percentIdentical;
}