showalign.cpp

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

  bool multiAlign = true;
  int firstFrame = 0;
  bool isFirstSeqalign = true;
  /*Make sure that in case of stdseg, the master must be in same frame for all pairwise for multialign view */
  for (CTypeConstIterator<CSeq_align> seqalign = ConstBegin(aln); seqalign; ++ seqalign){
    if (seqalign->GetSegs().Which() == CSeq_align::C_Segs::e_Denseg){
      break;
    } else if (seqalign->GetSegs().Which() == CSeq_align::C_Segs::e_Std){
      CTypeConstIterator<CStd_seg> ss = ConstBegin(*seqalign);     
      int curFrame = s_GetStdsegMasterFrame(*ss, scope);
      if (isFirstSeqalign){
	firstFrame = curFrame;
	isFirstSeqalign = false;
      } else {
	if (firstFrame != curFrame){
	  multiAlign = false;
	  break;
	}
       
      }
    }
   
  }
  return multiAlign;
}


static void getAlnScores(const CSeq_align& aln, int& score, double& bits, double& evalue){
  bool hasScore = false;
  //look for scores at seqalign level first
  hasScore = s_GetBlastScore(aln.GetScore(),  score, bits, evalue);

  //look at the seg level
  if(!hasScore){
    const CSeq_align::TSegs& seg = aln.GetSegs();
    if(seg.Which() == CSeq_align::C_Segs::e_Std){
      s_GetBlastScore(seg.GetStd().front()->GetScores(),  score, bits, evalue);
    } else if (seg.Which() == CSeq_align::C_Segs::e_Dendiag){
      s_GetBlastScore(seg.GetDendiag().front()->GetScores(),  score, bits, evalue);
    }  else if (seg.Which() == CSeq_align::C_Segs::e_Denseg){
       s_GetBlastScore(seg.GetDenseg().GetScores(),  score, bits, evalue);
    }
  }	
}

string CDisplaySeqalign::getUrl(const list<CRef<CSeq_id> >& ids, int row) const{
  string urlLink = NcbiEmptyString;
 
  char dopt[32], db[32];
  int gi = GetGiForSeqIdList(ids);
  string toolUrl= m_Reg->Get(m_BlastType, "TOOL_URL");
  if(toolUrl == NcbiEmptyString || (gi > 0 && toolUrl.find("dumpgnl.cgi") != string::npos)){ //use entrez or dbtag specified
    if(m_IsDbNa) {
      strcpy(dopt, "GenBank");
      strcpy(db, "Nucleotide");
    } else {
      strcpy(dopt, "GenPept");
      strcpy(db, "Protein");
    }    
 
    char urlBuf[1024];
    if (gi > 0) {
      sprintf(urlBuf, k_EntrezUrl.c_str(), db, gi, dopt, (m_AlignOption & eNewTargetWindow) ? "TARGET=\"EntrezView\"" : "");
      urlLink = urlBuf;
    } else {//seqid general, dbtag specified
      const CRef<CSeq_id> wid = FindBestChoice(ids, CSeq_id::WorstRank);
      if(wid->Which() == CSeq_id::e_General){
        const CDbtag& dtg = wid->GetGeneral();
        const string& dbName = dtg.GetDb();
        if(NStr::CompareNocase(dbName, "TI") == 0){
          sprintf(urlBuf, k_TraceUrl.c_str(), wid->GetSeqIdString().c_str());
          urlLink = urlBuf;
        } else { //future use

        }
      }
    }
  } else { //need to use url in configuration file
    string altUrl = NcbiEmptyString;
      urlLink = getDumpgnlLink(ids, row, altUrl);
  }
  return urlLink;
}

static string getNameInitials(string& name){
  vector<string> arr;
  string initials;
  NStr::Tokenize(name, " ", arr);

  for(vector<string>::iterator iter = arr.begin(); iter != arr.end(); iter ++){
    if (*iter != NcbiEmptyString){
      initials += (*iter)[0];
    }
  }
  return initials;
}

static void AddSpace(CNcbiOstream& out, int number){
  for(int i=0; i<number; i++){
    out<<" ";
  }
}


