C++ TopSegmentIteratorConstPtr::getTopSegment方法代码示例

void TopSegmentSequenceTest::checkCallBack(AlignmentConstPtr alignment)
{
  const Genome* ancGenome = alignment->openGenome("Anc0");
  TopSegmentIteratorConstPtr tsIt = ancGenome->getTopSegmentIterator(100);
  CuAssertTrue(_testCase, tsIt->getTopSegment()->getStartPosition() == 500);
  CuAssertTrue(_testCase, tsIt->getTopSegment()->getLength() == 9);
  string seq;
  tsIt->getString(seq);
  CuAssertTrue(_testCase, seq == "CACACATTC");
  tsIt->toReverse();
  tsIt->getString(seq);
  CuAssertTrue(_testCase, seq == "GAATGTGTG");
}

void SummarizeMutations::subsAndGapInserts(
  GappedTopSegmentIteratorConstPtr gappedTop, MutationsStats& stats)
{
  assert(gappedTop->getReversed() == false);
  hal_size_t numGaps = gappedTop->getNumGaps();
  if (numGaps > 0)
  {
    stats._gapInsertionLength.add(gappedTop->getNumGapBases(), numGaps);
  }

  string parent, child;
  TopSegmentIteratorConstPtr l = gappedTop->getLeft();
  TopSegmentIteratorConstPtr r = gappedTop->getRight();
  BottomSegmentIteratorConstPtr p = 
     l->getTopSegment()->getGenome()->getParent()->getBottomSegmentIterator();

  for (TopSegmentIteratorConstPtr i = l->copy(); 
       i->getTopSegment()->getArrayIndex() <= 
          r->getTopSegment()->getArrayIndex();
       i->toRight())
  {
    if (i->hasParent())
    {
      p->toParent(i);
      i->getString(child);
      p->getString(parent);
      assert(child.length() == parent.length());
      for (size_t j = 0; j < child.length(); ++j)
      {
        if (isTransition(child[j], parent[j]))
        {
          ++stats._transitions;
          ++stats._subs;
        }
        else if (isTransversion(child[j], parent[j]))
        {
          ++stats._transversions;
          ++stats._subs;
        }
        else if (isSubstitution(child[j], parent[j]))
        {
          ++stats._subs;
        }
        else if (!isMissingData(child[j]) && !isMissingData(parent[j]))
        {
          ++stats._matches;
        }
      }
    }
  }
}

void DefaultGappedBottomSegmentIterator::toRightNextUngapped(
  TopSegmentIteratorConstPtr ts) const
{
  while (ts->hasParent() == false &&
         ts->getLength() <= _gapThreshold)
  {
    if ((!ts->getReversed() && ts->getTopSegment()->isLast()) ||
         (ts->getReversed() && ts->getTopSegment()->isFirst()))
    {
      break;
    }
    ts->toRight();
  }
}

void DefaultRearrangement::resetStatus(TopSegmentIteratorConstPtr topSegment)
{
    _id = Invalid;
    assert(topSegment.get());
    _genome = topSegment->getTopSegment()->getGenome();
    _parent = _genome->getParent();
    assert(_parent != NULL);

    _cur->setLeft(topSegment);
    _next->copy(_cur);
    _left->copy(_cur);
    _right->copy(_left);
    assert(_cur->getGapThreshold() == _gapThreshold);
    assert(_next->getGapThreshold() == _gapThreshold);
    assert(_left->getGapThreshold() == _gapThreshold);
    assert(_right->getGapThreshold() == _gapThreshold);
    assert(_leftParent->getGapThreshold() == _gapThreshold);
    assert(_rightParent->getGapThreshold() == _gapThreshold);
    assert(_curParent->getGapThreshold() == _gapThreshold);
    assert(_cur->getAtomic() == _atomic);
    assert(_next->getAtomic() == _atomic);
    assert(_left->getAtomic() == _atomic);
    assert(_right->getAtomic() == _atomic);
    assert(_leftParent->getAtomic() == _atomic);
    assert(_rightParent->getAtomic() == _atomic);
    assert(_curParent->getAtomic() == _atomic);
}

