X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=refchimeratest.cpp;h=2f397e8f71ec0b8a1d47338b1f12cfc552642ed2;hb=006601d68abe8d0061f77e8d28323b160750e343;hp=246d4394a58ae7f5d65bf78c8b0d19a426ea30c6;hpb=5e5253ff472de3c6349e562d2580873287be0c65;p=mothur.git

diff --git a/refchimeratest.cpp b/refchimeratest.cpp
index 246d439..2f397e8 100644
--- a/refchimeratest.cpp
+++ b/refchimeratest.cpp
@@ -14,11 +14,10 @@ int MAXINT = numeric_limits<int>::max();
 
 //***************************************************************************************************************
 
-RefChimeraTest::RefChimeraTest(vector<Sequence>& refs, string chimeraReportFileName){
+RefChimeraTest::RefChimeraTest(vector<Sequence>& refs){
 
 	m = MothurOut::getInstance();
 
-	m->openOutputFile(chimeraReportFileName, chimeraReportFile);
 	numRefSeqs = refs.size();
 
 	referenceSeqs.resize(numRefSeqs);
@@ -30,14 +29,22 @@ RefChimeraTest::RefChimeraTest(vector<Sequence>& refs, string chimeraReportFileN
 	
 	alignLength = referenceSeqs[0].length();
 
-	chimeraReportFile << "queryName\tbestRef\tbestSequenceMismatch\tleftParentChi,rightParentChi\tbreakPointChi\tminMismatchToChimera\t";
-	chimeraReportFile << "leftParentChi,middleParentTri,rightParentChi\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;
@@ -54,30 +61,35 @@ int RefChimeraTest::analyzeQuery(string queryName, string querySeq){
 	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 minMismatchToTrimera = MAXINT;
+//	int leftParentTri, middleParentTri, rightParentTri, breakPointTriA, breakPointTriB;
 	
 	int nMera = 0;
 	string chimeraRefSeq = "";
 	
-	if(bestSequenceMismatch - minMismatchToChimera <= 3){
-		nMera = 1;
-		chimeraRefSeq = referenceSeqs[bestMatch];
-	}
-	else {
+	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);
-		}
+//		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;
+		chimeraRefSeq = referenceSeqs[bestMatch];
+	}
+	
+	
 	double distToChimera = calcDistToChimera(querySeq, chimeraRefSeq);
 	
 //	double loonIndex = calcLoonIndex(querySeq, referenceSeqs[leftParentBi], referenceSeqs[rightParentBi], breakPointBi, binMatrix);		
@@ -86,12 +98,12 @@ int RefChimeraTest::analyzeQuery(string queryName, string querySeq){
 	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;	
-	}
+//	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;
 		
@@ -107,8 +119,11 @@ int RefChimeraTest::getMismatches(string& querySeq, vector<vector<int> >& left,
 	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] != '.' && 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;
@@ -117,11 +132,13 @@ int RefChimeraTest::getMismatches(string& querySeq, vector<vector<int> >& left,
 		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] != '.' && 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;
@@ -251,8 +268,9 @@ double RefChimeraTest::calcDistToChimera(string& querySeq, string& chimeraRefSeq
 	int mismatch = 0;
 	
 	for(int i=0;i<alignLength;i++){
-		if(querySeq[i] != '.' && chimeraRefSeq[i] != '.'){
-			if(querySeq[i] == '-' && chimeraRefSeq[i] == '-'){	/*	do nothing	*/	}
+//		if(querySeq[i] != '.' && chimeraRefSeq[i] != '.'){
+		if(chimeraRefSeq[i] != '.' && querySeq[i] != '.'){
+			if(querySeq[i] == '-' && chimeraRefSeq[i] == '-' && chimeraRefSeq[i] != 'N'){	/*	do nothing	*/	}
 			else if(querySeq[i] == chimeraRefSeq[i]){
 				match++;
 			}