showalign.cpp

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

	else {
	  alnVector[0]->Set().push_back(convertedDs);
	}
      } else if((*alnIter)->GetSegs().Which() == CSeq_align::C_Segs::e_Denseg){
	alnVector[0]->Set().push_back(*alnIter);
      } else if((*alnIter)->GetSegs().Which() == CSeq_align::C_Segs::e_Dendiag){
	alnVector[0]->Set().push_back(CreateDensegFromDendiag(**alnIter));
      } else {
	NCBI_THROW(CException, eUnknown, "Input Seq-align should be Denseg, Stdseg or Dendiag!");
      }
    }
    
    for(int i = 0; i < (int)alnVector.size(); i ++){
      bool hasAln = false;
      for(CTypeConstIterator<CSeq_align> alnRef = ConstBegin(*alnVector[i]); alnRef; ++alnRef){
	CTypeConstIterator<CDense_seg> ds = ConstBegin(*alnRef);
        //*out2 << *ds;      
	try{
	  if (m_AlignOption & eTranslateNucToNucAlignment) {	 
	    mix[i]->Add(*ds, CAlnMix::fForceTranslation);
	  } else {
	    mix[i]->Add(*ds);
	  }
	} catch (CException& e){
	  continue;
	}
	 hasAln = true;
      }
      if(hasAln){
        //    *out2<<*alnVector[i];
	  mix[i]->Merge(CAlnMix::fGen2EST| CAlnMix::fMinGap | CAlnMix::fQuerySeqMergeOnly | CAlnMix::fFillUnalignedRegions);  
	//	*out2<<mix[i]->GetDenseg();
      }
    }
    
    int numDistinctFrames = 0;
    for(int i = 0; i < (int)alnVector.size(); i ++){
      if(!alnVector[i]->Get().empty()){
	numDistinctFrames ++;
      }
    }
    
    out<<endl;
    for(int i = 0; i < k_NumFrame; i ++){
      try{
	CRef<CAlnVec> avRef (new CAlnVec (mix[i]->GetDenseg(), m_Scope));
	m_AV = avRef;
	if(numDistinctFrames > 1){
	 out << "For reading frame " << k_FrameConversion[i] << " of query sequence:" << endl << endl;
	}
	DisplayAlnvec(out);
      } catch (CException e){
	continue;
      }
    } 
    delete [] mix;
  }
  if(m_AlignOption&eSequenceRetrieval && m_AlignOption&eHtml && m_IsDbGi){
    out<<"</form>\n";
    out<<GetSeqForm((char*)"submitterBottom", m_IsDbNa, m_QueryNumber);
  }
}

//compute number of identical and positive residues; set middle line accordingly
const void CDisplaySeqalign::fillIdentityInfo(const string& sequenceStandard, const string& sequence , int& match, int& positive, string& middleLine) {
  match = 0;
  positive = 0;
  int min_length=min<int>(sequenceStandard.size(), sequence.size());
  if(m_AlignOption & eShowMiddleLine){
    middleLine = sequence;
  }
  for(int i=0; i<min_length; i++){
    if(sequenceStandard[i]==sequence[i]){
      if(m_AlignOption & eShowMiddleLine){
	if(m_MidLineStyle == eBar ) {
	  middleLine[i] = '|';
	} else if (m_MidLineStyle == eChar){
	  middleLine[i] = sequence[i];
	}
      }
      match ++;
    } else {
      if ((m_AlignType&eProt) && m_Matrix[sequenceStandard[i]][sequence[i]] > 0){  
	positive ++;
	if(m_AlignOption & eShowMiddleLine){
	  if (m_MidLineStyle == eChar){
	    middleLine[i] = '+';
	  }
	}
      } else {
	if (m_AlignOption & eShowMiddleLine){
	  middleLine[i] = ' ';
	}
      }    
    }
  }  
}


