+++ /dev/null
-/*
- * 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;
-
-}
-
-//***************************************************************************************************************
projects / mothur.git / blobdiff