void TopSegmentStruct::compareTo(TopSegmentIteratorConstPtr it, 
                                 CuTest* testCase) const
{
  const TopSegment* seg = it->getTopSegment();
  CuAssertTrue(testCase, _length == seg->getLength());
  CuAssertTrue(testCase, _startPosition == seg->getStartPosition());
  CuAssertTrue(testCase, _nextParalogyIndex == seg->getNextParalogyIndex());
  CuAssertTrue(testCase, _parentIndex == seg->getParentIndex());
  CuAssertTrue(testCase, _bottomParseIndex == seg->getBottomParseIndex());
}

void hal::validateDuplications(const Genome* genome)
{
  const Genome* parent = genome->getParent();
  if (parent == NULL)
  {
    return;
  }
  TopSegmentIteratorConstPtr topIt = genome->getTopSegmentIterator();
  TopSegmentIteratorConstPtr endIt = genome->getTopSegmentEndIterator();
  vector<unsigned char> pcount(parent->getNumBottomSegments(), 0);
  for (; topIt != endIt; topIt->toRight())
  {
    if (topIt->hasParent())
    {
      if (pcount[topIt->getTopSegment()->getParentIndex()] < 250)
      {
        ++pcount[topIt->getTopSegment()->getParentIndex()];
      }
    }
  }
  for (topIt = genome->getTopSegmentIterator(); topIt != endIt; topIt->toRight())
  {
    if (topIt->hasParent())
    {
      size_t count = pcount[topIt->getTopSegment()->getParentIndex()];
      assert(count > 0);
      {
        if (topIt->hasNextParalogy() == false && count > 1)
        {
          stringstream ss;
          ss << "Top Segment " << topIt->getTopSegment()->getArrayIndex()
             << " in genome " << genome->getName() << " is not marked as a"
             << " duplication but it shares its parent " 
             << topIt->getTopSegment()->getArrayIndex() << " with at least " 
             << count - 1 << " other segments in the same genome";
          throw hal_exception(ss.str());
        }  
      }
    }
  }
}

void TopSegmentIteratorToSiteTest::checkGenome(const Genome* genome)
{
  TopSegmentIteratorConstPtr ti = genome->getTopSegmentIterator();
  for (hal_index_t pos = 0; 
       pos < (hal_index_t)genome->getSequenceLength(); ++pos)
  {
    ti->toSite(pos);
    CuAssertTrue(_testCase, ti->getStartPosition() == pos);
    CuAssertTrue(_testCase, ti->getLength() == 1);
    ti->toSite(pos, false);
    CuAssertTrue(_testCase, pos >= ti->getStartPosition() && 
                 pos < ti->getStartPosition() + (hal_index_t)ti->getLength());
    CuAssertTrue(_testCase, 
                 ti->getLength() == ti->getTopSegment()->getLength());
  }
}

void hal::validateSequence(const Sequence* sequence)
{
  // Verify that the DNA sequence doesn't contain funny characters
  DNAIteratorConstPtr dnaIt = sequence->getDNAIterator();
  hal_size_t length = sequence->getSequenceLength();
  for (hal_size_t i = 0; i < length; ++i)
  {
    char c = dnaIt->getChar();
    if (isNucleotide(c) == false)
    {
      stringstream ss;
      ss << "Non-nucleotide character discoverd at position " 
         << i << " of sequence " << sequence->getName() << ": " << c;
      throw hal_exception(ss.str());
    }
  }
  
  // Check the top segments
  if (sequence->getGenome()->getParent() != NULL)
  {
    hal_size_t totalTopLength = 0;
    TopSegmentIteratorConstPtr topIt = sequence->getTopSegmentIterator();
    hal_size_t numTopSegments = sequence->getNumTopSegments();
    for (hal_size_t i = 0; i < numTopSegments; ++i)
    {
      const TopSegment* topSegment = topIt->getTopSegment();
      validateTopSegment(topSegment);
      totalTopLength += topSegment->getLength();
      topIt->toRight();
    }
    if (totalTopLength != length)
    {
      stringstream ss;
      ss << "Sequence " << sequence->getName() << " has length " << length 
         << " but its top segments add up to " << totalTopLength;
      throw hal_exception(ss.str());
    }
  }

  // Check the bottom segments
  if (sequence->getGenome()->getNumChildren() > 0)
  {
    hal_size_t totalBottomLength = 0;
    BottomSegmentIteratorConstPtr bottomIt = 
       sequence->getBottomSegmentIterator();
    hal_size_t numBottomSegments = sequence->getNumBottomSegments();
    for (hal_size_t i = 0; i < numBottomSegments; ++i)
    {
      const BottomSegment* bottomSegment = bottomIt->getBottomSegment();
      validateBottomSegment(bottomSegment);
      totalBottomLength += bottomSegment->getLength();
      bottomIt->toRight();
    }
    if (totalBottomLength != length)
    {
      stringstream ss;
      ss << "Sequence " << sequence->getName() << " has length " << length 
         << " but its bottom segments add up to " << totalBottomLength;
      throw hal_exception(ss.str());
    }
  }
}

