C++ Sequence::getAligned方法代码示例

void DistanceDB::addSequence(Sequence seq) {
	try {
		//are the template sequences aligned
		if (!isAligned(seq.getAligned())) {
			templateAligned = false;
			m->mothurOut(seq.getName() + " is not aligned. Sequences must be aligned to use the distance method.");
			m->mothurOutEndLine(); 
		}
		
		if (templateSeqsLength == 0) { templateSeqsLength = seq.getAligned().length(); }
				
		data.push_back(seq);
	}
	catch(exception& e) {
		m->errorOut(e, "DistanceDB", "addSequence");
		exit(1);
	}	
}

int ChimeraPintailCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& outMPI, MPI_File& outAccMPI, vector<unsigned long long>& MPIPos){
	try {
				
		MPI_Status status; 
		int pid;
		MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
		
		for(int i=0;i<num;i++){
			
			if (m->control_pressed) {	return 1;	}
			
			//read next sequence
			int length = MPIPos[start+i+1] - MPIPos[start+i];
	
			char* buf4 = new char[length];
			MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
			
			string tempBuf = buf4;
			if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length);  }
			istringstream iss (tempBuf,istringstream::in);
			delete buf4;

			Sequence* candidateSeq = new Sequence(iss);  m->gobble(iss);
				
			if (candidateSeq->getName() != "") { //incase there is a commented sequence at the end of a file
				
				if	(candidateSeq->getAligned().length() != templateSeqsLength) {  //chimeracheck does not require seqs to be aligned
					m->mothurOut(candidateSeq->getName() + " is not the same length as the template sequences. Skipping."); m->mothurOutEndLine();
				}else{
					//find chimeras
					chimera->getChimeras(candidateSeq);
					
					if (m->control_pressed) {	delete candidateSeq; return 1;	}
		
					//print results
					chimera->print(outMPI, outAccMPI);
				}
			}
			delete candidateSeq;
			
			//report progress
			if((i+1) % 100 == 0){  cout << "Processing sequence: " << (i+1) << endl;			}
		}
		//report progress
		if(num % 100 != 0){		cout << "Processing sequence: " << num << endl;	 	}
		
				
		return 0;
	}
	catch(exception& e) {
		m->errorOut(e, "ChimeraPintailCommand", "driverMPI");
		exit(1);
	}
}

