C++ SiteContainer::getNumberOfSites方法代碼示例

void SiteContainerTools::removeGapOrUnresolvedOnlySites(SiteContainer& sites)
{
  size_t n = sites.getNumberOfSites();
  size_t i = n;
  while (i > 1)
  {
    ApplicationTools::displayGauge(n - i + 1, n);
    const Site* site = &sites.getSite(i - 1);
    if (SiteTools::isGapOnly(*site))
    {
      size_t end = i;
      while (SiteTools::isGapOrUnresolvedOnly(*site) && i > 1)
      {
        --i;
        site = &sites.getSite(i - 1);
      }
      sites.deleteSites(i, end - i);
    }
    else
    {
      --i;
    }
  }
  ApplicationTools::displayGauge(n, n);
  const Site* site = &sites.getSite(0);
  if (SiteTools::isGapOrUnresolvedOnly(*site))
    sites.deleteSite(0);
}

void RecursiveLikelihoodTree::initLikelihoods(const SiteContainer& sites, const SubstitutionProcess& process)
throw (Exception)
{
  if (sites.getNumberOfSequences() == 1)
    throw Exception("RecursiveLikelihoodTree::initLikelihoods. Only 1 sequence in data set.");
  if (sites.getNumberOfSequences() == 0)
    throw Exception("RecursiveLikelihoodTree::initLikelihoods. No sequence in data set.");
  if (!process.isCompatibleWith(sites))
    throw Exception("RecursiveLikelihoodTree::initLikelihoods. Data and model are not compatible.");
  alphabet_ = sites.getAlphabet();
  nbStates_ = process.getNumberOfStates();
  nbSites_  = sites.getNumberOfSites();
  unique_ptr<SitePatterns> patterns;

  if (usePatterns_)
  {
    patterns.reset(initLikelihoodsWithPatterns_(process.getTree().getRootNode(), sites, process));
    shrunkData_.reset(patterns->getSites());
    rootWeights_      = patterns->getWeights();
    rootPatternLinks_ = patterns->getIndices();
    nbDistinctSites_  = shrunkData_->getNumberOfSites();

    setPatterns(patternLinks_);
  }
  else
  {
    patterns.reset(new SitePatterns(&sites));
    shrunkData_.reset(patterns->getSites());
    rootWeights_      = patterns->getWeights();
    rootPatternLinks_ = patterns->getIndices();
    nbDistinctSites_  = shrunkData_->getNumberOfSites();
    initLikelihoodsWithoutPatterns_(process.getTree().getRootNode(), *shrunkData_, process);
  }
}

void SiteContainerTools::getSequencePositions(const SiteContainer& sites, Matrix<size_t>& positions)
{
  positions.resize(sites.getNumberOfSequences(), sites.getNumberOfSites());
  int gap = sites.getAlphabet()->getGapCharacterCode();
  for (size_t i = 0; i < sites.getNumberOfSequences(); ++i) {
    const Sequence& seq = sites.getSequence(i);
    unsigned int pos = 0;
    for (size_t j = 0; j < sites.getNumberOfSites(); ++j) {
      if (seq[j] != gap) {
        ++pos;
        positions(i, j) = pos;
      } else {
        positions(i, j) = 0;
      }
    }
  }
}

void SiteContainerTools::removeGapSites(SiteContainer& sites, double maxFreqGaps)
{
  for (size_t i = sites.getNumberOfSites(); i > 0; i--) {
    map<int, double> freq;
    SiteTools::getFrequencies(sites.getSite(i - 1), freq);
    if (freq[-1] > maxFreqGaps)
      sites.deleteSite(i - 1);
  }
}

