showalign.cpp

来自「ncbi源码」· C++ 代码 · 共 1,833 行 · 第 1/5 页

	int featSeqFrom = loc.GetInt().GetFrom();
	int featSeqTo = loc.GetInt().GetTo();
	int actualFeatSeqStart = 0, actualFeatSeqStop = 0;
	if(m_AV->IsPositiveStrand(row)){
	  if(featSeqFrom < (int)m_AV->GetSeqStart(row)){
	    actualFeatSeqStart = m_AV->GetSeqStart(row);
	  } else {
	    actualFeatSeqStart = featSeqFrom;
	  }

	  if(featSeqTo > (int)m_AV->GetSeqStop(row)){
	    actualFeatSeqStop = m_AV->GetSeqStop(row);
	  } else {
	    actualFeatSeqStop = featSeqTo;
	  }
	} else {
	  if(featSeqFrom < (int)m_AV->GetSeqStart(row)){
	    actualFeatSeqStart = featSeqFrom;
	  } else {
	    actualFeatSeqStart = m_AV->GetSeqStart(row); 
	  }

	  if(featSeqTo > (int)m_AV->GetSeqStop(row)){
	    actualFeatSeqStop = featSeqTo;
	  } else {
	    actualFeatSeqStop = m_AV->GetSeqStop(row);
	  }
	}

	int alnFrom = m_AV->GetAlnPosFromSeqPos(row, actualFeatSeqStart);
	int alnTo = m_AV->GetAlnPosFromSeqPos(row, actualFeatSeqStop);
	char featChar = ' ';
	if(choice == CSeqFeatData::e_Gene){
	  featChar = '^';
	} else if (choice == CSeqFeatData::e_Cdregion){
	  featChar = '~';
	}

	//need to construct the protein seq aligned to nucleotide seq 
	if (choice == CSeqFeatData::e_Cdregion){
	  string rawCdrProduct = NcbiEmptyString;
	  if(feat->IsSetProduct()){
	    const CSeq_loc& productLoc = feat->GetProduct(); 
	    //only show first k_FeatureIdLen letters
	    if(productLoc.IsWhole()){
	      const CSeq_id& productId = productLoc.GetWhole();
	      const CBioseq_Handle& productHandle = scope.GetBioseqHandle(productId );
	      featId = "CDS:" + GetTitle(productHandle).substr(0, k_FeatureIdLen);
	    }
	  }
	  //show protein product only if the row is plus strand
	  if(m_AV->IsPositiveStrand(row)){
	    string line(alnStop+1, ' ');  
	    CCdregion_translate::TranslateCdregion (rawCdrProduct, handle, loc, feat->GetData().GetCdregion(), true, false);
	    
	    bool firstBase = true;
	    char gapChar = m_AV->GetGapChar(row);
	    int marginAdjuster = 0;
	    int featStartSeqPos = 0;
	    int firstFeatStringPos = 0;
	    int numBase = 0;
	    
	    //put actual amino acid to the cdr product line in aln coord
	    for (int i = alnFrom; i < alnTo; i++){
	      if(sequence[i] != gapChar){
		numBase ++;
		if(firstBase){
		  firstBase = false;
		  featStartSeqPos = m_AV->GetSeqPosFromAlnPos(row, i);
		  const CCdregion& cdr = feat->GetData().GetCdregion();
		  int frame = 1;
		  if(cdr.IsSetFrame()){
		    frame = cdr.GetFrame();
		  }
		  int numBaseFromFeatStart = (featStartSeqPos - (featSeqFrom + (frame -1) )); //Number of bases between feature start and current base. adjust using frame
		  if(numBaseFromFeatStart % 3 == 0){ //this base is the 1st base
		    marginAdjuster = 0; //aa aligned to 2nd base of a condon
		    firstFeatStringPos = numBaseFromFeatStart / 3;  
		  } else if (numBaseFromFeatStart % 3 == 1) {
		    marginAdjuster = 1;
		    firstFeatStringPos = numBaseFromFeatStart / 3;
		  } else {
		    marginAdjuster = -1;
		    firstFeatStringPos = (numBaseFromFeatStart / 3) + 1;
		  }
		}
		if((numBase + marginAdjuster) % 3 == 2){
		  int stringPos = firstFeatStringPos + (numBase + marginAdjuster) / 3;
		  if(stringPos < (int)rawCdrProduct.size()){//should not need this check
		    line[i] = rawCdrProduct[stringPos];
		  }
		}
	      }
	    }
	    alternativeFeatStr = line;
	  }
	} 

	setFeatureInfo(featInfo, loc, alnFrom, alnTo, alnStop, featChar, featId, alternativeFeatStr);     
	feature.push_back(featInfo);
      
      }
    }
  }
}

