X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=refchimeratest.cpp;fp=refchimeratest.cpp;h=0000000000000000000000000000000000000000;hb=4a877efa127e56e81a21f53cfdbbfd3bfbe8c4ff;hp=2f397e8f71ec0b8a1d47338b1f12cfc552642ed2;hpb=a6cf29fa4dac0909c7582cb1094151d34093ee76;p=mothur.git

diff --git a/refchimeratest.cpp b/refchimeratest.cpp
deleted file mode 100644
index 2f397e8..0000000
--- a/refchimeratest.cpp
+++ /dev/null
@@ -1,294 +0,0 @@
-/*
- *  refchimeratest.cpp
- *  Mothur
- *
- *  Created by Pat Schloss on 1/31/11.
- *  Copyright 2011 Schloss Lab. All rights reserved.
- *
- */
-
-#include "refchimeratest.h"
-#include "mothur.h"
-
-int MAXINT = numeric_limits<int>::max();
-
-//***************************************************************************************************************
-
-RefChimeraTest::RefChimeraTest(vector<Sequence>& refs){
-
-	m = MothurOut::getInstance();
-
-	numRefSeqs = refs.size();
-
-	referenceSeqs.resize(numRefSeqs);
-	referenceNames.resize(numRefSeqs);
-	for(int i=0;i<numRefSeqs;i++){
-		referenceSeqs[i] = refs[i].getAligned();
-		referenceNames[i] = refs[i].getName();
-	}
-	
-	alignLength = referenceSeqs[0].length();
-
-
-}
-//***************************************************************************************************************
-
-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, ofstream& chimeraReportFile){
-	
-	vector<vector<int> > left(numRefSeqs);
-	vector<int> singleLeft, bestLeft;
-	vector<int> singleRight, bestRight;
-	
-	vector<vector<int> > right(numRefSeqs);
-	for(int i=0;i<numRefSeqs;i++){
-		left[i].assign(alignLength, 0);
-	}
-	right = left;
-	
-	int bestSequenceMismatch = getMismatches(querySeq, left, right, bestMatch);
-	
-	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;
-		chimeraRefSeq = referenceSeqs[bestMatch];
-	}
-	
-	
-	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;
-		
-	return nMera;
-}
-
-/**************************************************************************************************/
-
-int RefChimeraTest::getMismatches(string& querySeq, vector<vector<int> >& left, vector<vector<int> >& right, int& bestRefSeq){
-	
-	int bestSequenceMismatch = MAXINT;
-	
-	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] && 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] && referenceSeqs[i][l] != 'N'){
-				rDiffs++;
-			}			
-			right[i][index++] = rDiffs;
-		}
-		
-		if(lDiffs < bestSequenceMismatch){
-			bestSequenceMismatch = lDiffs;
-			bestRefSeq = i;
-		}
-	}
-	return bestSequenceMismatch;
-}
-
-/**************************************************************************************************/
-
-int RefChimeraTest::getChimera(vector<vector<int> >& left, vector<vector<int> >& right, int& leftParent, int& rightParent, int& breakPoint, vector<int>& singleLeft, vector<int>& bestLeft, vector<int>& singleRight, vector<int>& bestRight){
-	
-	singleLeft.resize(alignLength, MAXINT);
-	bestLeft.resize(alignLength, -1);
-	
-	for(int l=0;l<alignLength;l++){
-		for(int i=0;i<numRefSeqs;i++){
-			if(left[i][l] <= singleLeft[l]){
-				singleLeft[l] = left[i][l];
-				bestLeft[l] = i;
-			}
-		}
-	}
-	
-	singleRight.resize(alignLength, MAXINT);
-	bestRight.resize(alignLength, -1);
-	
-	for(int l=0;l<alignLength;l++){
-		for(int i=0;i<numRefSeqs;i++){
-			if(right[i][l] <= singleRight[l]){
-				singleRight[l] = right[i][l];
-				bestRight[l] = i;
-			}
-		}
-	}
-	
-	int bestChimeraMismatches = MAXINT;
-	leftParent = -1;
-	rightParent = -1;
-	breakPoint = -1;
-	
-	for(int l=0;l<alignLength;l++){
-		int chimera = singleLeft[l] + singleRight[alignLength - l - 1];
-		if(chimera < bestChimeraMismatches){
-			bestChimeraMismatches = chimera;
-			breakPoint = l;
-			leftParent = bestLeft[l];
-			rightParent = bestRight[alignLength - l - 1];
-		}
-	}
-	
-	return bestChimeraMismatches;
-}
-
-/**************************************************************************************************/
-
-int RefChimeraTest::getTrimera(vector<vector<int> >& left, vector<vector<int> >& right, int& leftParent, int& middleParent, int& rightParent, int& breakPointA, int& breakPointB, vector<int>& singleLeft, vector<int>& bestLeft, vector<int>& singleRight, vector<int>& bestRight){
-	
-	int bestTrimeraMismatches = MAXINT;
-	
-	leftParent = -1;
-	middleParent = -1;
-	rightParent = -1;
-	
-	breakPointA = -1;
-	breakPointB = -1;
-	
-	vector<vector<int> > minDelta(alignLength);
-	vector<vector<int> > minDeltaSeq(alignLength);
-	
-	for(int i=0;i<alignLength;i++){
-		minDelta[i].assign(alignLength, MAXINT);
-		minDeltaSeq[i].assign(alignLength, -1);
-	}
-	
-	for(int x=0;x<alignLength;x++){
-		for(int y=x;y<alignLength-1;y++){
-			
-			for(int i=0;i<numRefSeqs;i++){
-				int delta = left[i][y] - left[i][x];
-				
-				if(delta <= minDelta[x][y]){
-					minDelta[x][y] = delta;
-					minDeltaSeq[x][y] = i;					
-				}				
-			}
-			minDelta[x][y] += singleLeft[x] + singleRight[alignLength - y - 2];
-			
-			if(minDelta[x][y] < bestTrimeraMismatches){
-				bestTrimeraMismatches = minDelta[x][y];
-				
-				breakPointA = x;
-				breakPointB = y;
-				
-				leftParent = bestLeft[x];
-				middleParent = minDeltaSeq[x][y];
-				rightParent = bestRight[alignLength - y - 2];				
-			}
-		}		
-	}
-	return bestTrimeraMismatches;
-}
-
-/**************************************************************************************************/
-
-string RefChimeraTest::stitchBimera(int leftParent, int rightParent, int breakPoint){
-	
-	string chimeraRefSeq = referenceSeqs[leftParent].substr(0, breakPoint) + referenceSeqs[rightParent].substr(breakPoint);
-	return chimeraRefSeq;
-	
-}
-
-/**************************************************************************************************/
-
-string RefChimeraTest::stitchTrimera(int leftParent, int middleParent, int rightParent, int breakPointA, int breakPointB){
-	
-	string chimeraRefSeq = referenceSeqs[leftParent].substr(0, breakPointA) + referenceSeqs[middleParent].substr(breakPointA, breakPointB-breakPointA) + referenceSeqs[rightParent].substr(breakPointB);
-	
-	return chimeraRefSeq;
-}
-
-/**************************************************************************************************/
-
-double RefChimeraTest::calcDistToChimera(string& querySeq, string& chimeraRefSeq){
-	
-	int match = 0;
-	int mismatch = 0;
-	
-	for(int i=0;i<alignLength;i++){
-//		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++;
-			}
-			else{
-				mismatch++;
-			}			
-		}
-	}
-	
-	return (double)mismatch / (double)(mismatch + match);	
-}
-
-//***************************************************************************************************************
-
-int RefChimeraTest::getClosestRefIndex(){
-
-	return bestMatch;
-	
-}
-
-//***************************************************************************************************************