int main() {
  //ProteicAlphabet* alpha = new ProteicAlphabet;
  RNA* alpha = new RNA();
  SiteContainer* sites = new VectorSiteContainer(alpha);
  BasicSequence seq1("seq1", "----AUGCCG---GCGU----UUU----G--G-CCGACGUGUUUU--", alpha);
  BasicSequence seq2("seq2", "---GAAGGCG---G-GU----UUU----GC-GACCGACG--UUUU--", alpha);
  sites->addSequence(seq1, false);
  sites->addSequence(seq2, false);

  cout << sites->getNumberOfSites() << endl;
  cout << sites->toString("seq1") << endl;
  cout << sites->toString("seq2") << endl;
  SiteContainerTools::removeGapOnlySites(*sites);
  cout << endl;
  
  cout << sites->getNumberOfSites() << endl;
  cout << sites->toString("seq1") << endl;
  cout << sites->toString("seq2") << endl;

  return (sites->getNumberOfSites() == 30 ? 0 : 1);
}

void SiteContainerTools::changeUnresolvedCharactersToGaps(SiteContainer& sites)
{
  // NB: use iterators for a better algorithm?
  int gapCode = sites.getAlphabet()->getGapCharacterCode();
  for (unsigned int i = 0; i < sites.getNumberOfSites(); i++)
  {
    for (unsigned int j = 0; j < sites.getNumberOfSequences(); j++)
    {
      int* element = &sites(j, i);
      if (sites.getAlphabet()->isUnresolved(*element))
        *element = gapCode;
    }
  }
}

void SiteContainerTools::changeGapsToUnknownCharacters(SiteContainer& sites)
{
  // NB: use iterators for a better algorithm?
  int unknownCode = sites.getAlphabet()->getUnknownCharacterCode();
  for (unsigned int i = 0; i < sites.getNumberOfSites(); i++)
  {
    for (unsigned int j = 0; j < sites.getNumberOfSequences(); j++)
    {
      int* element = &sites(j, i);
      if (sites.getAlphabet()->isGap(*element))
        *element = unknownCode;
    }
  }
}

VectorSiteContainer::VectorSiteContainer(const SiteContainer& sc) :
  AbstractSequenceContainer(sc),
  sites_(0),
  names_(sc.getSequencesNames()),
  comments_(sc.getNumberOfSequences()),
  sequences_(sc.getNumberOfSequences())
{
  // Now try to add each site:
  for (size_t i = 0; i < sc.getNumberOfSites(); i++)
  {
    addSite(sc.getSite(i), false); // We assume that positions are correct.
  }
  // Seq comments:
  for (size_t i = 0; i < sc.getNumberOfSequences(); i++)
  {
    comments_[i] = new Comments(sc.getComments(i));
  }
}

void SiteContainerTools::merge(SiteContainer& seqCont1, const SiteContainer& seqCont2, bool leavePositionAsIs)
throw (AlphabetMismatchException, Exception)
{
  if (seqCont1.getAlphabet()->getAlphabetType() != seqCont2.getAlphabet()->getAlphabetType())
    throw AlphabetMismatchException("SiteContainerTools::merge.", seqCont1.getAlphabet(), seqCont2.getAlphabet());


  vector<string> seqNames1 = seqCont1.getSequencesNames();
  vector<string> seqNames2 = seqCont2.getSequencesNames();
  const SiteContainer* seqCont2bis = 0;
  bool del = false;
  if (seqNames1 == seqNames2)
  {
    seqCont2bis = &seqCont2;
  }
  else
  {
    // We shall reorder sequences first:
    SiteContainer* seqCont2ter = new VectorSiteContainer(seqCont2.getAlphabet());
    SequenceContainerTools::getSelectedSequences(seqCont2, seqNames1, *seqCont2ter);
    seqCont2bis = seqCont2ter;
    del = true;
  }

  if (leavePositionAsIs)
  {
    for (size_t i = 0; i < seqCont2bis->getNumberOfSites(); i++)
    {
      seqCont1.addSite(seqCont2bis->getSite(i), false);
    }
  }
  else
  {
    int offset = static_cast<int>(seqCont1.getNumberOfSites());
    for (size_t i = 0; i < seqCont2bis->getNumberOfSites(); i++)
    {
      seqCont1.addSite(seqCont2bis->getSite(i), offset + seqCont2bis->getSite(i).getPosition(), false);
    }
  }

  if (del)
    delete seqCont2bis;
}