const void CDisplaySeqalign::PrintDefLine(const CBioseq_Handle& bspHandle, CNcbiOstream& out) const
{
  if(bspHandle){
    const CRef<CSeq_id> wid = FindBestChoice(bspHandle.GetBioseqCore()->GetId(), CSeq_id::WorstRank);
 
    const CRef<CBlast_def_line_set> bdlRef = GetBlastDefline(bspHandle);
    const list< CRef< CBlast_def_line > >& bdl = bdlRef->Get();
    bool isFirst = true;
    int firstGi = 0;
  
    if(bdl.empty()){ //no blast defline struct, should be no such case now
      out << ">"; 
      wid->WriteAsFasta(out);
      out<<" ";
      s_WrapOutputLine(out, GetTitle(bspHandle));
      out << endl;
    } else {
      //print each defline 
      for(list< CRef< CBlast_def_line > >::const_iterator iter = bdl.begin(); iter != bdl.end(); iter++){
	string urlLink;
	if(isFirst){
	  out << ">";
   
	} else{
	  out << " ";
	}
	const CRef<CSeq_id> wid2 = FindBestChoice((*iter)->GetSeqid(), CSeq_id::WorstRank);
	int gi =  GetGiForSeqIdList((*iter)->GetSeqid());
	if(isFirst){
	  firstGi = gi;
	}
	if ((m_AlignOption&eSequenceRetrieval) && (m_AlignOption&eHtml) && m_IsDbGi && isFirst) {
	  char buf[512];
	  sprintf(buf, "<input type=\"checkbox\" name=\"getSeqGi\" value=\"%d\" onClick=\"synchronizeCheck(this.value, 'getSeqAlignment%d', 'getSeqGi', this.checked)\">", gi, m_QueryNumber);
	  out << buf;
	}
 
	if(m_AlignOption&eHtml){
     
	  urlLink = getUrl((*iter)->GetSeqid(), 1);    
	  out<<urlLink;
	}
    
	if(m_AlignOption&eShowGi && gi > 0){
	  out<<"gi|"<<gi<<"|";
	}       
   
	wid2->WriteAsFasta(out);
	if(m_AlignOption&eHtml){
	  if(urlLink != NcbiEmptyString){
	    out<<"</a>";
	  }
	  if(gi != 0){
	    out<<"<a name="<<gi<<"></a>";
	  } else {
	    out<<"<a name="<<wid2->GetSeqIdString()<<"></a>";
	  }
	  if(m_AlignOption&eLinkout){
	    out <<" ";
	    AddLinkout(*(bspHandle.GetBioseqCore()), (**iter), firstGi, gi, out);
	    if((int)bspHandle.GetBioseqCore()->GetInst().GetLength() > k_GetSubseqThreshhold){
	      string dumpGnlUrl = getDumpgnlLink((*iter)->GetSeqid(), 1, k_DumpGnlUrl);
	      out<<dumpGnlUrl<<"<img border=0 height=16 width=16 src=\"/blast/images/D.gif\" alt=\"Download subject sequence spanning the HSP\"></a>";
	    }
	  }
	}
 
	out <<" ";
	if((*iter)->IsSetTitle()){
	  s_WrapOutputLine(out, (*iter)->GetTitle());     
	}
	out<<endl;
	isFirst = false;
      }
    }
  }
}

//Output sequence and mask sequences if any
const void CDisplaySeqalign::OutputSeq(string& sequence, const CSeq_id& id, int start, int len, int frame, bool colorMismatch, CNcbiOstream& out) const {
  int actualSize = sequence.size();
  assert(actualSize > start);
  list<CRange<int> > actualSeqloc;
  string actualSeq = sequence.substr(start, len);
  
  if(id.Which() != CSeq_id::e_not_set){ /*only do this for sequence but not for others like middle line, features*/
    //go through seqloc containing mask info
    for (list<alnSeqlocInfo*>::const_iterator iter = m_Alnloc.begin();  iter != m_Alnloc.end(); iter++){
      int from=(*iter)->alnRange.GetFrom();
      int to=(*iter)->alnRange.GetTo();
      int locFrame = (*iter)->seqloc->frame;
      if(id.Match((*iter)->seqloc->seqloc->GetInt().GetId()) && locFrame == frame){
	bool isFirstChar = true;
	CRange<int> eachSeqloc(0, 0);
	//go through each residule and mask it
	for (int i=max<int>(from, start); i<=min<int>(to, start+len); i++){
	  //store seqloc start for font tag below
	  if ((m_AlignOption & eHtml) && isFirstChar){         
	    isFirstChar = false;
	    eachSeqloc.Set(i, eachSeqloc.GetTo());
	  }
	  if (m_SeqLocChar==eX){
	    actualSeq[i-start]='X';
	  } else if (m_SeqLocChar==eN){
	    actualSeq[i-start]='n';
	  } else if (m_SeqLocChar==eLowerCase){
	    actualSeq[i-start]=tolower(actualSeq[i-start]);
	  }
	  //store seqloc start for font tag below
	  if ((m_AlignOption & eHtml) && i == min<int>(to, start+len)){ 
	    eachSeqloc.Set(eachSeqloc.GetFrom(), i);
	  }
	}
	if(!(eachSeqloc.GetFrom()==0&&eachSeqloc.GetTo()==0)){
	  actualSeqloc.push_back(eachSeqloc);
	}
      }
    }
  }

  if(actualSeqloc.empty()){//no need to add font tag
    if((m_AlignOption & eColorDifferentBases) && (m_AlignOption & eHtml) && colorMismatch){
      //color the mismatches. Only for rows without mask.  Otherwise it may confilicts with mask font tag.
      x_ColorDifferentBases(actualSeq, k_IdentityChar, out);
    } else {
      out<<actualSeq;
    }
  } else {//now deal with font tag for mask for html display    
    bool endTag = false;
    bool numFrontTag = 0;
    for (int i = 0; i < (int)actualSeq.size(); i ++){
      for (list<CRange<int> >::iterator iter=actualSeqloc.begin();  iter!=actualSeqloc.end(); iter++){
        int from = (*iter).GetFrom() - start;
        int to = (*iter).GetTo() - start;
	//start tag
        if(from == i){
          out<<"<font color=\""+color[m_SeqLocColor]+"\">";
          numFrontTag = 1;
        }
	//need to close tag at the end of mask or end of sequence
        if(to == i || i == (int)actualSeq.size() - 1 ){
          endTag = true;
        }
      }
      out<<actualSeq[i];
      if(endTag && numFrontTag == 1){
        out<<"</font>";
        endTag = false;
        numFrontTag = 0;
      }
    }
  }
}

