showalign.cpp

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

/*
 * ===========================================================================
 * PRODUCTION $Log: showalign.cpp,v $
 * PRODUCTION Revision 1000.5  2004/06/01 19:41:04  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.36
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: showalign.cpp,v 1000.5 2004/06/01 19:41:04 gouriano Exp $
 * ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a "United States Government Work" under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's official duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 *
 * Author:  Jian Ye
 *
 * File Description:
 *   Sequence alignment display
 *
 */

#include <ncbi_pch.hpp>
#include <objtools/alnmgr/util/showalign.hpp>

#include <corelib/ncbiexpt.hpp>
#include <corelib/ncbiutil.hpp>
#include <corelib/ncbistre.hpp>
#include <corelib/ncbireg.hpp>

#include <util/range.hpp>
#include <util/md5.hpp>

#include <objects/general/Object_id.hpp>
#include <objects/general/User_object.hpp>
#include <objects/general/User_field.hpp>
#include <objects/general/Dbtag.hpp>

#include <serial/iterator.hpp>
#include <serial/objistr.hpp>
#include <serial/objostr.hpp>
#include <serial/serial.hpp>
#include <serial/objostrasnb.hpp> 
#include <serial/objistrasnb.hpp> 
#include <connect/ncbi_conn_stream.hpp>

#include <objmgr/object_manager.hpp>
#include <objmgr/scope.hpp>
#include <objmgr/feat_ci.hpp>
#include <objtools/data_loaders/genbank/gbloader.hpp>

#include <objmgr/util/sequence.hpp>
#include <objmgr/util/feature.hpp>

#include <objects/seqfeat/SeqFeatData.hpp>
#include <objects/seqfeat/Cdregion.hpp>

#include <objects/seq/Seq_descr.hpp>
#include <objects/seq/Seqdesc.hpp>
#include <objects/seq/Bioseq.hpp>

#include <objects/seqset/Seq_entry.hpp>

#include <objects/seqloc/Seq_id.hpp>
#include <objects/seqloc/Seq_interval.hpp>

#include <objects/seqalign/Seq_align_set.hpp>
#include <objects/seqalign/Score.hpp>
#include <objects/seqalign/Std_seg.hpp>
#include <objects/seqalign/Dense_diag.hpp>

#include <objtools/alnmgr/alnmix.hpp>
#include <objtools/alnmgr/alnvec.hpp>

#include <objects/blastdb/Blast_def_line.hpp>
#include <objects/blastdb/Blast_def_line_set.hpp>
#include <objects/blastdb/defline_extra.hpp>

#include <stdio.h>
#include <util/tables/raw_scoremat.h>


BEGIN_NCBI_SCOPE
BEGIN_SCOPE (objects)
USING_SCOPE (sequence);

static const char k_IdentityChar = '.';
static const int k_NumFrame = 6;
static const string k_FrameConversion[k_NumFrame] = {"+1", "+2", "+3", "-1", "-2", "-3"};
static const int k_GetSubseqThreshhold = 10000;
static const int k_ColorMismatchIdentity = 0;  /*threshhold to color mismatch.  98 means 98% */
static const string k_DumpGnlUrl = "/blast/dumpgnl.cgi";
static const int k_FeatureIdLen = 16;
static const int k_NumAsciiChar = 128;
const string color[]={"#000000", "#808080", "#FF0000"};
static const char k_PSymbol[CDisplaySeqalign::kPMatrixSize+1] = "ARNDCQEGHILKMFPSTWYVBZX";


/* url for linkout*/
static const string k_EntrezUrl = "<a href=\"http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=%s&list_uids=%d&dopt=%s\" %s>";

static const string k_TraceUrl = "<a href=\"http://www.ncbi.nlm.nih.gov/Traces/trace.cgi?cmd=retrieve&dopt=fasta&val=%s\">";

static const string k_LocusLinkUrl = "<a href=\"http://www.ncbi.nlm.nih.gov/LocusLink/list.cgi?Q=%d%s\"><img border=0 height=16 width=16 src=\"/blast/images/L.gif\" alt=\"LocusLink info\"></a>";

static const string k_UnigeneUrl = "<a href=\"http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=unigene&cmd=search&term=%d[Nucleotide+UID]\"><img border=0 height=16 width=16 src=\"/blast/images/U.gif\" alt=\"UniGene info\"></a>";

static const string k_StructureUrl = "<a href=\"http://www.ncbi.nlm.nih.gov/Structure/cblast/cblast.cgi?blast_RID=%s&blast_rep_gi=%d&hit=%d&blast_CD_RID=%s&blast_view=%s&hsp=0&taxname=%s&client=blast\"><img border=0 height=16 width=16 src=\"http://www.ncbi.nlm.nih.gov/Structure/cblast/str_link.gif\" alt=\"Related structures\"></a>";

