C++ Genome::setDimensions方法代码示例

void GenomeUpdateTest::createCallBack(Alignment *alignment) {
    hal_size_t alignmentSize = alignment->getNumGenomes();
    CuAssertTrue(_testCase, alignmentSize == 0);

    Genome *ancGenome = alignment->addRootGenome("AncGenome", 0);
    vector<Sequence::Info> seqVec(1);
    seqVec[0] = Sequence::Info("Sequence", 1000000, 5000, 700000);
    ancGenome->setDimensions(seqVec);
    seqVec[0] = Sequence::Info("Sequence", 10000005, 14000, 2000001);
    ancGenome->setDimensions(seqVec);
}

void TopSegmentSimpleIteratorTest::createCallBack(Alignment *alignment) {
    Genome *ancGenome = alignment->addRootGenome("Anc0", 0);
    size_t numChildren = 9;
    for (size_t i = 0; i < numChildren; ++i) {
        alignment->addLeafGenome("Leaf" + std::to_string(i), "Anc0", 0.1);
    }
    vector<Sequence::Info> seqVec(1);
    seqVec[0] = Sequence::Info("Sequence", 1000000, 5000, 10000);
    ancGenome->setDimensions(seqVec);

    CuAssertTrue(_testCase, ancGenome->getNumChildren() == numChildren);

    _topSegments.clear();
    for (size_t i = 0; i < ancGenome->getNumTopSegments(); ++i) {
        TopSegmentStruct topSeg;
        topSeg.setRandom();
        topSeg._length = ancGenome->getSequenceLength() / ancGenome->getNumTopSegments();
        topSeg._startPosition = i * topSeg._length;
        _topSegments.push_back(topSeg);
    }

    TopSegmentIteratorPtr tsIt = ancGenome->getTopSegmentIterator(0);
    for (size_t i = 0; not tsIt->atEnd(); tsIt->toRight(), ++i) {
        CuAssertTrue(_testCase, (size_t)tsIt->getTopSegment()->getArrayIndex() == i);
        _topSegments[i].applyTo(tsIt);
    }
}

void SequenceUpdateTest::createCallBack(Alignment *alignment) {
    hal_size_t alignmentSize = alignment->getNumGenomes();
    CuAssertTrue(_testCase, alignmentSize == 0);

    Genome *ancGenome = alignment->addRootGenome("AncGenome", 0);

    size_t numSequences = 1000;
    vector<Sequence::Info> seqVec;
    for (size_t i = 0; i < numSequences; ++i) {
        hal_size_t len = 1 + i * 5 + i;
        string name = "sequence" + std::to_string(i);
        seqVec.push_back(Sequence::Info(name, len, i, i * 2));
    }
    ancGenome->setDimensions(seqVec);
    alignment->closeGenome(ancGenome);
    ancGenome = alignment->openGenome("AncGenome");

    vector<Sequence::UpdateInfo> updateVec;
    for (size_t i = 0; i < numSequences / 2; ++i) {
        const Sequence *sequence = ancGenome->getSequence("sequence" + std::to_string(i));
        updateVec.push_back(Sequence::UpdateInfo(sequence->getName(), i * 7));
    }
    ancGenome->updateTopDimensions(updateVec);

    updateVec.clear();
    for (size_t i = 0; i < numSequences / 3; ++i) {
        const Sequence *sequence = ancGenome->getSequence("sequence" + std::to_string(i));
        updateVec.push_back(Sequence::UpdateInfo(sequence->getName(), i * 5));
    }
    ancGenome->updateBottomDimensions(updateVec);
}