void  CDisplaySeqalign::setFeatureInfo(alnFeatureInfo* featInfo, const CSeq_loc& seqloc, int alnFrom, int alnTo, int alnStop, char patternChar, string patternId, string& alternativeFeatStr) const{
  FeatureInfo* feat = new FeatureInfo;
  feat->seqloc = &seqloc;
  feat->featureChar = patternChar;
  feat->featureId = patternId;
    
  if(alternativeFeatStr != NcbiEmptyString){
    featInfo->featureString = alternativeFeatStr;
  } else {
    //fill feature string
    string line(alnStop+1, ' ');
    for (int j = alnFrom; j <= alnTo; j++){
      line[j] = feat->featureChar;
    }
    featInfo->featureString = line;
  }
 
  featInfo->alnRange.Set(alnFrom, alnTo); 
  featInfo->feature = feat;
}

//May need to add a "|" to the current insert for insert on next rows
static int addBar(string& seq, int insertAlnPos, int alnStart){
  int end = seq.size() -1 ;
  int barPos = insertAlnPos - alnStart + 1;
  string addOn;
  if(barPos - end > 1){
    string spacer(barPos - end - 1, ' ');
    addOn += spacer + "|";
  } else if (barPos - end == 1){
    addOn += "|";
  }
  seq += addOn;
  return max<int>((barPos - end), 0);
}

//Add new insert seq to the current insert seq and return the end position of the latest insert
static int adjustInsert(string& curInsert, string& newInsert, int insertAlnPos, int alnStart){
  int insertEnd = 0;
  int curInsertSize = curInsert.size();
  int insertLeftSpace = insertAlnPos - alnStart - curInsertSize + 2;  //plus2 because insert is put after the position
  if(curInsertSize > 0){
    assert(insertLeftSpace >= 2);
  }
  int newInsertSize = newInsert.size();  
  if(insertLeftSpace - newInsertSize >= 1){ //can insert with the end position right below the bar
    string spacer(insertLeftSpace - newInsertSize, ' ');
    curInsert += spacer + newInsert;
    
  } else { //Need to insert beyond the insert postion
    if(curInsertSize > 0){
      curInsert += " " + newInsert;
    } else {  //can insert right at the firt position
      curInsert += newInsert;
    }
  }
  insertEnd = alnStart + curInsert.size() -1 ; //-1 back to string position
  return insertEnd;
}

//recusively fill the insert
void CDisplaySeqalign::doFills(int row, CAlnMap::TSignedRange& alnRange, int  alnStart, list<insertInformation*>& insertList, list<string>& inserts) const {
  if(!insertList.empty()){
    string bar(alnRange.GetLength(), ' ');
    
    string seq;
    list<insertInformation*> leftOverInsertList;
    bool isFirstInsert = true;
    int curInsertAlnStart = 0;
    int prvsInsertAlnEnd = 0;
    
    //go through each insert and fills the seq if it can  be filled on the same line.  If not, go to the next line
    for(list<insertInformation*>::iterator iter = insertList.begin(); iter != insertList.end(); iter ++){
      curInsertAlnStart = (*iter)->alnStart;
      //always fill the first insert.  Also fill if there is enough space
      if(isFirstInsert || curInsertAlnStart - prvsInsertAlnEnd >= 1){
	bar[curInsertAlnStart-alnStart+1] = '|';  
	int seqStart = (*iter)->seqStart;
	int seqEnd = seqStart + (*iter)->insertLen - 1;
	string newInsert;
	newInsert = m_AV->GetSeqString(newInsert, row, seqStart, seqEnd);
	prvsInsertAlnEnd = adjustInsert(seq, newInsert, curInsertAlnStart, alnStart);
	isFirstInsert = false;
      } else { //if no space, save the chunk and go to next line 
	bar[curInsertAlnStart-alnStart+1] = '|';  //indicate insert goes to the next line
	prvsInsertAlnEnd += addBar(seq, curInsertAlnStart, alnStart);   //May need to add a bar after the current insert sequence to indicate insert goes to the next line.
	leftOverInsertList.push_back(*iter);    
      }
    }
    //save current insert.  Note that each insert has a bar and sequence below it
    inserts.push_back(bar);
    inserts.push_back(seq);
    //here recursively fill the chunk that don't have enough space
    doFills(row, alnRange, alnStart, leftOverInsertList, inserts);
  }
 
}

/*fill a list of inserts for a particular row*/
void CDisplaySeqalign::fillInserts(int row, CAlnMap::TSignedRange& alnRange, int alnStart, list<string>& inserts, string& insertPosString, list<insertInformation*>& insertList) const{

  string line(alnRange.GetLength(), ' ');
 
  ITERATE(list<insertInformation*>, iter, insertList){
    int from = (*iter)->alnStart;
    line[from - alnStart + 1] = '\\';
  }
  insertPosString = line; //this is the line with "\" right after each insert position
    
  //here fills the insert sequence
  doFills(row, alnRange, alnStart, insertList, inserts);
}

void CDisplaySeqalign::GetInserts(list<insertInformation*>& insertList, CAlnMap::TSeqPosList& insertAlnStart, CAlnMap::TSeqPosList& insertSeqStart, CAlnMap::TSeqPosList& insertLength, int lineAlnStop){

  while(!insertAlnStart.empty() && (int)insertAlnStart.front() < lineAlnStop){
    CDisplaySeqalign::insertInformation* insert = new CDisplaySeqalign::insertInformation;
    insert->alnStart = insertAlnStart.front() - 1; //Need to minus one as we are inserting after this position
    insert->seqStart = insertSeqStart.front();
    insert->insertLen = insertLength.front();
    insertList.push_back(insert);
    insertAlnStart.pop_front();
    insertSeqStart.pop_front();
    insertLength.pop_front();
  }
 
}