map<int, int> DeCalculator::trimSeqs(Sequence& query, vector<Sequence>& topMatches) {
	try {
		
		int frontPos = 0;  //should contain first position in all seqs that is not a gap character
		int rearPos = query.getAligned().length();
		
		//********find first position in topMatches that is a non gap character***********//
		//find first position all query seqs that is a non gap character
		for (int i = 0; i < topMatches.size(); i++) {
			
			string aligned = topMatches[i].getAligned();
			int pos = 0;
			
			//find first spot in this seq
			for (int j = 0; j < aligned.length(); j++) {
				if (isalpha(aligned[j])) {
					pos = j;
					break;
				}
			}
			
			//save this spot if it is the farthest
			if (pos > frontPos) { frontPos = pos; }
		}
		
		
		string aligned = query.getAligned();
		int pos = 0;
			
		//find first position in query that is a non gap character
		for (int j = 0; j < aligned.length(); j++) {
			if (isalpha(aligned[j])) {
				pos = j;
				break;
			}
		}
		
		//save this spot if it is the farthest
		if (pos > frontPos) { frontPos = pos; }
		
		
		//********find last position in topMatches that is a non gap character***********//
		for (int i = 0; i < topMatches.size(); i++) {
			
			string aligned = topMatches[i].getAligned();
			int pos = aligned.length();
			
			//find first spot in this seq
			for (int j = aligned.length()-1; j >= 0; j--) {
				if (isalpha(aligned[j])) {
					pos = j;
					break;
				}
			}
			
			//save this spot if it is the farthest
			if (pos < rearPos) { rearPos = pos; }
		}
		
		
		aligned = query.getAligned();
		pos = aligned.length();
		
		//find last position in query that is a non gap character
		for (int j = aligned.length()-1; j >= 0; j--) {
			if (isalpha(aligned[j])) {
				pos = j;
				break;
			}
		}
		
		//save this spot if it is the farthest
		if (pos < rearPos) { rearPos = pos; }
		
		map<int, int> trimmedPos;
		//check to make sure that is not whole seq
		if ((rearPos - frontPos - 1) <= 0) {  
			query.setAligned("");
			//trim topMatches
			for (int i = 0; i < topMatches.size(); i++) {
				topMatches[i].setAligned("");
			}
			
		}else {

			//trim query
			string newAligned = query.getAligned();
			newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
			query.setAligned(newAligned);
			
			//trim topMatches
			for (int i = 0; i < topMatches.size(); i++) {
				newAligned = topMatches[i].getAligned();
				newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1));
				topMatches[i].setAligned(newAligned);
			}
			
			for (int i = 0; i < newAligned.length(); i++) {
				trimmedPos[i] = i+frontPos;
			}
//.........这里部分代码省略.........

//***************************************************************************************************************
//gets closest matches to each end, since chimeras will most likely have different parents on each end
vector<Sequence> DeCalculator::findClosest(Sequence querySeq, vector<Sequence*>& thisTemplate, vector<Sequence*>& thisFilteredTemplate, int numWanted, int minSim) {
	try {
		//indexes.clear();
		
		vector<Sequence> seqsMatches;  
		
		vector<SeqDist> distsLeft;
		vector<SeqDist> distsRight;
		
		Dist* distcalculator = new eachGapDist();
		
		string queryUnAligned = querySeq.getUnaligned();
		int numBases = int(queryUnAligned.length() * 0.33);
		
		string leftQuery = ""; //first 1/3 of the sequence
		string rightQuery = ""; //last 1/3 of the sequence
		string queryAligned = querySeq.getAligned();
		
		//left side
		bool foundFirstBase = false;
		int baseCount = 0;
		int leftSpot = 0;
		int firstBaseSpot = 0;
		for (int i = 0; i < queryAligned.length(); i++) {
			//if you are a base
			if (isalpha(queryAligned[i])) {		
				baseCount++; 
				if (!foundFirstBase) {   foundFirstBase = true;  firstBaseSpot = i;  }
			}
			
			//eliminate opening .'s
			if (foundFirstBase) {   leftQuery += queryAligned[i];  }
			//if you have 1/3
			if (baseCount >= numBases) {  leftSpot = i; break; } //first 1/3
		}
		
		//right side - count through another 1/3, so you are at last third
		baseCount = 0;
		int rightSpot = 0;
		for (int i = leftSpot; i < queryAligned.length(); i++) {
			//if you are a base
			if (isalpha(queryAligned[i])) {		baseCount++;	}
			//if you have 1/3
			if (baseCount > numBases + 1) { rightSpot = i;  break; } //last 1/3
		}
		
		//trim end
		//find last position in query that is a non gap character
		int lastBaseSpot = queryAligned.length()-1;
		for (int j = queryAligned.length()-1; j >= 0; j--) {
			if (isalpha(queryAligned[j])) {
				lastBaseSpot = j;
				break;
			}
		}
		rightQuery = queryAligned.substr(rightSpot, (lastBaseSpot-rightSpot+1)); //sequence from pos spot to end
		
		Sequence queryLeft(querySeq.getName(), leftQuery);
		Sequence queryRight(querySeq.getName(), rightQuery);
		
//cout << querySeq->getName() << '\t' << leftSpot << '\t' << rightSpot << '\t' << firstBaseSpot << '\t' << lastBaseSpot << endl;
//cout << queryUnAligned.length() << '\t' << queryLeft.getUnaligned().length() << '\t' << queryRight.getUnaligned().length() << endl;
		for(int j = 0; j < thisFilteredTemplate.size(); j++){
			
			string dbAligned = thisFilteredTemplate[j]->getAligned();
			string leftDB = dbAligned.substr(firstBaseSpot, (leftSpot-firstBaseSpot+1)); //first 1/3 of the sequence
			string rightDB = dbAligned.substr(rightSpot, (lastBaseSpot-rightSpot+1)); //last 1/3 of the sequence

			Sequence dbLeft(thisFilteredTemplate[j]->getName(), leftDB);
			Sequence dbRight(thisFilteredTemplate[j]->getName(), rightDB);

			distcalculator->calcDist(queryLeft, dbLeft);
			float distLeft = distcalculator->getDist();
			
			distcalculator->calcDist(queryRight, dbRight);
			float distRight = distcalculator->getDist();

			SeqDist subjectLeft;
			subjectLeft.seq = NULL;
			subjectLeft.dist = distLeft;
			subjectLeft.index = j;
			
			distsLeft.push_back(subjectLeft);
			
			SeqDist subjectRight;
			subjectRight.seq = NULL;
			subjectRight.dist = distRight;
			subjectRight.index = j;
			
			distsRight.push_back(subjectRight);

		}
		
		delete distcalculator;
		
		//sort by smallest distance
		sort(distsRight.begin(), distsRight.end(), compareSeqDist);
		sort(distsLeft.begin(), distsLeft.end(), compareSeqDist);
//.........这里部分代码省略.........

//***************************************************************************************************************
vector<Sequence*> MothurChimera::readSeqs(string file) {
	try {
		
		vector<Sequence*> container;
		int count = 0;
		length = 0;
		unaligned = false;
		ReferenceDB* rdb = ReferenceDB::getInstance();
		
		if (file == "saved") {
			
			
			m->mothurOutEndLine();  m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory.");	m->mothurOutEndLine();
			
			for (int i = 0; i < rdb->referenceSeqs.size(); i++) {
				Sequence* temp = new Sequence(rdb->referenceSeqs[i].getName(), rdb->referenceSeqs[i].getAligned());
				
				if (count == 0) {  length = temp->getAligned().length();  count++;  } //gets first seqs length
				else if (length != temp->getAligned().length()) {	unaligned = true;	}
				
				if (temp->getName() != "") {  container.push_back(temp);  }
			}
			
			templateFileName = rdb->getSavedReference();
			
		}else {
			
			m->mothurOut("Reading sequences from " + file + "..."); cout.flush();
			
			
			ifstream in;
			m->openInputFile(file, in);
			
			//read in seqs and store in vector
			while(!in.eof()){
				
				if (m->control_pressed) { return container; }
				
				Sequence* current = new Sequence(in);  m->gobble(in);
				
				if (count == 0) {  length = current->getAligned().length();  count++;  } //gets first seqs length
				else if (length != current->getAligned().length()) {   unaligned = true;	}
							
				if (current->getName() != "") {  
					container.push_back(current);  
					if (rdb->save) { rdb->referenceSeqs.push_back(*current); }
				}
			}
			in.close();
		
			m->mothurOut("Done."); m->mothurOutEndLine();
			
			filterString = (string(container[0]->getAligned().length(), '1'));
		}
		
		return container;
	}
	catch(exception& e) {
		m->errorOut(e, "MothurChimera", "readSeqs");
		exit(1);
	}
}

//**********************************************************************************************************************
int AlignCommand::driver(linePair* filePos, string alignFName, string reportFName, string accnosFName, string filename){
	try {
		ofstream alignmentFile;
		m->openOutputFile(alignFName, alignmentFile);
		
		ofstream accnosFile;
		m->openOutputFile(accnosFName, accnosFile);
		
		NastReport report(reportFName);
		
		ifstream inFASTA;
		m->openInputFile(filename, inFASTA);

		inFASTA.seekg(filePos->start);

		bool done = false;
		int count = 0;
		
		//moved this into driver to avoid deep copies in windows paralellized version
		Alignment* alignment;
		int longestBase = templateDB->getLongestBase();
        if (m->debug) { m->mothurOut("[DEBUG]: template longest base = "  + toString(templateDB->getLongestBase()) + " \n"); }
		if(align == "gotoh")			{	alignment = new GotohOverlap(gapOpen, gapExtend, match, misMatch, longestBase);			}
		else if(align == "needleman")	{	alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);				}
		else if(align == "blast")		{	alignment = new BlastAlignment(gapOpen, gapExtend, match, misMatch);		}
		else if(align == "noalign")		{	alignment = new NoAlign();													}
		else {
			m->mothurOut(align + " is not a valid alignment option. I will run the command using needleman.");
			m->mothurOutEndLine();
			alignment = new NeedlemanOverlap(gapOpen, match, misMatch, longestBase);
		}
	
		while (!done) {
			
			if (m->control_pressed) {  break; }
			
			Sequence* candidateSeq = new Sequence(inFASTA);  m->gobble(inFASTA);
			report.setCandidate(candidateSeq);

			int origNumBases = candidateSeq->getNumBases();
			string originalUnaligned = candidateSeq->getUnaligned();
			int numBasesNeeded = origNumBases * threshold;
	
			if (candidateSeq->getName() != "") { //incase there is a commented sequence at the end of a file
				if (candidateSeq->getUnaligned().length()+1 > alignment->getnRows()) {
                    if (m->debug) { m->mothurOut("[DEBUG]: " + candidateSeq->getName() + " " + toString(candidateSeq->getUnaligned().length()) + " " + toString(alignment->getnRows()) + " \n"); }
					alignment->resize(candidateSeq->getUnaligned().length()+2);
				}
				Sequence temp = templateDB->findClosestSequence(candidateSeq);
				Sequence* templateSeq = new Sequence(temp.getName(), temp.getAligned());
				
				float searchScore = templateDB->getSearchScore();
								
				Nast* nast = new Nast(alignment, candidateSeq, templateSeq);
		
				Sequence* copy;
				
				Nast* nast2;
				bool needToDeleteCopy = false;  //this is needed in case you have you enter the ifs below
												//since nast does not make a copy of hte sequence passed, and it is used by the reporter below
												//you can't delete the copy sequence til after you report, but you may choose not to create it in the first place
												//so this bool tells you if you need to delete it
												
				//if there is a possibility that this sequence should be reversed
				if (candidateSeq->getNumBases() < numBasesNeeded) {
					
					string wasBetter =  "";
					//if the user wants you to try the reverse
					if (flip) {
				
						//get reverse compliment
						copy = new Sequence(candidateSeq->getName(), originalUnaligned);
						copy->reverseComplement();
                        
                        if (m->debug) { m->mothurOut("[DEBUG]: flipping "  + candidateSeq->getName() + " \n"); }
						
						//rerun alignment
						Sequence temp2 = templateDB->findClosestSequence(copy);
						Sequence* templateSeq2 = new Sequence(temp2.getName(), temp2.getAligned());
                        
                        if (m->debug) { m->mothurOut("[DEBUG]: closest template "  + temp2.getName() + " \n"); }
						
						searchScore = templateDB->getSearchScore();
						
						nast2 = new Nast(alignment, copy, templateSeq2);
                        
                        if (m->debug) { m->mothurOut("[DEBUG]: completed Nast2 "  + candidateSeq->getName() + " flipped numBases = " + toString(copy->getNumBases()) + " old numbases = " + toString(candidateSeq->getNumBases()) +" \n"); }
			
						//check if any better
						if (copy->getNumBases() > candidateSeq->getNumBases()) {
							candidateSeq->setAligned(copy->getAligned());  //use reverse compliments alignment since its better
                            delete templateSeq;
							templateSeq = templateSeq2;
							delete nast;
							nast = nast2;
							needToDeleteCopy = true;
							wasBetter = "\treverse complement produced a better alignment, so mothur used the reverse complement.";
						}else{  
							wasBetter = "\treverse complement did NOT produce a better alignment so it was not used, please check sequence.";
//.........这里部分代码省略.........

int ChimeraPintailCommand::driver(linePair* filePos, string outputFName, string filename, string accnos){
	try {
		ofstream out;
		m->openOutputFile(outputFName, out);
		
		ofstream out2;
		m->openOutputFile(accnos, out2);
		
		ifstream inFASTA;
		m->openInputFile(filename, inFASTA);

		inFASTA.seekg(filePos->start);

		bool done = false;
		int count = 0;
	
		while (!done) {
				
			if (m->control_pressed) {	return 1;	}
		
			Sequence* candidateSeq = new Sequence(inFASTA);  m->gobble(inFASTA);
				
			if (candidateSeq->getName() != "") { //incase there is a commented sequence at the end of a file
				
				if (candidateSeq->getAligned().length() != templateSeqsLength)  {  //chimeracheck does not require seqs to be aligned
					m->mothurOut(candidateSeq->getName() + " is not the same length as the template sequences. Skipping."); m->mothurOutEndLine();
				}else{
					//find chimeras
					chimera->getChimeras(candidateSeq);
					
					if (m->control_pressed) {	delete candidateSeq; return 1;	}
		
					//print results
					chimera->print(out, out2);
				}
				count++;
			}
			delete candidateSeq;
			
			#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
				unsigned long long pos = inFASTA.tellg();
				if ((pos == -1) || (pos >= filePos->end)) { break; }
			#else
				if (inFASTA.eof()) { break; }
			#endif
			
			//report progress
			if((count) % 100 == 0){	m->mothurOutJustToScreen("Processing sequence: " + toString(count) + "\n"); 		}
		}
		//report progress
		if((count) % 100 != 0){	m->mothurOutJustToScreen("Processing sequence: " + toString(count) + "\n"); 		}
		
		out.close();
		out2.close();
		inFASTA.close();
				
		return count;
	}
	catch(exception& e) {
		m->errorOut(e, "ChimeraPintailCommand", "driver");
		exit(1);
	}
}

//.........这里部分代码省略.........
				Sequence* templateSeq = &temp;
				
				float searchScore = templateDB->getSearchScore();
								
				Nast* nast = new Nast(alignment, candidateSeq, templateSeq);
				Sequence* copy;
				
				Nast* nast2;
				bool needToDeleteCopy = false;  //this is needed in case you have you enter the ifs below
												//since nast does not make a copy of hte sequence passed, and it is used by the reporter below
												//you can't delete the copy sequence til after you report, but you may choose not to create it in the first place
												//so this bool tells you if you need to delete it
												
				//if there is a possibility that this sequence should be reversed
				if (candidateSeq->getNumBases() < numBasesNeeded) {
					
					string wasBetter = "";
					//if the user wants you to try the reverse
					if (flip) {
						//get reverse compliment
						copy = new Sequence(candidateSeq->getName(), originalUnaligned);
						copy->reverseComplement();
						
						//rerun alignment
						Sequence temp2 = templateDB->findClosestSequence(copy);
						Sequence* templateSeq2 = &temp2;
						
						searchScore = templateDB->getSearchScore();
						
						nast2 = new Nast(alignment, copy, templateSeq2);
			
						//check if any better
						if (copy->getNumBases() > candidateSeq->getNumBases()) {
							candidateSeq->setAligned(copy->getAligned());  //use reverse compliments alignment since its better
							templateSeq = templateSeq2; 
							delete nast;
							nast = nast2;
							needToDeleteCopy = true;
							wasBetter = "\treverse complement produced a better alignment, so mothur used the reverse complement.";
						}else{  
							wasBetter = "\treverse complement did NOT produce a better alignment, please check sequence.";
							delete nast2;
							delete copy;	
						}
					}
					
					//create accnos file with names
					outputString = candidateSeq->getName() + wasBetter + "\n";
					
					//send results to parent
					int length = outputString.length();

					char* buf = new char[length];
					memcpy(buf, outputString.c_str(), length);
				
					MPI_File_write_shared(accnosFile, buf, length, MPI_CHAR, &statusAccnos);
					delete buf;
					MPIWroteAccnos = true;
				}
				
				report.setTemplate(templateSeq);
				report.setSearchParameters(search, searchScore);
				report.setAlignmentParameters(align, alignment);
				report.setNastParameters(*nast);
	
				outputString =  ">" + candidateSeq->getName() + "\n" + candidateSeq->getAligned() + "\n";

//**********************************************************************************************************************
void alignDriver(alignStruct* params) {
	try {
        NastReport report;
		
		ifstream inFASTA;
		params->util.openInputFile(params->inputFilename, inFASTA);

		inFASTA.seekg(params->filePos.start);

		bool done = false;
        
		long long count = 0;
        long long numFlipped_0 = 0;
        long long numFlipped_1 = 0;
		
		//moved this into driver to avoid deep copies in windows paralellized version
		Alignment* alignment;
		int longestBase = params->templateDB->getLongestBase();
        if (params->m->getDebug()) { params->m->mothurOut("[DEBUG]: template longest base = "  + toString(longestBase) + " \n"); }
		if(params->alignMethod == "gotoh")			{	alignment = new GotohOverlap(params->gapOpen, params->gapExtend, params->match, params->misMatch, longestBase);			}
		else if(params->alignMethod == "needleman")	{	alignment = new NeedlemanOverlap(params->gapOpen, params->match, params->misMatch, longestBase);				}
		else if(params->alignMethod == "blast")		{	alignment = new BlastAlignment(params->gapOpen, params->gapExtend, params->match, params->misMatch);		}
		else if(params->alignMethod == "noalign")		{	alignment = new NoAlign();													}
		else {
			params->m->mothurOut(params->alignMethod + " is not a valid alignment option. I will run the command using needleman.");
			params->m->mothurOutEndLine();
			alignment = new NeedlemanOverlap(params->gapOpen, params->match, params->misMatch, longestBase);
		}
	
		while (!done) {
			
			if (params->m->getControl_pressed()) {  break; }
			
			Sequence* candidateSeq = new Sequence(inFASTA);  params->util.gobble(inFASTA);
			report.setCandidate(candidateSeq);

			int origNumBases = candidateSeq->getNumBases();
			string originalUnaligned = candidateSeq->getUnaligned();
			int numBasesNeeded = origNumBases * params->threshold;
	
			if (candidateSeq->getName() != "") { //incase there is a commented sequence at the end of a file
				if (candidateSeq->getUnaligned().length()+1 > alignment->getnRows()) {
                    if (params->m->getDebug()) { params->m->mothurOut("[DEBUG]: " + candidateSeq->getName() + " " + toString(candidateSeq->getUnaligned().length()) + " " + toString(alignment->getnRows()) + " \n"); }
					alignment->resize(candidateSeq->getUnaligned().length()+2);
				}
                
                float searchScore;
				Sequence temp = params->templateDB->findClosestSequence(candidateSeq, searchScore);
				Sequence* templateSeq = new Sequence(temp.getName(), temp.getAligned());
								
				Nast* nast = new Nast(alignment, candidateSeq, templateSeq);
		
				Sequence* copy;
				
				Nast* nast2;
				bool needToDeleteCopy = false;  //this is needed in case you have you enter the ifs below
												//since nast does not make a copy of hte sequence passed, and it is used by the reporter below
												//you can't delete the copy sequence til after you report, but you may choose not to create it in the first place
												//so this bool tells you if you need to delete it
												
				//if there is a possibility that this sequence should be reversed
				if (candidateSeq->getNumBases() < numBasesNeeded) {
					numFlipped_1++;
					string wasBetter =  "";
					//if the user wants you to try the reverse
					if (params->flip) {
				
						//get reverse compliment
						copy = new Sequence(candidateSeq->getName(), originalUnaligned);
						copy->reverseComplement();
                        
                        if (params->m->getDebug()) { params->m->mothurOut("[DEBUG]: flipping "  + candidateSeq->getName() + " \n"); }
						
						//rerun alignment
						Sequence temp2 = params->templateDB->findClosestSequence(copy, searchScore);
						Sequence* templateSeq2 = new Sequence(temp2.getName(), temp2.getAligned());
                        
                        if (params->m->getDebug()) { params->m->mothurOut("[DEBUG]: closest template "  + temp2.getName() + " \n"); }
						
						nast2 = new Nast(alignment, copy, templateSeq2);
                        
                        if (params->m->getDebug()) { params->m->mothurOut("[DEBUG]: completed Nast2 "  + candidateSeq->getName() + " flipped numBases = " + toString(copy->getNumBases()) + " old numbases = " + toString(candidateSeq->getNumBases()) +" \n"); }
			
						//check if any better
						if (copy->getNumBases() > candidateSeq->getNumBases()) {
							candidateSeq->setAligned(copy->getAligned());  //use reverse compliments alignment since its better
                            delete templateSeq;
							templateSeq = templateSeq2;
							delete nast;
							nast = nast2;
							needToDeleteCopy = true;
							wasBetter = "\treverse complement produced a better alignment, so mothur used the reverse complement.";
                            numFlipped_0++;
						}else{  
							wasBetter = "\treverse complement did NOT produce a better alignment so it was not used, please check sequence.";
							delete nast2;
                            delete templateSeq2;
							delete copy;	
						}
//.........这里部分代码省略.........

//**********************************************************************************************************************
string ChopSeqsCommand::getChopped(Sequence seq, string& qualValues) {
	try {
		string temp = seq.getAligned();
		string tempUnaligned = seq.getUnaligned();
		
		if (countGaps) {
			//if needed trim sequence
			if (keep == "front") {//you want to keep the beginning
				int tempLength = temp.length();

				if (tempLength > numbases) { //you have enough bases to remove some
				
					int stopSpot = 0;
					int numBasesCounted = 0;
					
					for (int i = 0; i < temp.length(); i++) {
						//eliminate N's
                        if (!keepN) { if (toupper(temp[i]) == 'N') { temp[i] = '.'; } }
						
						numBasesCounted++; 
						
						if (numBasesCounted >= numbases) { stopSpot = i; break; }
					}
					
					if (stopSpot == 0) { temp = ""; }
					else {  temp = temp.substr(0, stopSpot+1);  }
							
				}else { 
					if (!Short) { temp = ""; } //sequence too short
				}
			}else { //you are keeping the back
				int tempLength = temp.length();
				if (tempLength > numbases) { //you have enough bases to remove some
					
					int stopSpot = 0;
					int numBasesCounted = 0;
					
					for (int i = (temp.length()-1); i >= 0; i--) {
						//eliminate N's
                        if (!keepN) { if (toupper(temp[i]) == 'N') { temp[i] = '.'; } }
						
						numBasesCounted++; 

						if (numBasesCounted >= numbases) { stopSpot = i; break; }
					}
				
					if (stopSpot == 0) { temp = ""; }
					else {  temp = temp.substr(stopSpot+1);  }
				}else { 
					if (!Short) { temp = ""; } //sequence too short
				}
			}

		}else{
				
			//if needed trim sequence
			if (keep == "front") {//you want to keep the beginning
				int tempLength = tempUnaligned.length();

				if (tempLength > numbases) { //you have enough bases to remove some
					
					int stopSpot = 0;
					int numBasesCounted = 0;
					
					for (int i = 0; i < temp.length(); i++) {
						//eliminate N's
                        if (!keepN) {
                            if (toupper(temp[i]) == 'N') {
                                temp[i] = '.';
                                tempLength--;
                                if (tempLength < numbases) { stopSpot = 0; break; }
                            }
                        }
						if(isalpha(temp[i])) { numBasesCounted++; }
						
						if (numBasesCounted >= numbases) { stopSpot = i; break; }
					}
					
					if (stopSpot == 0) { temp = ""; }
					else {  temp = temp.substr(0, stopSpot+1);  }
                    
					qualValues = seq.getName() +'\t' + toString(0) + '\t' + toString(stopSpot+1) + '\n';
                    
				}else { 
					if (!Short) { temp = ""; qualValues = seq.getName() +'\t' + toString(0) + '\t' + toString(0) + '\n'; } //sequence too short
                    else { qualValues = seq.getName() +'\t' + toString(0) + '\t' + toString(tempLength) + '\n'; }
				}				
			}else { //you are keeping the back
				int tempLength = tempUnaligned.length();
				if (tempLength > numbases) { //you have enough bases to remove some
					
					int stopSpot = 0;
					int numBasesCounted = 0;
					
					for (int i = (temp.length()-1); i >= 0; i--) {
                        if (!keepN) {
                            //eliminate N's
                            if (toupper(temp[i]) == 'N') {
                                temp[i] = '.';
//.........这里部分代码省略.........