std::map<size_t, size_t> SiteContainerTools::translateSequence(const SiteContainer& sequences, size_t i1, size_t i2)
{
  const Sequence* seq1 = &sequences.getSequence(i1);
  const Sequence* seq2 = &sequences.getSequence(i2);
  map<size_t, size_t> tln;
  size_t count1 = 0; // Sequence 1 counter
  size_t count2 = 0; // Sequence 2 counter
  int state1;
  int state2;
  for (size_t i = 0; i <  sequences.getNumberOfSites(); i++)
  {
    state1 = (*seq1)[i];
    if (state1 != -1)
      count1++;
    state2 = (*seq2)[i];
    if (state2 != -1)
      count2++;
    if (state1 != -1)
    {
      tln[count1] = (state2 == -1 ? 0 : count2);
    }
  }
  return tln;
}

VectorSiteContainer* SequenceApplicationTools::getSitesToAnalyse(
  const SiteContainer& allSites,
  map<string, string>& params,
  string suffix,
  bool suffixIsOptional,
  bool gapAsUnknown,
  bool verbose,
  int warn)
{
  // Fully resolved sites, i.e. without jokers and gaps:
  SiteContainer* sitesToAnalyse;
  VectorSiteContainer* sitesToAnalyse2;

  string option = ApplicationTools::getStringParameter("input.sequence.sites_to_use", params, "complete", suffix, suffixIsOptional, warn);
  if (verbose)
    ApplicationTools::displayResult("Sites to use", option);
  if (option == "all")
  {
    sitesToAnalyse = new VectorSiteContainer(allSites);
    string maxGapOption = ApplicationTools::getStringParameter("input.sequence.max_gap_allowed", params, "100%", suffix, suffixIsOptional, warn);

    if (maxGapOption[maxGapOption.size() - 1] == '%')
    {
      double gapFreq = TextTools::toDouble(maxGapOption.substr(0, maxGapOption.size() - 1)) / 100.;
      if (gapFreq < 1)
      {
        if (verbose)
          ApplicationTools::displayTask("Remove sites with gaps", true);
        for (size_t i = sitesToAnalyse->getNumberOfSites(); i > 0; --i)
        {
          if (verbose)
            ApplicationTools::displayGauge(sitesToAnalyse->getNumberOfSites() - i, sitesToAnalyse->getNumberOfSites() - 1, '=');
          map<int, double> freq;
          SiteTools::getFrequencies(sitesToAnalyse->getSite(i - 1), freq);
          if (freq[-1] > gapFreq)
            sitesToAnalyse->deleteSite(i - 1);
        }
        if (verbose)
          ApplicationTools::displayTaskDone();
      }
    }
    else
    {
      size_t gapNum = TextTools::to<size_t>(maxGapOption);
      if (gapNum < sitesToAnalyse->getNumberOfSequences())
      {
        if (verbose)
          ApplicationTools::displayTask("Remove sites with gaps", true);
        for (size_t i = sitesToAnalyse->getNumberOfSites(); i > 0; i--)
        {
          if (verbose)
            ApplicationTools::displayGauge(sitesToAnalyse->getNumberOfSites() - i, sitesToAnalyse->getNumberOfSites() - 1, '=');
          map<int, size_t> counts;
          SiteTools::getCounts(sitesToAnalyse->getSite(i - 1), counts);
          if (counts[-1] > gapNum)
            sitesToAnalyse->deleteSite(i - 1);
        }
        if (verbose)
          ApplicationTools::displayTaskDone();
      }
    }

    string maxUnresolvedOption = ApplicationTools::getStringParameter("input.sequence.max_unresolved_allowed", params, "100%", suffix, suffixIsOptional, warn);

    int sAlph = static_cast<int>(sitesToAnalyse->getAlphabet()->getSize());

    if (maxUnresolvedOption[maxUnresolvedOption.size() - 1] == '%')
    {
      double unresolvedFreq = TextTools::toDouble(maxUnresolvedOption.substr(0, maxUnresolvedOption.size() - 1)) / 100.;
      if (unresolvedFreq < 1)
      {
        if (verbose)
          ApplicationTools::displayTask("Remove unresolved sites", true);
        for (size_t i = sitesToAnalyse->getNumberOfSites(); i > 0; --i)
        {
          if (verbose)
            ApplicationTools::displayGauge(sitesToAnalyse->getNumberOfSites() - i, sitesToAnalyse->getNumberOfSites() - 1, '=');
          map<int, double> freq;
          SiteTools::getFrequencies(sitesToAnalyse->getSite(i - 1), freq);
          double x = 0;
          for (int l = 0; l < sAlph; ++l)
          {
            x += freq[l];
          }
          if (1 - x > unresolvedFreq)
            sitesToAnalyse->deleteSite(i - 1);
        }
        if (verbose)
          ApplicationTools::displayTaskDone();
      }
    }
    else
    {
      size_t nbSeq = sitesToAnalyse->getNumberOfSequences();
      size_t unresolvedNum = TextTools::to<size_t>(maxUnresolvedOption);
      if (unresolvedNum < nbSeq)
      {
        if (verbose)
          ApplicationTools::displayTask("Remove sites with gaps", true);
        for (size_t i = sitesToAnalyse->getNumberOfSites(); i > 0; i--)
//.........這裏部分代碼省略.........

//.........這裏部分代碼省略.........
          sc->addSequence(*seq);
        }
        delete sequences;
        sequences = sc;
      } else {
        VectorSiteContainer* sc = new VectorSiteContainer(sequences->getAlphabet());
        for (size_t i = 0; i < sequences->getNumberOfSequences(); ++i)
        {
          unique_ptr<Sequence> seq(sequences->getSequence(i).clone());
          SequenceTools::replaceStopsWithGaps(*seq, *gCode);
          sc->addSequence(*seq);
        }
        delete sequences;
        sequences = sc;
      }
    }

    // +--------------+
    // | Remove stops |
    // +--------------+
    else if (cmdName == "RemoveColumnsWithStops")
    {
      SiteContainer* sites = dynamic_cast<SiteContainer*>(sequences);
      if (!sites)
      {
        throw Exception("'RemoveColumnsWithStops' can only be used on alignment. You may consider using the 'CoerceToAlignment' command.");
      }
      if (!gCode.get()) {
        string codeDesc = ApplicationTools::getStringParameter("genetic_code", bppseqman.getParams(), "Standard", "", true, 1);
        ApplicationTools::displayResult("Genetic Code", codeDesc);
        gCode.reset(SequenceApplicationTools::getGeneticCode(codonAlphabet->getNucleicAlphabet(), codeDesc));
      }

      for (size_t i = sites->getNumberOfSites(); i > 0; i--)
      {
        if (CodonSiteTools::hasStop(sites->getSite(i-1), *gCode))
          sites->deleteSite(i - 1);
      }
    }

    // +---------+
    // | Get CDS |
    // +---------+
    else if (cmdName == "GetCDS")
    {
      if (!gCode.get()) {
        string codeDesc = ApplicationTools::getStringParameter("genetic_code", bppseqman.getParams(), "Standard", "", true, 1);
        ApplicationTools::displayResult("Genetic Code", codeDesc);
        gCode.reset(SequenceApplicationTools::getGeneticCode(codonAlphabet->getNucleicAlphabet(), codeDesc));
      }
      OrderedSequenceContainer* sc = 0;
      if (aligned) sc = new VectorSiteContainer(sequences->getAlphabet());
      else         sc = new VectorSequenceContainer(sequences->getAlphabet());
      for (size_t i = 0; i < sequences->getNumberOfSequences(); ++i)
      {
        BasicSequence seq = sequences->getSequence(i);
        size_t len = seq.size();
        SequenceTools::getCDS(seq, *gCode, false, true, true, false);
        if (aligned) {
          for (size_t c = seq.size(); c < len; ++c)
            seq.addElement(seq.getAlphabet()->getGapCharacterCode());
        }
        sc->addSequence(seq, false);
      }
      delete sequences;
      sequences = sc;