Changeset 3200


Ignore:
Timestamp:
May 3, 2014, 1:36:31 PM (7 years ago)
Author:
Peter
Message:

implement first version of CIGAR (refs #785). Note, this revision does not compile without samtools (should be fixed shortly).

Location:
trunk
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/Makefile.am

    r3190 r3200  
    4747  test/bam_header2.test \
    4848  test/bam_pair_iterator.test \
     49  test/cigar.test \
    4950  test/codon.test test/commandline.test \
    5051  test/concept.test \
  • trunk/yat/utility/Aligner.cc

    r3194 r3200  
    2626#include "Matrix.h"
    2727
     28
     29// FIXME fix this so it works also --without-samtools
     30#include <yat/omic/config_bam.h>
     31#include YAT_BAM_HEADER
     32
     33// FIXME we should include yat/omic/BamRead.h so we get #defines for
     34// e.g. bam_cigar_op also if not included in bam.h
     35
    2836#include <cassert>
    2937#include <limits>
     38
     39#include <iostream> // debug
    3040
    3141namespace theplu {
     
    118128
    119129
     130  const Aligner::Cigar Aligner::cigar(size_t row, size_t column) const
     131  {
     132    Cigar cig;
     133    while (row || column) {
     134      switch (alignment(row, column)) {
     135      case(diagonal):
     136        assert(row);
     137        --row;
     138        assert(column);
     139        --column;
     140        cig.push_front(BAM_CMATCH);
     141        break;
     142      case (right):
     143        assert(column);
     144        --column;
     145        cig.push_front(BAM_CINS);
     146        break;
     147      case (down):
     148        assert(row);
     149        --row;
     150        cig.push_front(BAM_CDEL);
     151        break;
     152      default:
     153        return cig;
     154      }
     155    }
     156    return cig;
     157  }
     158
     159
    120160  Aligner::direction& Aligner::directions(size_t i, size_t j)
    121161  {
     
    123163  }
    124164
     165
     166  // ================ Cigar =======================
     167
     168  void Aligner::Cigar::clear(void)
     169  {
     170    cigar_.clear();
     171  }
     172
     173
     174  uint8_t Aligner::Cigar::op(size_t i) const
     175  {
     176    assert(i<cigar_.size());
     177    return bam_cigar_op(cigar_[i]);
     178  }
     179
     180
     181  char Aligner::Cigar::opchr(size_t i) const
     182  {
     183    assert(i<cigar_.size());
     184    return bam_cigar_opchr(cigar_[i]);
     185  }
     186
     187
     188  uint32_t Aligner::Cigar::oplen(size_t i) const
     189  {
     190    assert(i<cigar_.size());
     191    return bam_cigar_oplen(cigar_[i]);
     192  }
     193
     194
     195  void Aligner::Cigar::pop_back(uint32_t len)
     196  {
     197    while (len) {
     198      assert(size());
     199      size_t i = size()-1;
     200      if (len < oplen(i)) {
     201        cigar_[i] -= (len<<BAM_CIGAR_SHIFT);
     202        return;
     203      }
     204      else {
     205        len -= oplen(i);
     206        cigar_.pop_back();
     207      }
     208    }
     209  }
     210
     211
     212  void Aligner::Cigar::pop_front(uint32_t len)
     213  {
     214    while (len) {
     215      assert(size());
     216      if (len < oplen(0)) {
     217        cigar_[0] -= (len<<BAM_CIGAR_SHIFT);
     218        return;
     219      }
     220      else {
     221        len -= oplen(0);
     222        cigar_.pop_front();
     223      }
     224    }
     225  }
     226
     227
     228  void Aligner::Cigar::push_back(uint8_t op, uint32_t len)
     229  {
     230    if (cigar_.empty() || this->op(size()-1)!=op)
     231      cigar_.push_front(bam_cigar_gen(len, op));
     232    else
     233      cigar_[size()-1] += (len<<BAM_CIGAR_SHIFT);
     234  }
     235
     236
     237  void Aligner::Cigar::push_front(uint8_t op, uint32_t len)
     238  {
     239    if (cigar_.empty() || this->op(0)!=op)
     240      cigar_.push_front(bam_cigar_gen(len, op));
     241    else
     242      cigar_[0] += (len<<BAM_CIGAR_SHIFT);
     243  }
     244
     245
     246  size_t Aligner::Cigar::size(void) const
     247  {
     248    return cigar_.size();
     249  }
     250
     251
     252  uint32_t Aligner::Cigar::operator[](size_t i) const
     253  {
     254    assert(i<cigar_.size());
     255    return cigar_[i];
     256  }
     257
     258
     259  std::ostream& operator<<(std::ostream& os, const Aligner::Cigar& cigar)
     260  {
     261    for (size_t i=0; i<cigar.size(); ++i)
     262      os << cigar.oplen(i) << cigar.opchr(i);
     263    return os;
     264  }
     265
    125266}}} // of namespace utility, yat, and theplu
  • trunk/yat/utility/Aligner.h

    r3195 r3200  
    2727#include <boost/cstdint.hpp>
    2828
     29#include <deque>
    2930#include <cstddef>
    3031#include <iosfwd>
     
    142143
    143144      /**
     145         \return operation i
     146       */
     147      uint8_t op(size_t i) const;
     148
     149      /**
     150         \return character describing operation i
     151       */
     152      char opchr(size_t i) const;
     153
     154      /**
     155         \return length of element i
     156       */
     157      uint32_t oplen(size_t i) const;
     158
     159      /**
     160         Erase last operation.
     161       */
     162      void pop_back(uint32_t len=1);
     163
     164      /**
     165         Erase first operation.
     166       */
     167      void pop_front(uint32_t len=1);
     168
     169      /**
    144170         Add an operation \a op of length \a len at back of cigar.
    145171       */
     
    150176       */
    151177      void push_front(uint8_t op, uint32_t len=1);
    152 
    153       /**
    154          Erase last operation.
    155        */
    156       void pop_back(void);
    157 
    158       /**
    159          Erase first operation.
    160        */
    161       void pop_front(void);
    162 
    163       /**
    164          \return operation i
    165        */
    166       uint8_t op(size_t i) const;
    167 
    168       /**
    169          \return character describing operation i
    170        */
    171       char opchr(size_t i) const;
    172 
    173       /**
    174          \return length of element i
    175        */
    176       uint32_t oplen(size_t i) const;
    177178
    178179      /**
     
    186187      size_t size(void) const;
    187188    private:
     189      std::deque<uint32_t> cigar_;
     190
    188191      // using compiler generated copy
    189192      // Cigar(const Cigar& other);
Note: See TracChangeset for help on using the changeset viewer.