int CDisplaySeqalign::getNumGaps() {
  int gap = 0;
  for (int row=0; row<m_AV->GetNumRows(); row++) {
    CRef<CAlnMap::CAlnChunkVec> chunk_vec = m_AV->GetAlnChunks(row, m_AV->GetSeqAlnRange(0));
    for (int i=0; i<chunk_vec->size(); i++) {
      CConstRef<CAlnMap::CAlnChunk> chunk = (*chunk_vec)[i];
      if (chunk->IsGap()) {
        gap += (chunk->GetAlnRange().GetTo() - chunk->GetAlnRange().GetFrom() + 1);
      }
    }
  }
  return gap;
}


const CRef<CBlast_def_line_set>  CDisplaySeqalign::GetBlastDefline (const CBioseq_Handle& handle) const {
  CRef<CBlast_def_line_set> bdls(new CBlast_def_line_set());
  if(handle.IsSetDescr()){
    const CSeq_descr& desc = handle.GetDescr();
    const list< CRef< CSeqdesc > >& descList = desc.Get();
    for (list<CRef< CSeqdesc > >::const_iterator iter = descList.begin(); iter != descList.end(); iter++){
      
      if((*iter)->IsUser()){
        const CUser_object& uobj = (*iter)->GetUser();
        const CObject_id& uobjid = uobj.GetType();
        if(uobjid.IsStr()){
   
          const string& label = uobjid.GetStr();
          if (label == kAsnDeflineObjLabel){
           const vector< CRef< CUser_field > >& usf = uobj.GetData();
           string buf;
          
	   if(usf.front()->GetData().IsOss()){ //only one user field
             typedef const CUser_field::TData::TOss TOss;
             const TOss& oss = usf.front()->GetData().GetOss();
             int size = 0;
             //determine the octet string length
             ITERATE (TOss, iter3, oss) {
	       size += (**iter3).size();
             }
             
             int i =0;
             char* temp = new char[size];
             //retrive the string
             ITERATE (TOss, iter3, oss) {
      
               for(vector< char >::iterator iter4 = (**iter3).begin(); iter4 !=(**iter3).end(); iter4++){
                 temp[i] = *iter4;
                 i++;
               }
             }            
	    
	     CConn_MemoryStream stream;
             stream.write(temp, i);
             auto_ptr<CObjectIStream> ois(CObjectIStream::Open(eSerial_AsnBinary, stream));
             *ois >> *bdls;
	     delete [] temp;
           }         
          }
        }
      }
    }
  }
  return bdls;
}

static string GetTaxNames(const CBioseq& cbsp, int taxid){
  string name;
  if(cbsp.IsSetDescr()){  
    const CSeq_descr& desc = cbsp.GetDescr();
    const list< CRef< CSeqdesc > >& descList = desc.Get();   
    for (list<CRef< CSeqdesc > >::const_iterator iter = descList.begin(); iter != descList.end(); iter++){
      if((*iter)->IsUser()){
        const CUser_object& uobj = (*iter)->GetUser();
        const CObject_id& uobjid = uobj.GetType();
        if(uobjid.IsStr()){   
          const string& label = uobjid.GetStr();
          if (label == kTaxDataObjLabel){
            const vector< CRef< CUser_field > >& usf = uobj.GetData();        
            for (vector< CRef< CUser_field > >::const_iterator iter2 = usf.begin(); iter2 != usf.end(); iter2 ++){
              const CObject_id& oid = (*iter2)->GetLabel();
              if (oid.GetId() == taxid){
                (**iter2).GetData().Which();
                name = (**iter2).GetData().GetStrs().front();
                break;
              }
            }
          }
        }
      }
    }
  }
  return name;
}

void CDisplaySeqalign::getFeatureInfo(list<alnFeatureInfo*>& feature, CScope& scope, CSeqFeatData::E_Choice choice, int row, string& sequence) const {
  //Only fetch features for seq that has a gi
  CRef<CSeq_id> id = GetSeqIdByType(m_AV->GetBioseqHandle(row).GetBioseqCore()->GetId(), CSeq_id::e_Gi);
  if(!(id.Empty())){
    const CBioseq_Handle& handle = scope.GetBioseqHandle(*id);
    //cds feature
    for  (CFeat_CI feat (handle, m_AV->GetSeqStart(row), m_AV->GetSeqStop(row), choice); feat;  ++feat) {
      
      const CSeq_loc& loc = feat->GetLocation();
      string featLable = NcbiEmptyString;
      string featId;
      string alternativeFeatStr = NcbiEmptyString;

      feature::GetLabel(feat->GetOriginalFeature(), &featLable, feature::eBoth, &(m_AV->GetScope()));
      featId = featLable.substr(0, k_FeatureIdLen); //default
	  
      int alnStop = m_AV->GetAlnStop();      
      if(loc.IsInt()){
	alnFeatureInfo* featInfo = new alnFeatureInfo;