static const string k_StructureOverviewUrl = "<a href=\"http://www.ncbi.nlm.nih.gov/Structure/cblast/cblast.cgi?blast_RID=%s&blast_rep_gi=%d&hit=%d&blast_CD_RID=%s&blast_view=%s&hsp=0&taxname=%s&client=blast\">Related Structures</a>";

static const string k_GeoUrl =  "<a href=\"http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=geo&term=%d[gi]\"><img border=0 height=16 width=16 src=\"/blast/images/G.gif\" alt=\"Geo\"></a>";


//Constructor
CDisplaySeqalign::CDisplaySeqalign(const CSeq_align_set& seqalign, list <SeqlocInfo*>& maskSeqloc, list <FeatureInfo*>& externalFeature, const int matrix[][kPMatrixSize], CScope& scope) : m_SeqalignSetRef(&seqalign), m_Seqloc(maskSeqloc), m_QueryFeature(externalFeature), m_Scope(scope) {
  m_AlignOption = 0;
  m_SeqLocChar = eX;
  m_SeqLocColor = eBlack;
  m_LineLen = 60;
  m_IdStartMargin = 2;
  m_StartSequenceMargin = 2;
  m_SeqStopMargin = 2;
  m_IsDbNa = true;
  m_IsQueryNa = true;
  m_IsDbGi = false;
  m_DbName = NcbiEmptyString;
  m_NumAlignToShow = 10000;
  m_AlignType = eNotSet;
  m_Rid = "0";
  m_CddRid = "0";
  m_EntrezTerm = NcbiEmptyString;
  m_QueryNumber = 0;
  m_BlastType = NcbiEmptyString;
  m_MidLineStyle = eBar;
  m_ConfigFile = NULL;
  m_Reg = NULL;

  SNCBIFullScoreMatrix blosumMatrix;
  NCBISM_Unpack(&NCBISM_Blosum62, &blosumMatrix);
 
  int** temp = new int*[k_NumAsciiChar];
  for(int i = 0; i<k_NumAsciiChar; ++i) {
    temp[i] = new int[k_NumAsciiChar];
  }
  for (int i=0; i<k_NumAsciiChar; i++){
    for (int j=0; j<k_NumAsciiChar; j++){
      temp[i][j] = -1000;
    }
  }
  for(int i = 0; i < kPMatrixSize; ++i){
    for(int j = 0; j < kPMatrixSize; ++j){
      if(matrix){
	temp[(size_t)k_PSymbol[i]][(size_t)k_PSymbol[j]] = matrix[i][j];
      } else {
	temp[(size_t)k_PSymbol[i]][(size_t)k_PSymbol[j]] = blosumMatrix.s[(size_t)k_PSymbol[i]][(size_t)k_PSymbol[j]];
      }
     
    }
  }
  for(int i = 0; i < kPMatrixSize; ++i) {
    temp[(size_t)k_PSymbol[i]]['*'] = temp['*'][(size_t)k_PSymbol[i]] = -4;
  }
  temp['*']['*'] = 1;
 
  m_Matrix = temp;
 
}

//Destructor
CDisplaySeqalign::~CDisplaySeqalign(){
  for(int i = 0; i<k_NumAsciiChar; ++i) {
    delete [] m_Matrix[i];
  }
  delete [] m_Matrix;
  if (m_ConfigFile) {
    delete m_ConfigFile;
  } 
  if (m_Reg) {
    delete m_Reg;
  }

}

static void AddSpace(CNcbiOstream& out, int number);

static string GetTaxNames(const CBioseq& cbsp, int taxid);

static string getNameInitials(string& name);
static string GetSeqForm(char* formName, bool dbIsNa, int queryNumber);
static const string GetSeqIdStringByFastaOrder(const CSeq_id& id, CScope& sp, bool with_version);
static int GetGiForSeqIdList (const list<CRef<CSeq_id> >& ids);

static string MakeURLSafe(char* src);
static void getAlnScores(const CSeq_align& aln, int& score, double& bits, double& evalue);
template<class container> static bool s_GetBlastScore(const container& scoreList,  int& score, double& bits, double& evalue);
static bool s_canDoMultiAlign(const CSeq_align_set& aln, CScope& scope);
static CRef<CSeq_id> GetSeqIdByType(const list<CRef<CSeq_id> >& ids, CSeq_id::E_Choice choice);
static int s_getFrame (int start, ENa_strand strand, const CSeq_id& id, CScope& sp);
static CRef<CSeq_align> CreateDensegFromDendiag(const CSeq_align& aln);
static void x_ColorDifferentBases(string& seq, char identityChar, CNcbiOstream& out);
static void s_WrapOutputLine(CNcbiOstream& out, const string& str);