void CDisplaySeqalign::AddLinkout(const CBioseq& cbsp, const CBlast_def_line& bdl, int firstGi, int gi, CNcbiOstream& out) const{

  char molType[8]={""};
  if(cbsp.IsAa()){
    sprintf(molType, "[pgi]");
  }
  else {
    sprintf(molType, "[ngi]");
  }
 
  if (bdl.IsSetLinks()){
    for (list< int >::const_iterator iter = bdl.GetLinks().begin(); iter != bdl.GetLinks().end(); iter ++){
      char buf[1024];
      if((*iter) & eLocuslink){
        sprintf(buf, k_LocusLinkUrl.c_str(), gi, molType);
        out << buf;
      }
      if ((*iter) & eUnigene) {
	sprintf(buf, k_UnigeneUrl.c_str(),  gi);
	out << buf;
      }
      if ((*iter) & eStructure){
        sprintf(buf, k_StructureUrl.c_str(), m_Rid.c_str(), firstGi, gi, m_CddRid.c_str(), "onepair", (m_EntrezTerm == NcbiEmptyString) ? "none":((char*) m_EntrezTerm.c_str()));
	out << buf;
      }
      if ((*iter) & eGeo){
         sprintf(buf, k_GeoUrl.c_str(), gi);
         out << buf;
      }
    }
  }
}

//return the get sequence table for html display
static string GetSeqForm(char* formName, bool dbIsNa, int queryNumber){
  char buf[2048] = {""};
  if(formName){
    sprintf(buf, "<table border=\"0\"><tr><td><FORM  method=\"post\" action=\"http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?SUBMIT=y\" name=\"%s\"><input type=button value=\"Get selected sequences\" onClick=\"finalSubmit(%d, 'getSeqAlignment%d', 'getSeqGi', '%s')\"><input type=\"hidden\" name=\"db\" value=\"\"><input type=\"hidden\" name=\"term\" value=\"\"><input type=\"hidden\" name=\"doptcmdl\" value=\"docsum\"><input type=\"hidden\" name=\"cmd\" value=\"search\"></form></td><td><FORM><input type=\"button\" value=\"Select all\" onClick=\"handleCheckAll('select', 'getSeqAlignment%d', 'getSeqGi')\"></form></td><td><FORM><input type=\"button\" value=\"Deselect all\" onClick=\"handleCheckAll('deselect', 'getSeqAlignment%d', 'getSeqGi')\"></form></td></tr></table>", formName, dbIsNa?1:0, queryNumber, formName, queryNumber, queryNumber);
   
  }
  return buf;
}

//Return the seqid in "fasta" orders
static const string GetSeqIdStringByFastaOrder(const CSeq_id& id, CScope& sp, bool with_version){
  string idString = NcbiEmptyString;
  static const int total_seqid_types=19;
  static int fasta_order[total_seqid_types];
  //fasta order.  See seqidwrite() in C library
  fasta_order[CSeq_id::e_not_set]=33;
  fasta_order[CSeq_id::e_Local]=20;
  fasta_order[CSeq_id::e_Gibbsq]=15;
  fasta_order[CSeq_id::e_Gibbmt]=16;
  fasta_order[CSeq_id::e_Giim]=30;
  fasta_order[CSeq_id::e_Genbank]=10;   
  fasta_order[CSeq_id::e_Embl]=10;
  fasta_order[CSeq_id::e_Pir]=10;
  fasta_order[CSeq_id::e_Swissprot]=10;
  fasta_order[CSeq_id::e_Patent]=15;
  fasta_order[CSeq_id::e_Other]=12;
  fasta_order[CSeq_id::e_General]=13;
  fasta_order[CSeq_id::e_Gi]=255;
  fasta_order[CSeq_id::e_Ddbj]=10;
  fasta_order[CSeq_id::e_Prf]=10;
  fasta_order[CSeq_id::e_Pdb]=12;
  fasta_order[CSeq_id::e_Tpg]=10;
  fasta_order[CSeq_id::e_Tpe]=10;
  fasta_order[CSeq_id::e_Tpd]=10;
 
  CRef<CSeq_id> idRef;
  const list<CRef<CSeq_id> >& ids=sp.GetBioseqHandle(id).GetBioseqCore()->GetId();
 
  for (CBioseq::TId::const_iterator iter = ids.begin(); iter != ids.end(); iter ++){
    if(iter == ids.begin()){
      idRef = *iter;
    }
    else if(fasta_order[(*iter)->Which()]<fasta_order[idRef->Which()]){
      idRef = *iter;
    }
    
  }
  if(!(ids.empty())){
    idString = idRef->GetSeqIdString(with_version);
  }
  return idString;
}

