Changeset 3351


Ignore:
Timestamp:
Nov 21, 2014, 1:20:55 PM (8 years ago)
Author:
Peter
Message:

make BamRead? compile with htslib. refs #794

Location:
trunk/yat/omic
Files:
2 edited

Legend:

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

    r3324 r3351  
    2525#include "config_bam.h"
    2626
    27 #include YAT_BAM_HEADER
     27#include YAT_SAM_HEADER
    2828
    2929#include "yat/utility/Exception.h"
     
    4141#endif
    4242
     43// these #defines are no longer available in htslib, define them using
     44// htslib version
     45#ifndef bam1_seq
     46#define bam1_seq(b) (bam_get_seq(b))
     47#endif
     48
     49#ifndef bam1_seqi
     50#define bam1_seqi(b, i) (bam_seqi(b, i))
     51#endif
     52
     53#ifndef bam1_cigar
     54#define bam1_cigar(b) (bam_get_cigar(b))
     55#endif
     56
     57#ifndef bam1_qname
     58#define bam1_qname(b) (bam_get_qname(b))
     59#endif
     60
     61#ifndef bam1_qual
     62#define bam1_qual(b) (bam_get_qual(b))
     63#endif
     64
    4365namespace theplu {
    4466namespace yat {
     
    7395  {
    7496    assert(bam_);
     97#if YAT_HAVE_HTSLIB
     98    return bam_get_aux(bam_);
     99#else
    75100    return bam1_aux(bam_);
     101#endif
    76102  }
    77103
     
    107133  {
    108134    assert(bam_);
     135#if YAT_HAVE_HTSLIB
     136    return bam_get_l_aux(bam_);
     137#else
    109138    return bam_->l_aux;
     139#endif
    110140  }
    111141
     
    143173  const uint32_t* BamRead::cigar(void) const
    144174  {
     175#if YAT_HAVE_HTSLIB
     176    return bam_get_cigar(bam_);
     177#else
    145178    return bam1_cigar(bam_);
     179#endif
    146180  }
    147181
     
    150184  {
    151185    assert(i<core().n_cigar);
    152     return bam1_cigar(bam_)[i];
     186    return cigar()[i];
    153187  }
    154188
     
    172206    std::ostringstream os;
    173207    for (uint32_t i=0; i<core().n_cigar; ++i) {
    174       os << bam_cigar_oplen(bam1_cigar(bam_)[i])
    175          << bam_cigar_opchr(bam1_cigar(bam_)[i]);
     208      os << bam_cigar_oplen(cigar(i))
     209         << bam_cigar_opchr(cigar(i));
    176210    }
    177211    return os.str();
     
    192226  {
    193227    int offset = 4*c.size() - 4*core().n_cigar;
    194     reserve(bam_->data_len + offset);
     228   
     229    reserve(data_size() + offset);
    195230    /*
    196231      data come as qname, cigar, rest
     
    199234    if (offset)
    200235      memmove(bam1_seq(bam_)+offset, bam1_seq(bam_),
    201               bam_->data_len - 4*core().n_cigar - core().l_qname);
     236              data_size() - 4*core().n_cigar - core().l_qname);
    202237
    203238    // copy new cigar
     
    205240      bam1_cigar(bam_)[i] = c[i];
    206241
    207     bam_->data_len += offset;
     242    data_size() += offset;
    208243    core().n_cigar = c.size();
    209     assert(bam_->data_len <= bam_->m_data);
     244    assert(data_size() <= bam_->m_data);
    210245  }
    211246
     
    219254
    220255
     256  int& BamRead::data_size(void)
     257  {
     258#if YAT_HAVE_HTSLIB
     259    return bam_->l_data;
     260#else
     261    return bam_->data_len;
     262#endif
     263  }
     264
     265
     266  const int& BamRead::data_size(void) const
     267  {
     268#if YAT_HAVE_HTSLIB
     269    return bam_->l_data;
     270#else
     271    return bam_->data_len;
     272#endif
     273  }
     274
     275
    221276  int32_t BamRead::end(void) const
    222277  {
     
    230285    return bam_->core.flag;
    231286  }
     287
    232288
    233289  const char* BamRead::name(void) const
     
    239295    int offset = n.size() + 1 - core().l_qname;
    240296    assert(bam_);
    241     reserve(bam_->data_len + offset);
     297    reserve(data_size() + offset);
    242298    // move remaining data
    243299    if (offset)
    244300      memmove(bam_->data + core().l_qname + offset,
    245301              bam_->data + core().l_qname,
    246               bam_->data_len - core().l_qname);
     302              data_size() - core().l_qname);
    247303    core().l_qname += offset;
    248304    // copy new name
    249305    memcpy(bam1_qname(bam_), n.c_str(), n.size()+1);
    250     bam_->data_len += offset;
    251     assert(bam_->data_len <= bam_->m_data);
     306    data_size() += offset;
     307    assert(data_size() <= bam_->m_data);
    252308  }
    253309
     
    278334  {
    279335    std::string result(sequence_length(), ' ');
    280 #ifdef HAVE_BAM_NT16_REV_TABLE
     336#if HAVE_BAM_NT16_REV_TABLE
    281337    for (size_t i=0; i<result.size(); ++i)
    282338      result[i] = bam_nt16_rev_table[sequence(i)];
     
    298354
    299355    int aux_pos = bam1_seq(bam_) - bam_->data + (seq.size()+1)/2 + qual.size();
    300     int len = aux_pos + bam_->l_aux;
     356    int len = aux_pos + aux_size();
    301357    reserve(len);
    302358    // move aux to its new location
    303359    memmove(static_cast<void*>(bam_->data+aux_pos),
    304             static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
     360            static_cast<const void*>(aux()), aux_size());
    305361    // copy sequence
    306362    for (size_t i=0; i<seq.size(); ++i)
     363#if YAT_HAVE_HTSLIB
     364      sequence(i, seq_nt16_table[static_cast<size_t>(seq[i])]);
     365#else
    307366      sequence(i, bam_nt16_table[static_cast<size_t>(seq[i])]);
     367#endif
    308368    core().l_qseq = seq.size();
    309369    // copy quality
     
    311371      memcpy(static_cast<void*>(bam1_qual(bam_)),
    312372             static_cast<const void*>(&qual[0]), qual.size());
    313     bam_->data_len = len;
    314     assert(bam_->data_len <= bam_->m_data);
     373    data_size() = len;
     374    assert(data_size() <= bam_->m_data);
    315375  }
    316376
  • trunk/yat/omic/BamRead.h

    r3350 r3351  
    2424
    2525#include "config_bam.h"
    26 #include YAT_SAM_HEADER
    2726
    2827#include "yat/utility/Aligner.h"
     
    3130#include "yat/utility/deprecate.h"
    3231
     32#include YAT_SAM_HEADER
     33
    3334#include <functional>
    3435#include <string>
     
    309310    friend class BamReadIterator;
    310311    uint32_t calend(const bam1_core_t *c, const uint32_t *cigar) const;
     312    // access data length, current length of data
     313    int& data_size(void);
     314    const int& data_size(void) const;
    311315  };
    312316
Note: See TracChangeset for help on using the changeset viewer.