C++ Alignment::align方法代码示例

bool CellTrack::tracking()
{	
	assert(m_numFrames >= 2);
	/*
	 * 0. reset containers
	 */
	cout<<endl<<"start tracking ..."<<endl;
	m_frames.clear();
	m_frames.resize(m_numFrames);
	if(! m_tracks.empty())
	{
		vector<Track>::iterator it = m_tracks.begin();
		while(it != m_tracks.end())
		{
			Cell* p = (*it).entry;
			assert(p!=NULL);
			Cell* next = p;
			while(p->next)
			{
				next = p->next;
				delete p;
				p = next;
			}
			delete p;
			it++;
		}
		m_tracks.clear();
		m_numTracks = 0;
	}
	/*
	 * 1. local variables and member vaiable getParams
	 */
	vector<Frame> t_frames;
	Alignment alg;
	t_frames.resize(2*m_numFrames - 2);
	double * row = NULL; // the delete operation of row is done by alg

	/*
	 * 2. the first round of comparison
	 */
	
	// 2.1 load the first tree

	int i = 0; // the entry of t_frames
	int t = 1; // the entry of time
	// 2.2 go to the loop of m_numFrames - 1 times comparison
	//     and m_numFrames -1 matches will bed set
	
	m_matches.resize(m_numFrames - 1);
	cout<<"\tThe first round alignment "<<endl;
	for(int t = 1; t < m_numFrames; t++)
	{
		cout<<"\talign tree"<<t<<" and tree"<<(t+1)<<endl;
		
		alg.align(m_trees[t-1], m_trees[t]);
		
		row = alg.result();

		create_match(row, m_trees[t-1].nodeCount(), m_trees[t].nodeCount(), m_matches[t-1]);

		create_frames(m_matches[t-1],t_frames[i],t_frames[i+1]);
				
		alg.clear();
		
		i = i + 2;
	}
	
	/*
	 * 3. the second round of alignment to merge nodes from tim 1 ~ m_numFrames -2
	 */
	cout<<endl<<"\tThe second round alignment "<<endl;
	i = 0;
	t = 0;
	m_frames[t++] = t_frames[i++];
	
	while( t <= m_numFrames - 2)
	{
		cout<<"\talign tree"<<(t+1)<<endl;
		
		alg.align(m_trees[t],t_frames[i], t_frames[i+1]);
		
		row = alg.result();
		
		m_frames[t++] = merge_frames(m_trees[t], row , t_frames[i], t_frames[i+1]);
		
		alg.clear();
				
		i = i+2;
	}
	// now t = m_numFrames - 1,  i+1 = 2*m_numFrames -2
	m_frames[t] = t_frames[i];
	
	t_frames.clear();
	
	/*
	 * 4. set tracks and trackId for each cell
	 */
	t = 0;
	int trackId = 1;         //the cell id starts from 1, the cell id = 0 is the background
	while(t < m_numFrames)
//.........这里部分代码省略.........

/********************************************************************/
TrimOligos::~TrimOligos() {}
//*******************************************************************/
int TrimOligos::stripBarcode(Sequence& seq, QualityScores& qual, int& group){
	try {
		
		string rawSequence = seq.getUnaligned();
		int success = bdiffs + 1;	//guilty until proven innocent
		
		//can you find the barcode
		for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
			string oligo = it->first;
			if(rawSequence.length() < oligo.length()){	//let's just assume that the barcodes are the same length
				success = bdiffs + 10;					//if the sequence is shorter than the barcode then bail out
				break;	
			}
			
			if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
				group = it->second;
				seq.setUnaligned(rawSequence.substr(oligo.length()));
				
				if(qual.getName() != ""){
					qual.trimQScores(oligo.length(), -1);
				}
				
				success = 0;
				break;
			}
		}
		
		//if you found the barcode or if you don't want to allow for diffs
		if ((bdiffs == 0) || (success == 0)) { return success;  }
		
		else { //try aligning and see if you can find it
			
			int maxLength = 0;
			
			Alignment* alignment;
			if (barcodes.size() > 0) {
				map<string,int>::iterator it=barcodes.begin();
				
				for(it;it!=barcodes.end();it++){
					if(it->first.length() > maxLength){
						maxLength = it->first.length();
					}
				}
				alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1));  
				
			}else{ alignment = NULL; } 
			
			//can you find the barcode
			int minDiff = 1e6;
			int minCount = 1;
			int minGroup = -1;
			int minPos = 0;
			
			for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
				string oligo = it->first;
				//				int length = oligo.length();
				
				if(rawSequence.length() < maxLength){	//let's just assume that the barcodes are the same length
					success = bdiffs + 10;
					break;
				}
				
				//use needleman to align first barcode.length()+numdiffs of sequence to each barcode
				alignment->align(oligo, rawSequence.substr(0,oligo.length()+bdiffs));
				oligo = alignment->getSeqAAln();
				string temp = alignment->getSeqBAln();
				
				int alnLength = oligo.length();
				
				for(int i=oligo.length()-1;i>=0;i--){
					if(oligo[i] != '-'){	alnLength = i+1;	break;	}
				}
				oligo = oligo.substr(0,alnLength);
				temp = temp.substr(0,alnLength);
				
				int numDiff = countDiffs(oligo, temp);
				
				if(numDiff < minDiff){
					minDiff = numDiff;
					minCount = 1;
					minGroup = it->second;
					minPos = 0;
					for(int i=0;i<alnLength;i++){
						if(temp[i] != '-'){
							minPos++;
						}
					}
				}
				else if(numDiff == minDiff){
					minCount++;
				}
				
			}
			
			if(minDiff > bdiffs)	{	success = minDiff;		}	//no good matches
			else if(minCount > 1)	{	success = bdiffs + 100;	}	//can't tell the difference between multiple barcodes
			else{													//use the best match
//.........这里部分代码省略.........