showalign.cpp

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

  }  //end of preparing row data
  
  bool colorMismatch = false; //color the mismatches
  //output identities info 
  if(m_AlignOption&eShowBlastInfo && !(m_AlignOption&eMultiAlign)) {
    int match = 0;
    int positive = 0;
    int gap = 0;
    int identity = 0;
    fillIdentityInfo(sequence[0], sequence[1],  match,  positive, middleLine);
    identity = (match*100)/(aln_stop+1);
    if(identity >= k_ColorMismatchIdentity && identity <100){
      colorMismatch = true;
    }
    out<<" Identities = "<<match<<"/"<<(aln_stop+1)<<" ("<<identity<<"%"<<")";
    if(m_AlignType&eProt) {
      out<<", Positives = "<<(positive + match)<<"/"<<(aln_stop+1)<<" ("<<(((positive + match)*100)/(aln_stop+1))<<"%"<<")";
    }
    gap = getNumGaps();
    out<<", Gaps = "<<gap<<"/"<<(aln_stop+1)<<" ("<<((gap*100)/(aln_stop+1))<<"%"<<")"<<endl;
    if (m_AlignType&eNuc){ 
      out<<" Strand="<<(m_AV->StrandSign(0)==1 ? "Plus" : "Minus")<<"/"<<(m_AV->StrandSign(1)==1? "Plus" : "Minus")<<endl;
    }
    if(frame[0] != 0 && frame[1] != 0) {
      out <<" Frame = " << ((frame[0] > 0) ? "+" : "") << frame[0] <<"/"<<((frame[1] > 0) ? "+" : "") << frame[1]<<endl;
    } else if (frame[0] != 0){
      out <<" Frame = " << ((frame[0] > 0) ? "+" : "") << frame[0] <<endl;
    }  else if (frame[1] != 0){
      out <<" Frame = " << ((frame[1] > 0) ? "+" : "") << frame[1] <<endl;
    } 
    out<<endl;
  }
  //output rows
  for(int j=0; j<=aln_stop; j+=m_LineLen){
    //output according to aln coordinates
    if(aln_stop-j+1<m_LineLen) {
      actualLineLen=aln_stop-j+1;
    } else {
      actualLineLen=m_LineLen;
    }
    CAlnMap::TSignedRange curRange(j, j+actualLineLen-1);
    //here is each row
    for (int row=0; row<rowNum; row++) {
      bool hasSequence = true;
   
      hasSequence = curRange.IntersectingWith(rowRng[row]);
     
      //only output rows that have sequence
      if (hasSequence){
	int start = seqStarts[row].front() + 1;  //+1 for 1 based
	int end = seqStops[row].front() + 1;
	list<string> inserts;
	string insertPosString;  //the one with "\" to indicate insert
	if(m_AlignOption & eMasterAnchored){
	  list<insertInformation*> insertList;
	  GetInserts(insertList, insertAlnStart[row], insertStart[row], insertLength[row],  j + m_LineLen);
	  fillInserts(row, curRange, j, inserts, insertPosString, insertList);
	  ITERATE(list<insertInformation*>, iterINsert, insertList){
	    delete *iterINsert;
	  }
	}
        if(row == 0&&(m_AlignOption&eHtml)&&(m_AlignOption&eMultiAlign) && (m_AlignOption&eSequenceRetrieval && m_IsDbGi)){
          char checkboxBuf[200];
          sprintf(checkboxBuf, "<input type=\"checkbox\" name=\"getSeqMaster\" value=\"\" onClick=\"uncheckable('getSeqAlignment%d', 'getSeqMaster')\">", m_QueryNumber);
          out << checkboxBuf;
        }
        string urlLink;
        //setup url link for seqid
        if(row>0&&(m_AlignOption&eHtml)&&(m_AlignOption&eMultiAlign)){
        
          int gi = GetGiForSeqIdList(m_AV->GetBioseqHandle(row).GetBioseqCore()->GetId());
          if(gi > 0){
            out<<"<a name="<<gi<<"></a>";
          } else {
            out<<"<a name="<<seqidArray[row]<<"></a>";
          }
          //get sequence checkbox
          if(m_AlignOption&eSequenceRetrieval && m_IsDbGi){
            char checkBoxBuf[512];
            sprintf(checkBoxBuf, "<input type=\"checkbox\" name=\"getSeqGi\" value=\"%d\" onClick=\"synchronizeCheck(this.value, 'getSeqAlignment%d', 'getSeqGi', this.checked)\">", gi, m_QueryNumber);
            out << checkBoxBuf;        
          }
          urlLink = getUrl(m_AV->GetBioseqHandle(row).GetBioseqCore()->GetId(), row);         
          out << urlLink;
         
        }
        out<<seqidArray[row]; 
        if(row>0&& m_AlignOption&eHtml && m_AlignOption&eMultiAlign && urlLink != NcbiEmptyString){
          out<<"</a>";
         
        }
        //adjust space between id and start
        AddSpace(out, maxIdLen-seqidArray[row].size()+m_IdStartMargin);
	out << start;
        startLen=NStr::IntToString(start).size();
        AddSpace(out, maxStartLen-startLen+m_StartSequenceMargin);
	if (row>0 && m_AlignOption & eShowIdentity){
	  for (int index = j; index < j + actualLineLen && index < (int)sequence[row].size(); index ++){
	    if (sequence[row][index] == sequence[0][index] && isalpha(sequence[row][index])) {
	      sequence[row][index] = k_IdentityChar;           
	    }         
	  }
	}
	OutputSeq(sequence[row], m_AV->GetSeqId(row), j, actualLineLen, frame[row], (row > 0 && colorMismatch)?true:false, out);
        AddSpace(out, m_SeqStopMargin);
	out << end;
        
        out<<endl;
     
	//display inserts for anchored type
	if(m_AlignOption & eMasterAnchored){
	  bool insertAlready = false;
	  for(list<string>::iterator iter = inserts.begin(); iter != inserts.end(); iter ++){
	   
	    if(!insertAlready){
	      if((m_AlignOption&eHtml)&&(m_AlignOption&eMultiAlign) && (m_AlignOption&eSequenceRetrieval && m_IsDbGi)){
		char checkboxBuf[200];
		sprintf(checkboxBuf, "<input type=\"checkbox\" name=\"getSeqMaster\" value=\"\" onClick=\"uncheckable('getSeqAlignment%d', 'getSeqMaster')\">", m_QueryNumber);
		out << checkboxBuf;
	      }
	      AddSpace(out, maxIdLen+m_IdStartMargin+maxStartLen+m_StartSequenceMargin);
	      out << insertPosString<<endl;
	    }
	    if((m_AlignOption&eHtml)&&(m_AlignOption&eMultiAlign) && (m_AlignOption&eSequenceRetrieval && m_IsDbGi)){
	      char checkboxBuf[200];
	      sprintf(checkboxBuf, "<input type=\"checkbox\" name=\"getSeqMaster\" value=\"\" onClick=\"uncheckable('getSeqAlignment%d', 'getSeqMaster')\">", m_QueryNumber);
	      out << checkboxBuf;
	    }
	    AddSpace(out, maxIdLen+m_IdStartMargin+maxStartLen+m_StartSequenceMargin);
	    out<<*iter<<endl;
	    
	    insertAlready = true;
	  }
	} 
	//display feature. Feature, if set, will be displayed for query regardless

        CSeq_id no_id;
	for (list<alnFeatureInfo*>::iterator iter=bioseqFeature[row].begin();  iter != bioseqFeature[row].end(); iter++){
	  if ( curRange.IntersectingWith((*iter)->alnRange)){  
	    if((m_AlignOption&eHtml)&&(m_AlignOption&eMultiAlign) && (m_AlignOption&eSequenceRetrieval && m_IsDbGi)){
	      char checkboxBuf[200];
	      sprintf(checkboxBuf, "<input type=\"checkbox\" name=\"getSeqMaster\" value=\"\" onClick=\"uncheckable('getSeqAlignment%d', 'getSeqMaster')\">", m_QueryNumber);
	      out << checkboxBuf;
	    }
	    out<<(*iter)->feature->featureId;
	    AddSpace(out, maxIdLen+m_IdStartMargin+maxStartLen+m_StartSequenceMargin-(*iter)->feature->featureId.size());
	    OutputSeq((*iter)->featureString, no_id, j, actualLineLen, 0, false, out);
	    out<<endl;
	  }
	}
	
	//display middle line
	if (row == 0 && ((m_AlignOption & eShowMiddleLine)) && !(m_AlignOption&eMultiAlign)) {
	  AddSpace(out, maxIdLen+m_IdStartMargin+maxStartLen+m_StartSequenceMargin);
	  OutputSeq(middleLine, no_id, j, actualLineLen, 0, false, out);
	  out<<endl;
	}
      }
      if(!seqStarts[row].empty()){ //shouldn't need this check
	seqStarts[row].pop_front();
      }
      if(!seqStops[row].empty()){
	seqStops[row].pop_front();
      }
    }
    out<<endl;
  }//end of displaying rows
 
  //free allocation
  for(int i = 0; i < rowNum; i ++){
    for (list<alnFeatureInfo*>::iterator iter=bioseqFeature[i].begin();  iter != bioseqFeature[i].end(); iter++){
      delete (*iter)->feature;
      delete (*iter);
    }
  } 

  for (list<alnSeqlocInfo*>::const_iterator iter = alnLocList.begin();  iter != alnLocList.end(); iter++){
    delete (*iter);
  }

  delete [] bioseqFeature;
  delete [] seqidArray;
  delete [] rowRng;
  delete [] seqStarts;
  delete [] seqStops;
  delete [] frame;
  delete [] insertStart;
  delete [] insertAlnStart;
  delete [] insertLength;
}

