//***************************************************************************************************************
-RefChimeraTest::RefChimeraTest(vector<Sequence>& refs, string chimeraReportFileName){
+RefChimeraTest::RefChimeraTest(vector<Sequence>& refs, bool aligned) : aligned(aligned){
m = MothurOut::getInstance();
- m->openOutputFile(chimeraReportFileName, chimeraReportFile);
numRefSeqs = refs.size();
referenceSeqs.resize(numRefSeqs);
alignLength = referenceSeqs[0].length();
- chimeraReportFile << "queryName\tbestRef\tbestSequenceMismatch\tleftParentChi,rightParentChi\tbreakPointChi\tminMismatchToChimera\tdistToBestMera\tnumParents" << endl;
-// chimeraReportFile << "leftParentTri,middleParentTri,rightParentTri\tbreakPointTriA,breakPointTriB\tminMismatchToTrimera\tdistToBestMera\tnMera" << endl;
-
}
+//***************************************************************************************************************
+int RefChimeraTest::printHeader(ofstream& chimeraReportFile){
+ try {
+ chimeraReportFile << "queryName\tbestRef\tbestSequenceMismatch\tleftParentChi,rightParentChi\tbreakPointChi\tminMismatchToChimera\tdistToBestMera\tnumParents" << endl;
+ return 0;
+ }catch(exception& e) {
+ m->errorOut(e, "RefChimeraTest", "execute");
+ exit(1);
+ }
+}
//***************************************************************************************************************
-int RefChimeraTest::analyzeQuery(string queryName, string querySeq){
+int RefChimeraTest::analyzeQuery(string queryName, string querySeq, ofstream& chimeraReportFile){
vector<vector<int> > left(numRefSeqs);
vector<int> singleLeft, bestLeft;
int leftParentBi, rightParentBi, breakPointBi;
int minMismatchToChimera = getChimera(left, right, leftParentBi, rightParentBi, breakPointBi, singleLeft, bestLeft, singleRight, bestRight);
-// int minMismatchToTrimera = MAXINT;
-// int leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB;
-
int nMera = 0;
string chimeraRefSeq = "";
if(bestSequenceMismatch - minMismatchToChimera >= 3){// || (minMismatchToChimera == 0 && bestSequenceMismatch != 0)){
-
nMera = 2;
chimeraRefSeq = stitchBimera(leftParentBi, rightParentBi, breakPointBi);
-
-// minMismatchToTrimera = getTrimera(left, right, leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB, singleLeft, bestLeft, singleRight, bestRight);
-//
-// if(minMismatchToChimera - minMismatchToTrimera <= 3){
-// nMera = 2;
-// chimeraRefSeq = stitchBimera(leftParentBi, rightParentBi, breakPointBi);
-// }
-// else{
-// nMera = 3;
-// chimeraRefSeq = stitchTrimera(leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB);
-// }
-
}
else{
nMera = 1;
double distToChimera = calcDistToChimera(querySeq, chimeraRefSeq);
-// double loonIndex = calcLoonIndex(querySeq, referenceSeqs[leftParentBi], referenceSeqs[rightParentBi], breakPointBi, binMatrix);
-
chimeraReportFile << queryName << '\t' << referenceNames[bestMatch] << '\t' << bestSequenceMismatch << '\t';
chimeraReportFile << referenceNames[leftParentBi] << ',' << referenceNames[rightParentBi] << '\t' << breakPointBi << '\t';
chimeraReportFile << minMismatchToChimera << '\t';
-
-// if(nMera == 1){
-// chimeraReportFile << "NA" << '\t' << "NA" << '\t' << "NA";
-// }
-// else{
-// chimeraReportFile << referenceNames[leftParentTri] << ',' << referenceNames[middleParentTri] << ',' << referenceNames[rightParentTri] << '\t' << breakPointTriA << ',' << breakPointTriB << '\t' << minMismatchToTrimera;
-// }
-
- chimeraReportFile << '\t' << distToChimera << '\t' << nMera << endl;
+ chimeraReportFile << '\t' << distToChimera << '\t' << nMera << endl;
return nMera;
}
for(int i=0;i<numRefSeqs;i++){
int lDiffs = 0;
+
for(int l=0;l<alignLength;l++){
-// if(querySeq[l] != '.' && querySeq[l] != referenceSeqs[i][l]){
- if(querySeq[l] != '.' && referenceSeqs[i][l] != '.' && querySeq[l] != referenceSeqs[i][l]){
+ if(querySeq[l] != '.' && referenceSeqs[i][l] != '.' && querySeq[l] != referenceSeqs[i][l] && referenceSeqs[i][l] != 'N'){
lDiffs++;
}
left[i][l] = lDiffs;
int rDiffs = 0;
int index = 0;
for(int l=alignLength-1;l>=0;l--){
-// if(querySeq[l] != '.' && querySeq[l] != referenceSeqs[i][l]){
- if(querySeq[l] != '.' && referenceSeqs[i][l] != '.' && querySeq[l] != referenceSeqs[i][l]){
+ if(querySeq[l] != '.' && referenceSeqs[i][l] != '.' && querySeq[l] != referenceSeqs[i][l] && referenceSeqs[i][l] != 'N'){
rDiffs++;
}
right[i][index++] = rDiffs;
}
+
if(lDiffs < bestSequenceMismatch){
bestSequenceMismatch = lDiffs;
bestRefSeq = i;
for(int i=0;i<alignLength;i++){
// if(querySeq[i] != '.' && chimeraRefSeq[i] != '.'){
- if(chimeraRefSeq[i] != '.'){
- if(querySeq[i] == '-' && chimeraRefSeq[i] == '-'){ /* do nothing */ }
+ if(chimeraRefSeq[i] != '.' && querySeq[i] != '.'){
+ if(querySeq[i] == '-' && chimeraRefSeq[i] == '-' && chimeraRefSeq[i] != 'N'){ /* do nothing */ }
else if(querySeq[i] == chimeraRefSeq[i]){
match++;
}