void hal::validateBottomSegment(const BottomSegment* bottomSegment)
{
  const Genome* genome = bottomSegment->getGenome();
  hal_index_t index = bottomSegment->getArrayIndex();
  if (index < 0 || index >= (hal_index_t)genome->getSequenceLength())
  {
    stringstream ss;
    ss << "Bottom segment out of range " << index << " in genome "
       << genome->getName();
    throw hal_exception(ss.str());
  }
  
  if (bottomSegment->getLength() < 1)
  {
    stringstream ss;
    ss << "Bottom segment " << index  << " in genome " << genome->getName()
       << " has length 0 which is not currently supported";
    throw hal_exception(ss.str());
  }

  hal_size_t numChildren = bottomSegment->getNumChildren();
  for (hal_size_t child = 0; child < numChildren; ++child)
  {
    const Genome* childGenome = genome->getChild(child);
    const hal_index_t childIndex = bottomSegment->getChildIndex(child);
    if (childGenome != NULL && childIndex != NULL_INDEX)
    {
      if (childIndex >= (hal_index_t)childGenome->getNumTopSegments())
      {
        stringstream ss;
        ss << "Child " << child << " index " <<childIndex << " of segment "
           << bottomSegment->getArrayIndex() << " out of range in genome "
           << childGenome->getName();
        throw hal_exception(ss.str());
      }
      TopSegmentIteratorConstPtr topSegmentIteratr = 
         childGenome->getTopSegmentIterator(childIndex);
      const TopSegment* childSegment = topSegmentIteratr->getTopSegment();
      if (childSegment->getLength() != bottomSegment->getLength())
      {
        stringstream ss;
        ss << "Child " << child << " with index " 
           << childSegment->getArrayIndex()
           << " and start position " << childSegment->getStartPosition() 
           << " and sequence " << childSegment->getSequence()->getName()
           << " has length " << childSegment->getLength()
           << " but parent with index " << bottomSegment->getArrayIndex() 
           << " and start position " << bottomSegment->getStartPosition()
           << " in sequence " << bottomSegment->getSequence()->getName() 
           << " has length " << bottomSegment->getLength();
        throw hal_exception(ss.str());
      }
      if (childSegment->getNextParalogyIndex() == NULL_INDEX &&
          childSegment->getParentIndex() != bottomSegment->getArrayIndex())
      {
        throw hal_exception("parent / child index mismatch (parent=" +
                            genome->getName() + " child=" +
                            childGenome->getName());
      }
      if (childSegment->getParentReversed() != 
          bottomSegment->getChildReversed(child))
      {
        throw hal_exception("parent / child reversal mismatch (parent=" +
                            genome->getName() + " child=" +
                            childGenome->getName());
      }
    }
  }

  const hal_index_t parseIndex = bottomSegment->getTopParseIndex();
  if (parseIndex == NULL_INDEX)
  {
    if (genome->getParent() != NULL)
    {
      stringstream ss;
      ss << "Bottom segment " << bottomSegment->getArrayIndex() << " in genome "
         << genome->getName() << " has null parse index";
      throw hal_exception(ss.str());
    }
  }
  else
  {
    if (parseIndex >= (hal_index_t)genome->getNumTopSegments())
    {
      stringstream ss;
      ss << "BottomSegment " << bottomSegment->getArrayIndex() << " in genome "
         << genome->getName() << " has parse index " << parseIndex 
         << " greater than the number of top segments, " 
         << (hal_index_t)genome->getNumTopSegments();
      throw hal_exception(ss.str());
    }
    TopSegmentIteratorConstPtr parseIterator = 
       genome->getTopSegmentIterator(parseIndex);
    const TopSegment* parseSegment = parseIterator->getTopSegment();
    hal_offset_t parseOffset = bottomSegment->getTopParseOffset();
    if (parseOffset >= parseSegment->getLength())
    {
      stringstream ss;
      ss << "BottomSegment " << bottomSegment->getArrayIndex() << " in genome "
         << genome->getName() << " has parse offset, " << parseOffset 
//.........这里部分代码省略.........

void hal::validateTopSegment(const TopSegment* topSegment)
{
  const Genome* genome = topSegment->getGenome();
  hal_index_t index = topSegment->getArrayIndex();
  if (index < 0 || index >= (hal_index_t)genome->getSequenceLength())
  {
    stringstream ss;
    ss << "Segment out of range " << index << " in genome "
       << genome->getName();
    throw hal_exception(ss.str());
  }

  if (topSegment->getLength() < 1)
  {
    stringstream ss;
    ss << "Top segment " << index  << " in genome " << genome->getName()
       << " has length 0 which is not currently supported";
    throw hal_exception(ss.str());
  }

  const Genome* parentGenome = genome->getParent();
  const hal_index_t parentIndex = topSegment->getParentIndex();
  if (parentGenome != NULL && parentIndex != NULL_INDEX)
  {
    if (parentIndex >= (hal_index_t)parentGenome->getNumBottomSegments())
    {
      stringstream ss;
      ss << "Parent index " << parentIndex << " of segment "
         << topSegment->getArrayIndex() << " out of range in genome "
         << parentGenome->getName();
      throw hal_exception(ss.str());
    }
    BottomSegmentIteratorConstPtr bottomSegmentIterator = 
       parentGenome->getBottomSegmentIterator(parentIndex);
    const BottomSegment* parentSegment = 
       bottomSegmentIterator->getBottomSegment();
    if (topSegment->getLength() != parentSegment->getLength())
    {
      stringstream ss;
      ss << "Parent length of segment " << topSegment->getArrayIndex() 
         << " in genome " << genome->getName() << " has length "
         << parentSegment->getLength() << " which does not match "
         << topSegment->getLength();
      throw hal_exception(ss.str());
    }
  }

  const hal_index_t parseIndex = topSegment->getBottomParseIndex();
  if (parseIndex == NULL_INDEX)
  {
    if (genome->getNumChildren() != 0)
    {
      stringstream ss;
      ss << "Top Segment " << topSegment->getArrayIndex() << " in genome "
         << genome->getName() << " has null parse index";
      throw hal_exception(ss.str());
    }
  }
  else
  {
    if (parseIndex >= (hal_index_t)genome->getNumBottomSegments())
    {
      stringstream ss;
      ss << "Top Segment " << topSegment->getArrayIndex() << " in genome "
         << genome->getName() << " has parse index out of range";
      throw hal_exception(ss.str());
    }
    hal_offset_t parseOffset = topSegment->getBottomParseOffset();
    BottomSegmentIteratorConstPtr bottomSegmentIterator =
       genome->getBottomSegmentIterator(parseIndex);
    const BottomSegment* parseSegment = 
       bottomSegmentIterator->getBottomSegment();
    if (parseOffset >= parseSegment->getLength())
    {
      stringstream ss;
      ss << "Top Segment " << topSegment->getArrayIndex() << " in genome "
         << genome->getName() << " has parse offset out of range";
      throw hal_exception(ss.str());
    }
    if ((hal_index_t)parseOffset + parseSegment->getStartPosition() != 
        topSegment->getStartPosition())
    {
      throw hal_exception("parse index broken in top segment in genome " +
                          genome->getName());
                          
    }
  }

  const hal_index_t paralogyIndex = topSegment->getNextParalogyIndex();
  if (paralogyIndex != NULL_INDEX)
  {
    TopSegmentIteratorConstPtr pti = 
       genome->getTopSegmentIterator(paralogyIndex);
    if (pti->getTopSegment()->getParentIndex() != topSegment->getParentIndex())
    {
      stringstream ss;
      ss << "Top segment " << topSegment->getArrayIndex() 
         << " has parent index "
         << topSegment->getParentIndex() << ", but next paraglog " 
         << topSegment->getNextParalogyIndex() << " has parent Index " 
//.........这里部分代码省略.........

void 
GappedSegmentSimpleIteratorTest::checkCallBack(AlignmentConstPtr alignment)
{
  const Genome* child = alignment->openGenome("child");
  const Genome* parent = alignment->openGenome("parent");

  GappedTopSegmentIteratorConstPtr gtsIt = 
     child->getGappedTopSegmentIterator(0, 9999999);
  GappedBottomSegmentIteratorConstPtr gbsIt = 
     parent->getGappedBottomSegmentIterator(0, 0, 9999999);
  GappedTopSegmentIteratorConstPtr gtsItRev = 
     child->getGappedTopSegmentIterator(0, 9999999);
  gtsItRev->toReverse();
  GappedBottomSegmentIteratorConstPtr gbsItRev = 
     parent->getGappedBottomSegmentIterator(0, 0, 9999999);
  gbsItRev->toReverse();

  for (size_t i = 0; i < child->getNumTopSegments(); ++i)
  {
    TopSegmentIteratorConstPtr tsIt = gtsIt->getLeft();
    CuAssertTrue(_testCase, tsIt->equals(gtsIt->getRight()));
    CuAssertTrue(_testCase, 
                 (size_t)tsIt->getTopSegment()->getArrayIndex() == i);
    gtsIt->toRight();

    BottomSegmentIteratorConstPtr bsIt = gbsIt->getLeft();
    CuAssertTrue(_testCase, bsIt->equals(gbsIt->getRight()));
    CuAssertTrue(_testCase, 
                 (size_t)bsIt->getBottomSegment()->getArrayIndex() == i);
    gbsIt->toRight();

    TopSegmentIteratorConstPtr tsItRev = gtsItRev->getLeft();
    CuAssertTrue(_testCase, tsItRev->equals(gtsItRev->getRight()));
    CuAssertTrue(_testCase, 
                 (size_t)tsItRev->getTopSegment()->getArrayIndex() == i);
    gtsItRev->toLeft();

    BottomSegmentIteratorConstPtr bsItRev = gbsItRev->getLeft();
    CuAssertTrue(_testCase, bsItRev->equals(gbsItRev->getRight()));
    CuAssertTrue(_testCase, 
                 (size_t)bsItRev->getBottomSegment()->getArrayIndex() == i);
    gbsItRev->toLeft();
  }

  gtsIt = child->getGappedTopSegmentIterator(
    child->getNumTopSegments() - 1, 9999999);
  gbsIt = parent->getGappedBottomSegmentIterator(
    child->getNumTopSegments() - 1, 0, 9999999);
  gtsItRev = child->getGappedTopSegmentIterator(
    child->getNumTopSegments() - 1, 9999999);
  gtsItRev->toReverse();
  gbsItRev = parent->getGappedBottomSegmentIterator(
    child->getNumTopSegments() - 1, 0, 9999999);
  gbsItRev->toReverse();

  for (hal_index_t i = child->getNumTopSegments() - 1; i >= 0; --i)
  {
    TopSegmentIteratorConstPtr tsIt = gtsIt->getLeft();
    CuAssertTrue(_testCase, tsIt->equals(gtsIt->getRight()));
    CuAssertTrue(_testCase, tsIt->getTopSegment()->getArrayIndex() == i);
    CuAssertTrue(_testCase, gtsIt->getReversed() == false);
    gtsIt->toLeft();

    BottomSegmentIteratorConstPtr bsIt = gbsIt->getLeft();
    CuAssertTrue(_testCase, bsIt->equals(gbsIt->getRight()));
    CuAssertTrue(_testCase, bsIt->getBottomSegment()->getArrayIndex() == i);
    CuAssertTrue(_testCase, gbsIt->getReversed() == false);
    gbsIt->toLeft();

    TopSegmentIteratorConstPtr tsItRev = gtsItRev->getLeft();
    CuAssertTrue(_testCase, tsItRev->equals(gtsItRev->getRight()));
    CuAssertTrue(_testCase, tsItRev->getTopSegment()->getArrayIndex() == i);
    CuAssertTrue(_testCase, gtsItRev->getReversed() == true);
    gtsItRev->toRight();

    BottomSegmentIteratorConstPtr bsItRev = gbsItRev->getLeft();
    CuAssertTrue(_testCase, bsItRev->equals(gbsItRev->getRight()));
    CuAssertTrue(_testCase, bsItRev->getBottomSegment()->getArrayIndex() == i);
    CuAssertTrue(_testCase, gbsItRev->getReversed() == true);
    gbsItRev->toRight();
  }

}