void MappedSegmentMapDupeTest::createCallBack(AlignmentPtr alignment)
{
  vector<Sequence::Info> seqVec(1);
  
  BottomSegmentIteratorPtr bi;
  BottomSegmentStruct bs;
  TopSegmentIteratorPtr ti;
  TopSegmentStruct ts;
  
  // setup simple case were there is an edge from a parent to 
  // child and it is reversed
  Genome* parent = alignment->addRootGenome("parent");
  Genome* child1 = alignment->addLeafGenome("child1", "parent", 1);
  Genome* child2 = alignment->addLeafGenome("child2", "parent", 1);
  seqVec[0] = Sequence::Info("Sequence", 3, 0, 1);
  parent->setDimensions(seqVec);
  seqVec[0] = Sequence::Info("Sequence", 9, 3, 0);
  child1->setDimensions(seqVec);
  seqVec[0] = Sequence::Info("Sequence", 9, 3, 0);
  child2->setDimensions(seqVec);

  parent->setString("CCC");
  child1->setString("CCCTACGTG");
  child2->setString("CCCTACGTG");

  bi = parent->getBottomSegmentIterator();
  bs.set(0, 3);
  bs._children.push_back(pair<hal_size_t, bool>(0, true));
  bs._children.push_back(pair<hal_size_t, bool>(0, false));
  bs.applyTo(bi);
     
  ti = child1->getTopSegmentIterator();
  ts.set(0, 3, 0, true, NULL_INDEX, 1);
  ts.applyTo(ti);
  ti->toRight();
  ts.set(3, 3, 0, true, NULL_INDEX, 2);
  ts.applyTo(ti);
  ti->toRight();
  ts.set(6, 3, 0, true, NULL_INDEX, 0);
  ts.applyTo(ti);

  ti = child2->getTopSegmentIterator();
  ts.set(0, 3, 0, false);
  ts.applyTo(ti);
  ti->toRight();
  ts.set(3, 3, NULL_INDEX, true);
  ts.applyTo(ti);
  ti->toRight();
  ts.set(6, 3, NULL_INDEX, false);
  ts.applyTo(ti);
}

void LodExtract::writeDimensions(
  const map<const Sequence*, hal_size_t>& segmentCounts, 
  const string& parentName,
  const vector<string>& childNames)
{
  // initialize a dimensions list for each (input) genome
  map<const Genome*, vector<Sequence::Info> > dimMap;
  map<const Genome*, vector<Sequence::Info> >::iterator dimMapIt;
  vector<string> newGenomeNames = childNames;
  newGenomeNames.push_back(parentName);
 
  for (size_t i = 0; i < newGenomeNames.size(); ++i)
  {
    const Genome* inGenome = _inAlignment->openGenome(newGenomeNames[i]);
    pair<const Genome*, vector<Sequence::Info> > newEntry;
    newEntry.first = inGenome;
    
    // it's important we keep the sequences in the output genome
    // in the same order as the sequences in the input genome since
    // we always use global coordinates!
    SequenceIteratorConstPtr seqIt = inGenome->getSequenceIterator();
    SequenceIteratorConstPtr seqEnd = inGenome->getSequenceEndIterator();
    for (; seqIt != seqEnd; seqIt->toNext())
    {
      const Sequence* inSequence = seqIt->getSequence();
      map<const Sequence*, hal_size_t>::const_iterator segMapIt;
      segMapIt = segmentCounts.find(inSequence);
      // we skip empty sequences for now with below check
      if (segMapIt != segmentCounts.end())
      {
        vector<Sequence::Info>& segDims = newEntry.second;
        hal_size_t nTop = 
           inGenome->getName() == parentName ? 0 : segMapIt->second;
        hal_size_t nBot = 
           inGenome->getName() != parentName ? 0 : segMapIt->second;
        segDims.push_back(Sequence::Info(inSequence->getName(),
                                         inSequence->getSequenceLength(),
                                         nTop,
                                         nBot));
      }
    }

    // note potential bug here for genome with no data
    dimMap.insert(newEntry);
  }
  
  // now that we have the dimensions for each genome, update them in
  // the output alignment
  for (dimMapIt = dimMap.begin(); dimMapIt != dimMap.end(); ++dimMapIt)
  {
    Genome* newGenome = _outAlignment->openGenome(dimMapIt->first->getName());
    assert(newGenome != NULL);
    vector<Sequence::Info>& segDims = dimMapIt->second;
    // ROOT 
    if (newGenome->getName() == _outAlignment->getRootName())
    {
      assert(newGenome->getName() == parentName);
      newGenome->setDimensions(segDims, _keepSequences);
    }
    // LEAF
    else if (newGenome->getName() != parentName)
    {
      newGenome->setDimensions(segDims, _keepSequences);
    }
    // INTERNAL NODE
    else
    {
      vector<Sequence::UpdateInfo> updateInfo;
      for (size_t i = 0; i < segDims.size(); ++i)
      {
        updateInfo.push_back(
          Sequence::UpdateInfo(segDims[i]._name,
                               segDims[i]._numBottomSegments));
      }
      newGenome->updateBottomDimensions(updateInfo);
    }
  }
}