Ignore:
Timestamp:
Apr 27, 2013, 5:06:52 PM (9 years ago)
Author:
Peter
Message:

refs #746. Simplify cigar function using logic similar to sequence(2). Avoid bug as pointers to data get invalidated when calling reserve. Forgot to set data_len and avoid compiler warning using char as index to array

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/omic/BamRead.cc

    r3031 r3032  
    151151  void BamRead::cigar(const std::vector<uint32_t>& c)
    152152  {
    153     // special case when cigar size is not changed
    154     if (c.size() == core().n_cigar) {
    155       if (c.size())
    156         memcpy(bam1_cigar(bam_), &c[0], 4*c.size());
    157       return;
    158     }
    159 
    160     int m_data = bam_->m_data;
    161     int data_len = bam_->data_len + 4*c.size() - 4*core().n_cigar;
    162     // double capacity if needed
    163     if (m_data < data_len) {
    164       m_data = data_len;
    165       kroundup32(m_data);
    166     }
    167     uint8_t* data = (uint8_t*)calloc(m_data, 1);
     153    int offset = 4*c.size() - 4*core().n_cigar;
     154    reserve(bam_->data_len + offset);
    168155    /*
    169156      data comes as qname, cigar, rest
    170157     */
    171     // first copy qname
    172     memcpy(data, bam_->data, core().l_qname);
     158    // first move remaining data to new location
     159    if (offset)
     160      memmove(bam1_seq(bam_)+offset, bam1_seq(bam_),
     161              bam_->data_len - 4*core().n_cigar - core().l_qname);
     162
    173163    // copy new cigar
    174164    if (c.size())
    175165      // each cigar element is 4 bytes
    176       memcpy(data + core().l_qname, &c[0], 4*c.size());
    177     // copy remaining data
    178     memcpy(data + core().l_qname + 4*c.size(), bam1_seq(bam_),
    179            data_len - 4*c.size() - core().l_qname);
    180     bam_->m_data = m_data;
    181     bam_->data_len = data_len;
     166      memcpy(bam1_cigar(bam_), &c[0], 4*c.size());
     167    bam_->data_len += offset;
    182168    core().n_cigar = c.size();
    183     free(bam_->data);
    184     bam_->data = data;
     169    assert(bam_->data_len <= bam_->m_data);
    185170  }
    186171
     
    242227    // |-- qname -->-- cigar -->-- seq -->-- qual -->-- aux -->
    243228
    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;
     229    int aux_pos = bam1_seq(bam_) - bam_->data + (seq.size()+1)/2 + qual.size();
     230    int len = aux_pos + bam_->l_aux;
    247231    reserve(len);
    248232    // move aux to its new location
    249     memmove(static_cast<void*>(new_aux),
     233    memmove(static_cast<void*>(bam_->data+aux_pos),
    250234            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
    251235    // copy sequence
    252236    for (size_t i=0; i<seq.size(); ++i)
    253       sequence(i, bam_nt16_table[seq[i]]);
     237      sequence(i, bam_nt16_table[static_cast<size_t>(seq[i])]);
    254238    core().l_qseq = seq.size();
    255239    // copy quality
     
    257241      memcpy(static_cast<void*>(bam1_qual(bam_)),
    258242             static_cast<const void*>(&qual[0]), qual.size());
     243    bam_->data_len = len;
     244    assert(bam_->data_len <= bam_->m_data);
    259245  }
    260246
Note: See TracChangeset for help on using the changeset viewer.