Changeset 3031


Ignore:
Timestamp:
Apr 27, 2013, 1:24:13 PM (9 years ago)
Author:
Peter
Message:

new function BamRead::sequence(2). refs #746

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/bam.cc

    r3028 r3031  
    120120}
    121121
     122
     123void test_sequence(test::Suite& suite, const BamRead& b)
     124{
     125  suite.out() << "test sequence\n";
     126  BamRead bam(b);
     127  std::string seq = b.sequence();
     128  std::vector<uint8_t> qual;
     129  qual.reserve(seq.size());
     130  for (size_t i=0; i<seq.size(); ++i)
     131    qual.push_back(bam.qual(i));
     132  assert(seq.size()==100);
     133  bam.sequence(seq, qual);
     134
     135  if (bam.sequence()!=seq) {
     136    suite.add(false);
     137    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
     138    suite.err() << "expected: " << seq << "\n";
     139  }
     140
     141  // just for consistency
     142  std::vector<uint32_t> cig;
     143  bam.cigar(cig);
     144  bam.core().flag &= ~BAM_FUNMAP;
     145
     146  // trim off one base
     147  seq.resize(seq.size()-1);
     148  qual.resize(qual.size()-1);
     149  bam.sequence(seq, qual);
     150
     151  if (bam.sequence()!=seq) {
     152    suite.add(false);
     153    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
     154    suite.err() << "expected: " << seq << "\n";
     155  }
     156
     157  // extend sequence
     158  seq.resize(120, 'A');
     159  qual.resize(120, 0);
     160  bam.sequence(seq, qual);
     161  if (bam.sequence()!=seq) {
     162    suite.add(false);
     163    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
     164    suite.err() << "expected: " << seq << "\n";
     165  }
     166}
     167
     168
    122169void test1(test::Suite& suite)
    123170{
     
    165212  test_cigar(suite, bam, in.header());
    166213  test_aux(suite, bam);
    167 }
    168 #endif
     214  test_sequence(suite, bam);
     215}
     216#endif
  • trunk/yat/omic/BamRead.cc

    r3030 r3031  
    206206
    207207
     208  void BamRead::reserve(int n)
     209  {
     210    assert(bam_);
     211    if (bam_->m_data >= n)
     212      return;
     213    // at least double capacity
     214    bam_->m_data = std::max(bam_->m_data<<1, n);
     215    bam_->data = (uint8_t*)realloc(bam_->data, bam_->m_data);
     216    if (!bam_->data) {
     217      throw utility::runtime_error("BamRead::reserve failed");
     218    }
     219  }
     220
     221
    208222  std::string BamRead::sequence(void) const
    209223  {
     
    221235
    222236
     237  void BamRead::sequence(const std::string& seq,
     238                         const std::vector<uint8_t>& qual)
     239  {
     240    assert(seq.size()==qual.size());
     241    // data consist of
     242    // |-- qname -->-- cigar -->-- seq -->-- qual -->-- aux -->
     243
     244    uint8_t* new_aux = bam1_seq(bam_) + (seq.size()+1)/2 + qual.size();
     245
     246    int len = new_aux + bam_->l_aux - bam_->data;
     247    reserve(len);
     248    // move aux to its new location
     249    memmove(static_cast<void*>(new_aux),
     250            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
     251    // copy sequence
     252    for (size_t i=0; i<seq.size(); ++i)
     253      sequence(i, bam_nt16_table[seq[i]]);
     254    core().l_qseq = seq.size();
     255    // copy quality
     256    if (qual.size()) // FIXME is this needed or does memcpy work with n=0?
     257      memcpy(static_cast<void*>(bam1_qual(bam_)),
     258             static_cast<const void*>(&qual[0]), qual.size());
     259  }
     260
     261
    223262  void BamRead::sequence(size_t i, uint8_t seq)
    224263  {
  • trunk/yat/omic/BamRead.h

    r3029 r3031  
    289289
    290290    /**
     291       \brief set sequence and quality
     292
     293       \since New in yat 0.11
     294     */
     295    void sequence(const std::string& seq, const std::vector<uint8_t>& qual);
     296
     297    /**
    291298       \see core().l_qseq
    292299     */
     
    306313
    307314  private:
     315    // ensure capacity of data pointer is (at least) n
     316    void reserve(int n);
     317
    308318    bam1_t* bam_;
    309319
Note: See TracChangeset for help on using the changeset viewer.