static CRef<CSeq_id> GetSeqIdByType(const list<CRef<CSeq_id> >& ids, CSeq_id::E_Choice choice) {
  CRef<CSeq_id> cid;

  for (CBioseq::TId::const_iterator iter = ids.begin(); iter != ids.end(); iter ++){
    if ((*iter)->Which() == choice){
      cid = *iter;
      break;
    }
  }
 
  return cid;
}

static int GetGiForSeqIdList (const list<CRef<CSeq_id> >& ids){
  int gi = 0;
  CRef<CSeq_id> id = GetSeqIdByType(ids, CSeq_id::e_Gi);
  if (!(id.Empty())){
    return id->GetGi();
  }
  return gi;
}

//To display the seqalign represented by internal alnvec
void CDisplaySeqalign::DisplayAlnvec(CNcbiOstream& out){ 

  int maxIdLen=0;
  int maxStartLen=0;
  int startLen=0;
  int actualLineLen=0;
  int aln_stop=m_AV->GetAlnStop();
  const int rowNum=m_AV->GetNumRows();
 
  if(m_AlignOption & eMasterAnchored){
    m_AV->SetAnchor(0);
  }
  m_AV->SetGapChar('-');
  m_AV->SetEndChar(' ');
  vector<string> sequence(rowNum);

  CAlnMap::TSeqPosList* seqStarts = new CAlnMap::TSeqPosList[rowNum];
  CAlnMap::TSeqPosList* seqStops = new CAlnMap::TSeqPosList[rowNum];
  CAlnMap::TSeqPosList* insertStart = new CAlnMap::TSeqPosList[rowNum];
  CAlnMap::TSeqPosList* insertAlnStart = new CAlnMap::TSeqPosList[rowNum];
  CAlnMap::TSeqPosList* insertLength = new CAlnMap::TSeqPosList[rowNum];

  string* seqidArray=new string[rowNum];
  string middleLine;
  list<alnFeatureInfo*>* bioseqFeature= new list<alnFeatureInfo*>[rowNum];
  CAlnMap::TSignedRange* rowRng = new CAlnMap::TSignedRange[rowNum];
  int* frame = new int[rowNum];

  //conver to aln coordinates for mask seqloc
  list<alnSeqlocInfo*> alnLocList;
  for (list<SeqlocInfo*>::iterator iter=m_Seqloc.begin();  iter!=m_Seqloc.end(); iter++){
    alnSeqlocInfo* alnloc = new alnSeqlocInfo;    
   
    for (int i=0; i<rowNum; i++){
     
      if((*iter)->seqloc->GetInt().GetId().Match(m_AV->GetSeqId(i))){
	int actualAlnStart = 0, actualAlnStop = 0;
	if(m_AV->IsPositiveStrand(i)){
	  actualAlnStart = m_AV->GetAlnPosFromSeqPos(i, (*iter)->seqloc->GetInt().GetFrom());
	  actualAlnStop = m_AV->GetAlnPosFromSeqPos(i, (*iter)->seqloc->GetInt().GetTo());
	} else {
	  actualAlnStart = m_AV->GetAlnPosFromSeqPos(i, (*iter)->seqloc->GetInt().GetTo());
	  actualAlnStop = m_AV->GetAlnPosFromSeqPos(i, (*iter)->seqloc->GetInt().GetFrom());
	}
        alnloc->alnRange.Set(actualAlnStart, actualAlnStop);      
	break;
      }
    }
    alnloc->seqloc = *iter;   
    alnLocList.push_back(alnloc);
  }
  m_Alnloc = alnLocList;
  //Add external query feature info such as phi blast pattern
  for (list<FeatureInfo*>::iterator iter=m_QueryFeature.begin();  iter!=m_QueryFeature.end(); iter++){
    for(int i = 0; i < rowNum; i++){
      if((*iter)->seqloc->GetInt().GetId().Match(m_AV->GetSeqId(i))){
	int actualSeqStart = 0, actualSeqStop = 0;
	if(m_AV->IsPositiveStrand(i)){
	  if((*iter)->seqloc->GetInt().GetFrom() < m_AV->GetSeqStart(i)){
	    actualSeqStart = m_AV->GetSeqStart(i);
	  } else {
	    actualSeqStart = (*iter)->seqloc->GetInt().GetFrom();
	  }

	  if((*iter)->seqloc->GetInt().GetTo() > m_AV->GetSeqStop(i)){
	    actualSeqStop = m_AV->GetSeqStop(i);
	  } else {
	    actualSeqStop = (*iter)->seqloc->GetInt().GetTo();
	  }
	} else {
	  if((*iter)->seqloc->GetInt().GetFrom() < m_AV->GetSeqStart(i)){
	    actualSeqStart = (*iter)->seqloc->GetInt().GetFrom();
	  } else {
	    actualSeqStart = m_AV->GetSeqStart(i);
	  }
	  
	  if((*iter)->seqloc->GetInt().GetTo() > m_AV->GetSeqStop(i)){
	    actualSeqStop = (*iter)->seqloc->GetInt().GetTo();
	  } else {
	    actualSeqStop = m_AV->GetSeqStop(i);
	  }
	}
	int alnFrom = m_AV->GetAlnPosFromSeqPos(i, actualSeqStart);
	int alnTo = m_AV->GetAlnPosFromSeqPos(i, actualSeqStop);
	
	alnFeatureInfo* featInfo = new alnFeatureInfo;
	string tempFeat = NcbiEmptyString;
	setFeatureInfo(featInfo, *((*iter)->seqloc), alnFrom, alnTo, aln_stop, (*iter)->featureChar, (*iter)->featureId, tempFeat);    
	bioseqFeature[i].push_back(featInfo);
      }
    }
  }

  //prepare data for each row
  for (int row=0; row<rowNum; row++) {
    rowRng[row] = m_AV->GetSeqAlnRange(row);
    frame[row] = (m_AV->GetWidth(row) == 3 ? s_getFrame(m_AV->IsPositiveStrand(row) ? m_AV->GetSeqStart(row) : m_AV->GetSeqStop(row), m_AV->IsPositiveStrand(row) ? eNa_strand_plus : eNa_strand_minus, m_AV->GetSeqId(row), m_Scope) : 0);
 
       
    //make sequence
    m_AV->GetWholeAlnSeqString(row, sequence[row],  &insertAlnStart[row], &insertStart[row], &insertLength[row], m_LineLen, &seqStarts[row], &seqStops[row]);
    
    //make feature.  Only for pairwise, non-query-anchored and untranslated
    if(!(m_AlignOption & eMasterAnchored) && !(m_AlignOption & eMultiAlign) && m_AV->GetWidth(row) != 3){
      if(m_AlignOption & eShowCdsFeature){
	getFeatureInfo(bioseqFeature[row], *m_featScope, CSeqFeatData::e_Cdregion, row, sequence[row]);
      }
      if(m_AlignOption & eShowGeneFeature){
	getFeatureInfo(bioseqFeature[row], *m_featScope, CSeqFeatData::e_Gene, row, sequence[row]);
      }
    }
 
    //make id
    if(m_AlignOption & eShowBlastStyleId) {
       if(row==0){//query
         seqidArray[row]="Query";
       } else {//hits
         if (!(m_AlignOption&eMultiAlign)){
         //hits for pairwise 
           seqidArray[row]="Sbjct";
         } else {
           int gi = GetGiForSeqIdList(m_AV->GetBioseqHandle(row).GetBioseqCore()->GetId());
           if (m_AlignOption&eShowGi && gi > 0){
             seqidArray[row]=NStr::IntToString(gi);
           } else {
             const CRef<CSeq_id> wid = FindBestChoice(m_AV->GetBioseqHandle(row).GetBioseqCore()->GetId(), CSeq_id::WorstRank);
             seqidArray[row]=wid->GetSeqIdString();
           }           
         }
       }

    } else {
      int gi = GetGiForSeqIdList(m_AV->GetBioseqHandle(row).GetBioseqCore()->GetId());
      if (m_AlignOption&eShowGi && gi > 0){
        seqidArray[row]=NStr::IntToString(gi);
      } else {
       
        const CRef<CSeq_id> wid = FindBestChoice(m_AV->GetBioseqHandle(row).GetBioseqCore()->GetId(), CSeq_id::WorstRank);
        seqidArray[row]=wid->GetSeqIdString();
      }      
    }

    //max id length
    maxIdLen=max<int>(seqidArray[row].size(), maxIdLen);
    //max start length
    int maxCood=max<int>(m_AV->GetSeqStart(row), m_AV->GetSeqStop(row));
    maxStartLen = max<int>(NStr::IntToString(maxCood).size(), maxStartLen);
  }
  //adjust max id length for feature id 
  for(int i = 0; i < rowNum; i ++){
    for (list<alnFeatureInfo*>::iterator iter=bioseqFeature[i].begin();  iter != bioseqFeature[i].end(); iter++){
      maxIdLen=max<int>((*iter)->feature->featureId.size(), maxIdLen );
    }