//To display the seqalign
void CDisplaySeqalign::DisplaySeqalign(CNcbiOstream& out){
  
  if(m_SeqalignSetRef->Get().empty()){
    return;
  }

  //scope for feature fetching
  if(!(m_AlignOption & eMasterAnchored) && (m_AlignOption & eShowCdsFeature || m_AlignOption & eShowGeneFeature)){
    m_FeatObj = new CObjectManager();
    m_FeatObj->RegisterDataLoader(*new CGBDataLoader("ID", NULL, 2), CObjectManager::eDefault); 
    m_featScope = new CScope(*m_FeatObj);  //for seq feature fetch
    m_featScope->AddDefaults();	     
  }	
 
  setDbGi(); //for whether to add get sequence feature
  if(m_AlignOption & eHtml){
    //set config file
    m_ConfigFile = new CNcbiIfstream(".ncbirc");
    m_Reg = new CNcbiRegistry(*m_ConfigFile);
    out<<"<script src=\"blastResult.js\"></script>";
  }
   //get sequence 
  if(m_AlignOption&eSequenceRetrieval && m_AlignOption&eHtml && m_IsDbGi){ 
        out<<GetSeqForm((char*)"submitterTop", m_IsDbNa, m_QueryNumber);
        out<<"<form name=\"getSeqAlignment"<<m_QueryNumber<<"\">\n";
      }

  //begin to display
  int num_align = 0;
  string toolUrl = NcbiEmptyString;
  if(m_AlignOption & eHtml){
    toolUrl = m_Reg->Get(m_BlastType, "TOOL_URL");
  }
  auto_ptr<CObjectOStream> out2(CObjectOStream::Open(eSerial_AsnText, out));
  //*out2 << *m_SeqalignSetRef;
  if(!(m_AlignOption&eMultiAlign)){/*pairwise alignment. Note we can't just show each alnment as we go because we will need seg information form all hsp's with the same id for genome url link.  As a result we show hsp's with the same id as a group*/
    list<alnInfo*> avList;        
    CConstRef<CSeq_id> previousId, subid;
    bool isFirstAln = true;
    for (CSeq_align_set::Tdata::const_iterator iter = m_SeqalignSetRef->Get().begin(); iter != m_SeqalignSetRef->Get().end()&&num_align<m_NumAlignToShow; iter++, num_align++) {
      
      //make alnvector
      CRef<CAlnVec> avRef;
      CRef<CSeq_align> finalAln;
      if((*iter)->GetSegs().Which() == CSeq_align::C_Segs::e_Std){
	CRef<CSeq_align> densegAln = (*iter)->CreateDensegFromStdseg();
	if (m_AlignOption & eTranslateNucToNucAlignment) { 
	  finalAln = densegAln->CreateTranslatedDensegFromNADenseg();
	} else {
	  finalAln = densegAln;
	}

      } else if((*iter)->GetSegs().Which() == CSeq_align::C_Segs::e_Denseg){
	if (m_AlignOption & eTranslateNucToNucAlignment) { 
	  finalAln = (*iter)->CreateTranslatedDensegFromNADenseg();
	} else {
	  finalAln = (*iter);
	}

      } else if((*iter)->GetSegs().Which() == CSeq_align::C_Segs::e_Dendiag){
	CRef<CSeq_align> densegAln = CreateDensegFromDendiag(**iter);
	if (m_AlignOption & eTranslateNucToNucAlignment) { 
	  finalAln = densegAln->CreateTranslatedDensegFromNADenseg();
	} else {
	  finalAln = densegAln;
	}

      } else {
	NCBI_THROW(CException, eUnknown, "Seq-align should be Denseg, Stdseg or Dendiag!");
      }
      CRef<CDense_seg> finalDenseg(new CDense_seg);
      const CTypeIterator<CDense_seg> ds = Begin(*finalAln);
      if((ds->IsSetStrands() && ds->GetStrands().front()==eNa_strand_minus) && !(ds->IsSetWidths() && ds->GetWidths()[0] == 3)){
	//show plus strand if master is minus for non-translated case
	memcpy(&*finalDenseg, &(*ds), sizeof(CDense_seg));
	finalDenseg->Reverse();
	avRef = new CAlnVec(*finalDenseg, m_Scope);
	
      } else {
	avRef = new CAlnVec(*ds, m_Scope);
      }
    
      if(!(avRef.Empty())){
	try{
	  const CBioseq_Handle& handle = avRef->GetBioseqHandle(1);	
	  if(handle){
	    subid=&(avRef->GetSeqId(1));
	    
	    if(!isFirstAln && !subid->Match(*previousId)) {//this aln is a new id, show result for previous id
	      x_DisplayAlnvecList(out, avList);
	    
	      for(list<alnInfo*>::iterator iterAv = avList.begin(); iterAv != avList.end(); iterAv ++){
		delete(*iterAv);
	      }
	      avList.clear();
	      
	    }
	    //save the current alnment regardless
	    alnInfo* alnvecInfo = new alnInfo;
	    getAlnScores(**iter, alnvecInfo->score, alnvecInfo->bits, alnvecInfo->eValue);
	    alnvecInfo->alnVec = avRef;
	    avList.push_back(alnvecInfo);
	    int gi = GetGiForSeqIdList(handle.GetBioseqCore()->GetId());
	    if(!(toolUrl == NcbiEmptyString || (gi > 0 && toolUrl.find("dumpgnl.cgi") != string::npos)) || (m_AlignOption & eLinkout)){ /*need to construct segs for dumpgnl and get sub-sequence for long sequences*/
	      string idString = avRef->GetSeqId(1).GetSeqIdString();
	      if(m_Segs.count(idString) > 0){ 	//already has seg, concatenate
		/*Note that currently it's not necessary to use map to store this information.  But I already implemented this way for previous version.  Will keep this way as it's more flexible if we change something*/
	
		m_Segs[idString] += "," + NStr::IntToString(avRef->GetSeqStart(1)) + "-" + NStr::IntToString(avRef->GetSeqStop(1));
	      } else {//new segs
		m_Segs.insert(map<string, string>::value_type(idString, NStr::IntToString(avRef->GetSeqStart(1)) + "-" + NStr::IntToString(avRef->GetSeqStop(1))));
	      }
	    }	    
	    isFirstAln = false;
	    previousId = subid;
	  }
	 
	} catch (CException& e){
	continue;
	}
	
      }
    } 
  
    //Show here for the last one 
    if(!avList.empty()){
      x_DisplayAlnvecList(out, avList);
      for(list<alnInfo*>::iterator iterAv = avList.begin(); iterAv != avList.end(); iterAv ++){
	delete(*iterAv);
      }
      avList.clear();
    }
  	     
  } else if(m_AlignOption&eMultiAlign){ //multiple alignment
       
    CRef<CAlnMix>* mix = new CRef<CAlnMix>[k_NumFrame]; //each for one frame for translated alignment
    for(int i = 0; i < k_NumFrame; i++){
      mix[i] = new CAlnMix(m_Scope);
    }

    num_align = 0;
    vector<CRef<CSeq_align_set> > alnVector(k_NumFrame);
    for(int i = 0; i <  k_NumFrame; i ++){
      alnVector[i] = new CSeq_align_set;
    }
    for (CSeq_align_set::Tdata::const_iterator alnIter = m_SeqalignSetRef->Get().begin(); alnIter != m_SeqalignSetRef->Get().end()&&num_align<m_NumAlignToShow; alnIter ++, num_align++) {
      //need to convert to denseg for stdseg
      if((*alnIter)->GetSegs().Which() == CSeq_align::C_Segs::e_Std) {
	CTypeConstIterator<CStd_seg> ss = ConstBegin(**alnIter); 
	CRef<CSeq_align> convertedDs = (*alnIter)->CreateDensegFromStdseg();
	if((convertedDs->GetSegs().GetDenseg().IsSetWidths() && convertedDs->GetSegs().GetDenseg().GetWidths()[0] == 3) || m_AlignOption & eTranslateNucToNucAlignment){//only do this for translated master
	  int frame = s_GetStdsegMasterFrame(*ss, m_Scope);
	  switch(frame){
	  case 1:
	    alnVector[0]->Set().push_back(convertedDs);
	    break;
	  case 2:
	    alnVector[1]->Set().push_back(convertedDs);
	    break;
	  case 3:
	    alnVector[2]->Set().push_back(convertedDs);
	    break;
	  case -1:
	    alnVector[3]->Set().push_back(convertedDs);
	    break;
	  case -2:
	    alnVector[4]->Set().push_back(convertedDs);
	    break;
	  case -3:
	    alnVector[5]->Set().push_back(convertedDs);
	    break;
	  default:
	    break;
	  }
	}