Changeset 4003


Ignore:
Timestamp:
Oct 15, 2020, 4:00:02 AM (8 months ago)
Author:
Peter
Message:

Add a constructor Aligner::Cigar(string) (closes #962)

Add a function cigar_table, which is a wrapper around bam_cigar_table
from htslib, translating operation characters in a CIGAR such as 'M'
or 'I' into the value as they are stored in a Cigar object (or a bam
entry). Function, bam_cigar_table, not available in older versions of
htslib or when building --without-htslib, so implement our own version
and decide at configure time whether using homebrew or call htslib.

Location:
trunk
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/m4/yat_check_htslib.m4

    r3890 r4003  
    11## $Id$
    22#
    3 # serial 7 (yat 0.18)
     3# serial 8 (yat 0.19)
    44#
    55#   Copyright (C) 2016, 2018, 2020 Peter Johansson
     
    4545  _YAT_CHECK_ZLIB
    4646  _YAT_CHECK_HTSLIB
     47  YAT_VAR_BAM_CIGAR_TABLE([
     48    AC_DEFINE([HAVE_BAM_CIGAR_TABLE], [1],
     49              [Define to 1 if bam_cigar_table is available])
     50  ])
    4751])
    4852_YAT_PROG_SAMTOOLS
     
    110114                [$2])
    111115]) # YAT_LIB_HTS
     116
     117
     118# YAT_VAR_BAM_CIGAR_TABLE([action-if-found], [action-if-not-found])
     119# =====================
     120# Check if variable bam_cigar_table is available.
     121AC_DEFUN([YAT_VAR_BAM_CIGAR_TABLE],
     122[
     123  AC_CACHE_CHECK([if bam_cigar_table is available],
     124    [yat_cv_var_bam_cigar_table],
     125    [AC_COMPILE_IFELSE([
     126        AC_LANG_PROGRAM([@%:@include <htslib/sam.h>],
     127                        [int8_t x = bam_cigar_table@<:@0@:>@;])
     128      ], [
     129        yat_cv_var_bam_cigar_table=yes
     130      ], [
     131        yat_cv_var_bam_cigar_table=no
     132      ])
     133    ])
     134  AS_VAR_IF([yat_cv_var_bam_cigar_table], [yes], [$1], [$2])
     135])
    112136
    113137
  • trunk/test/cigar.cc

    r3840 r4003  
    22
    33/*
    4   Copyright (C) 2014, 2019 Peter Johansson
     4  Copyright (C) 2014, 2019, 2020 Peter Johansson
    55
    66  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    2525
    2626#include "yat/utility/Aligner.h"
    27 
     27#include "yat/utility/Exception.h"
    2828#include "yat/utility/Matrix.h"
     29#include "yat/utility/utility.h"
    2930
    3031#include <sstream>
     
    8384  }
    8485  suite.out() << "\n";
     86  utility::Aligner::Cigar cig2(os.str());
     87  if (std::distance(cig.begin(), cig.end())!=
     88      std::distance(cig2.begin(), cig2.end()) ||
     89      !std::equal(cig.begin(), cig.end(), cig2.begin())) {
     90    suite.add(false);
     91    suite.err() << "error: incorrect construction from Cigar(string)\n";
     92    suite.err() << ": " << cig2 << "\nexpected:" << cig << "\n";
     93  }
     94  try {
     95    utility::Aligner::Cigar cig3("BANANA");
     96    suite.add(false);
     97    suite.err() << "error: Cigar(\"BANANA\") did not throw\n";
     98    suite.err() << "cigar: " << cig3 << "\n";
     99  }
     100  catch (utility::invalid_argument& e) {
     101    suite.out() << "expected exception: ";
     102    utility::print_what(e, suite.out());
     103    suite.out() << "\n";
     104  }
    85105
    86106  if (!suite.add(cig.length()==17))
  • trunk/yat/utility/Aligner.cc

    r3840 r4003  
    2828#include <algorithm>
    2929#include <cassert>
     30#include <cctype>
    3031#include <limits>
    3132#include <ostream>
    3233
     34#include <iostream>
    3335namespace theplu {
    3436namespace yat {
     
    158160  // ================ Cigar =======================
    159161
     162  Aligner::Cigar::Cigar(const std::string& s)
     163  {
     164    try {
     165      std::cerr << "Cigar: " << s << "\n";
     166      std::istringstream is(s);
     167      while (parse_element(is));
     168    }
     169    catch (invalid_argument& e) {
     170      std::ostringstream msg;
     171      msg << "Aligner::Cigar: " << s;
     172      std::throw_with_nested(invalid_argument(msg.str()));
     173    }
     174  }
     175
     176
    160177  Aligner::Cigar::const_iterator Aligner::Cigar::begin(void) const
    161178  {
     
    210227    assert(i<cigar_.size());
    211228    return bam_cigar_oplen(cigar_[i]);
     229  }
     230
     231
     232  bool Aligner::Cigar::parse_element(std::istream& is)
     233  {
     234    char c = is.get();
     235    if (!is.good())
     236      return false;
     237    if (!isdigit(c)) {
     238      std::string msg("not a digit: '");
     239      msg += c;
     240      msg += "'";
     241      throw invalid_argument(msg);
     242    }
     243    uint32_t len = c - '0';
     244
     245    while (is.get(c).good()) {
     246      if (isdigit(c))
     247        len = 10 * len + (c - '0');
     248      else {
     249        push_back(cigar_table(static_cast<uint8_t>(c)), len);
     250        return true;
     251      }
     252    }
     253    throw invalid_argument("Cigar::parse_element: missing operator");
     254    return true;
    212255  }
    213256
  • trunk/yat/utility/Aligner.h

    r3847 r4003  
    190190
    191191      /**
     192         Default constructor
     193       */
     194      Cigar(void) = default;
     195
     196      /**
     197         \param s string in same format as output from operator<<,
     198         i.e., an concatenated string of integer followed by one of
     199         "MIDNSHP=XB".
     200
     201         \since New in yat 0.19
     202       */
     203      explicit Cigar(const std::string& s);
     204
     205      /**
    192206         \return iterator pointing to first element
    193207
     
    302316    private:
    303317      std::deque<uint32_t> cigar_;
     318
     319      // parse a string of type 12M and push_back into cigar. Format
     320      // is is first an integer followed by a 1-char describing the
     321      // operation (see Cigar.h)
     322      bool parse_element(std::istream& is);
    304323
    305324      // calculate length only counting operations whose type is set
  • trunk/yat/utility/Cigar.h

    r3999 r4003  
    8585#endif // end of backport
    8686
     87namespace theplu {
     88namespace yat {
     89namespace utility {
     90
     91  /**
     92     The inverse of BAM_CIGAR_STR[x], where BAM_CIGAR_STR is
     93     "MIDNSHP=XB".
     94
     95     \throws if input is not one of "MIDNSHP=XB".
     96
     97     \since New in yat 0.19
     98  */
     99  int8_t cigar_table(unsigned char c);
     100
     101}}} // end of namespace utility, yat and theplu
     102
    87103#endif // end of file
  • trunk/yat/utility/Exception.cc

    r2919 r4003  
    3535namespace utility {
    3636
     37  invalid_argument::invalid_argument(const std::string& message)
     38    : std::invalid_argument(message)
     39  {}
     40
     41
    3742  runtime_error::runtime_error(std::string message)
    3843    : std::runtime_error(message)
  • trunk/yat/utility/Exception.h

    r2526 r4003  
    3131namespace utility {
    3232
    33  
     33  /**
     34     \brief Class used for all invalid arguments detected within yat library.
     35   */
     36  class invalid_argument : public std::invalid_argument
     37  {
     38  public:
     39    /**
     40       \brief Constructor
     41
     42       \param message message to be displayed using function what().
     43     */
     44    invalid_argument(const std::string& message);
     45  };
     46
    3447
    3548  /**
  • trunk/yat/utility/Makefile.am

    r3999 r4003  
    2323  yat/utility/Alignment.cc \
    2424  yat/utility/BasicMatrix.cc \
     25  yat/utility/Cigar.cc \
    2526  yat/utility/ColumnStream.cc \
    2627  yat/utility/CommandLine.cc \
Note: See TracChangeset for help on using the changeset viewer.