//segments starts and stops used for map viewer 
string CDisplaySeqalign::getSegs(int row) const {
  string segs = NcbiEmptyString;
  if(m_AlignOption & eMultiAlign){ //only show this hsp
    segs = NStr::IntToString(m_AV->GetSeqStart(row)) + "-" + NStr::IntToString(m_AV->GetSeqStop(row));
  } else { //for all segs
    string idString = m_AV->GetSeqId(1).GetSeqIdString();
    map<string, string>::const_iterator iter = m_Segs.find(idString);
    if ( iter != m_Segs.end() ){
      segs = iter->second;
    }
  }
  return segs;
}



/* transforms a string so that it becomes safe to be used as part of URL
 * the function converts characters with special meaning (such as
 * semicolon -- protocol separator) to escaped hexadecimal (%xx)
 */
static string MakeURLSafe(char* src){
  static char HEXDIGS[] = "0123456789ABCDEF";
  char* buf;
  size_t len;
  char* p;
  char c;
  string url = NcbiEmptyString;

  if (src){
    /* first pass to calculate required buffer size */
    for (p = src, len = 0; (c = *(p++)) != '\0'; ) {
      switch (c) {
      default:
	if (c < '0' || (c > '9' && c < 'A') ||
	    (c > 'Z' && c < 'a') || c > 'z') {
	  len += 3;
	  break;
	}
      case '-': case '_': case '.': case '!': case '~':
      case '*': case '\'': case '(': case ')':
	++len;
      }
    }
    buf = new char[len + 1];
    /* second pass -- conversion */
    for (p = buf; (c = *(src++)) != '\0'; ) {
      switch (c) {
      default:
	if (c < '0' || (c > '9' && c < 'A') ||
	    (c > 'Z' && c < 'a') || c > 'z') {
	  *(p++) = '%';
	  *(p++) = HEXDIGS[(c >> 4) & 0xf];
	  *(p++) = HEXDIGS[c & 0xf];
	  break;
	}
      case '-': case '_': case '.': case '!': case '~':
      case '*': case '\'': case '(': case ')':
	*(p++) = c;
      }
    }
    *p = '\0';
    url = buf;
    delete [] buf;
  }
  return url;
}

//make url for dumpgnl.cgi
string CDisplaySeqalign::getDumpgnlLink(const list<CRef<CSeq_id> >& ids, int row, const string& alternativeUrl)const {
  string link = NcbiEmptyString;  
  string toolUrl= m_Reg->Get(m_BlastType, "TOOL_URL");
  string passwd = m_Reg->Get(m_BlastType, "PASSWD");
  bool nodb_path =  false;
  CRef<CSeq_id> idGeneral = GetSeqIdByType(ids, CSeq_id::e_General);
  CRef<CSeq_id> idOther = GetSeqIdByType(ids, CSeq_id::e_Other);
  const CRef<CSeq_id> idAccession = FindBestChoice(ids, CSeq_id::WorstRank);
  string segs = getSegs(row);
  int gi = GetGiForSeqIdList(ids);
  if(!idGeneral.Empty() && idGeneral->AsFastaString().find("gnl|BL_ORD_ID")){ /* We do need to make security protected link to BLAST gnl */
    return NcbiEmptyString;
}
  if(alternativeUrl != NcbiEmptyString){ 
    toolUrl = alternativeUrl;
  }
  /* dumpgnl.cgi need to use path  */
  if (toolUrl.find("dumpgnl.cgi") ==string::npos){
    nodb_path = true;
  }  
  int length = m_DbName.size();
  string str;
  char  *chptr, *dbtmp;
  Char tmpbuff[256];
  char* dbname = new char[sizeof(char)*length + 2];
  strcpy(dbname, m_DbName.c_str());
  if(nodb_path) {
    int i, j;
    dbtmp = new char[sizeof(char)*length + 2]; /* aditional space and NULL */
    memset(dbtmp, '\0', sizeof(char)*length + 2);
    for(i = 0; i < length; i++) {            
       if(isspace(dbname[i]) || dbname[i] == ',') {/* Rolling spaces */
	 continue;
       }
       j = 0;
       while (!isspace(dbname[i]) && j < 256  && i < length) { 
	 tmpbuff[j] = dbname[i];
	 j++; i++;
	 if(dbname[i] == ',') { /* Comma is valid delimiter */
	   break;
	 }
       }
       tmpbuff[j] = '\0';
       if((chptr = strrchr(tmpbuff, '/')) != NULL) { 
	 strcat(dbtmp, (char*)(chptr+1));
       } else {
	 strcat(dbtmp, tmpbuff);
       }
       strcat(dbtmp, " ");            
     }
   } else {
     dbtmp = dbname;
   }
  
  const CSeq_id* bestid = NULL;
  if (idGeneral.Empty()){
    bestid = idOther;
    if (idOther.Empty()){
      bestid = idAccession;
    }