static void s_WrapOutputLine(CNcbiOstream& out, const string& str){
  const int line_len = 60;
  const int front_space = 12;
  bool do_wrap = false;
  int length = (int) str.size();
  if (length > line_len) {
    for (int i = 0; i < length; i ++){
      if(i > 0 && i % line_len == 0){
	do_wrap = true;
      }   
      out << str[i];
      if(do_wrap && isspace(str[i])){
	out << endl;  
	AddSpace(out, front_space);
	do_wrap = false;
      }
    }
  } else {
    out << str;
  }
}


//To add color to bases other than identityChar
static void x_ColorDifferentBases(string& seq, char identityChar, CNcbiOstream& out){
  string base_color = "#FF0000";
  bool tagOpened = false;
  for(int i = 0; i < (int)seq.size(); i ++){
    if(seq[i] != identityChar){
      if(!tagOpened){
	out << "<font color=\""+base_color+"\"><b>";
	tagOpened =  true;
      }
     
    } else {
      if(tagOpened){
	out << "</b></font>";
	tagOpened = false;
      }
    }
    out << seq[i];
    if(tagOpened && i == (int)seq.size() - 1){
      out << "</b></font>";
      tagOpened = false;
    }
  } 
}


template<class container> static bool s_GetBlastScore(const container&  scoreList,  int& score, double& bits, double& evalue){
  bool hasScore = false;
  ITERATE (typename container, iter, scoreList) {
    const CObject_id& id=(*iter)->GetId();
    if (id.IsStr()) {
      hasScore = true;
      if (id.GetStr()=="score"){
	score = (*iter)->GetValue().GetInt();
      
      } else if (id.GetStr()=="bit_score"){
	bits = (*iter)->GetValue().GetReal();
	
      } else if (id.GetStr()=="e_value" || id.GetStr()=="sum_e") {
	evalue = (*iter)->GetValue().GetReal();
      } 
    }
  }
  return hasScore;
}
static CRef<CSeq_align> CreateDensegFromDendiag(const CSeq_align& aln) {
  CRef<CSeq_align> sa(new CSeq_align);
  if ( !aln.GetSegs().IsDendiag()) {
    NCBI_THROW(CException, eUnknown, "Input Seq-align should be Dendiag!");
  }
  
  if(aln.IsSetType()){
    sa->SetType(aln.GetType());
  }
  if(aln.IsSetDim()){
    sa->SetDim(aln.GetDim());
  }
  if(aln.IsSetScore()){
    sa->SetScore() = aln.GetScore();
  }
  if(aln.IsSetBounds()){
    sa->SetBounds() = aln.GetBounds();
  }
  
  CDense_seg& ds = sa->SetSegs().SetDenseg();

  int counter = 0;
  ds.SetNumseg() = 0;
  ITERATE (CSeq_align::C_Segs::TDendiag, iter, aln.GetSegs().GetDendiag()){
     
    if(counter == 0){//assume all dendiag segments have same dim and ids
      if((*iter)->IsSetDim()){
	ds.SetDim((*iter)->GetDim());
      }
      if((*iter)->IsSetIds()){
	ds.SetIds() = (*iter)->GetIds();
      }
    }
    ds.SetNumseg() ++;
    if((*iter)->IsSetStarts()){
      ITERATE(CDense_diag::TStarts, iterStarts, (*iter)->GetStarts()){
	ds.SetStarts().push_back(*iterStarts);
      }
    }
    if((*iter)->IsSetLen()){
      ds.SetLens().push_back((*iter)->GetLen());
    }
    if((*iter)->IsSetStrands()){
      ITERATE(CDense_diag::TStrands, iterStrands, (*iter)->GetStrands()){
	ds.SetStrands().push_back(*iterStrands);
      }
    }
    if((*iter)->IsSetScores()){
      ITERATE(CDense_diag::TScores, iterScores, (*iter)->GetScores()){
	ds.SetScores().push_back(*iterScores); //this might not have right meaning
      }
    }
    counter ++;
  }
  
  return sa;
}
static int s_GetStdsegMasterFrame(const CStd_seg& ss, CScope& scope){
  const CRef<CSeq_loc> slc = ss.GetLoc().front();
  ENa_strand strand = GetStrand(*slc);
  int frame = s_getFrame(strand ==  eNa_strand_plus ? GetStart(*slc) : GetStop(*slc), strand ==  eNa_strand_plus? eNa_strand_plus : eNa_strand_minus, *(ss.GetIds().front()), scope);
  
  return frame;
}

/*Note that start is zero bases.  It returns frame +/1(1-3).  ) indicates erro*/
static int s_getFrame (int start, ENa_strand strand, const CSeq_id& id, CScope& sp) {
  int frame = 0;
  if (strand == eNa_strand_plus) {
    frame = (start % 3) + 1;
  } else if (strand == eNa_strand_minus) {
    frame = -((sp.GetBioseqHandle(id).GetBioseqCore()->GetInst().GetLength() - start - 1) % 3 + 1);
   
  }
  return frame;
}
static bool s_canDoMultiAlign(const CSeq_align_set& aln, CScope& scope){