yat/omic/BamRead.h

Code
Comments
Other
Rev Date Author Line
2883 03 Dec 12 peter 1 #ifndef theplu_yat_omic_bam_read
2883 03 Dec 12 peter 2 #define theplu_yat_omic_bam_read
2883 03 Dec 12 peter 3
2883 03 Dec 12 peter 4 // $Id$
2883 03 Dec 12 peter 5
2993 03 Mar 13 peter 6 /*
4359 23 Aug 23 peter 7   Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2023 Peter Johansson
2993 03 Mar 13 peter 8
2993 03 Mar 13 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
2993 03 Mar 13 peter 10
2993 03 Mar 13 peter 11   The yat library is free software; you can redistribute it and/or
2993 03 Mar 13 peter 12   modify it under the terms of the GNU General Public License as
2993 03 Mar 13 peter 13   published by the Free Software Foundation; either version 3 of the
2993 03 Mar 13 peter 14   License, or (at your option) any later version.
2993 03 Mar 13 peter 15
2993 03 Mar 13 peter 16   The yat library is distributed in the hope that it will be useful,
2993 03 Mar 13 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
2993 03 Mar 13 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2993 03 Mar 13 peter 19   General Public License for more details.
2993 03 Mar 13 peter 20
2993 03 Mar 13 peter 21   You should have received a copy of the GNU General Public License
3752 17 Oct 18 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
2993 03 Mar 13 peter 23 */
2993 03 Mar 13 peter 24
3295 25 Jul 14 peter 25 #include "yat/utility/Aligner.h"
3212 05 May 14 peter 26 // This file has to be included to keep compatibility with yat 0.11
3295 25 Jul 14 peter 27 #include "yat/utility/Cigar.h"
3295 25 Jul 14 peter 28 #include "yat/utility/deprecate.h"
3202 03 May 14 peter 29
3883 24 Mar 20 peter 30 #include <htslib/sam.h>
3351 21 Nov 14 peter 31
2883 03 Dec 12 peter 32 #include <functional>
2883 03 Dec 12 peter 33 #include <string>
2883 03 Dec 12 peter 34 #include <vector>
2883 03 Dec 12 peter 35
2883 03 Dec 12 peter 36 namespace theplu {
2883 03 Dec 12 peter 37 namespace yat {
2883 03 Dec 12 peter 38 namespace omic {
2883 03 Dec 12 peter 39
2883 03 Dec 12 peter 40   /**
2883 03 Dec 12 peter 41      \brief Class holding a bam query
2883 03 Dec 12 peter 42
2883 03 Dec 12 peter 43      This is a wrapper around bam1_t and most of its information is
2883 03 Dec 12 peter 44      held in the core struct. A BamRead is typically created from a
2883 03 Dec 12 peter 45      InBamFile.
2883 03 Dec 12 peter 46
2883 03 Dec 12 peter 47      \see samtools
2884 04 Dec 12 peter 48
2884 04 Dec 12 peter 49      \since New in yat 0.10
2883 03 Dec 12 peter 50    */
2883 03 Dec 12 peter 51   class BamRead
2883 03 Dec 12 peter 52   {
2883 03 Dec 12 peter 53   public:
2883 03 Dec 12 peter 54     /**
2883 03 Dec 12 peter 55        \brief default constructor
2883 03 Dec 12 peter 56
2883 03 Dec 12 peter 57        Constructed object contains no data and most operations are not
2883 03 Dec 12 peter 58        defined
2883 03 Dec 12 peter 59      */
2883 03 Dec 12 peter 60     BamRead(void);
2883 03 Dec 12 peter 61
2883 03 Dec 12 peter 62     /**
2883 03 Dec 12 peter 63        \brief Copy constructor
2883 03 Dec 12 peter 64      */
2883 03 Dec 12 peter 65     BamRead(const BamRead& other);
2883 03 Dec 12 peter 66
2883 03 Dec 12 peter 67     /**
3583 19 Jan 17 peter 68        \brief move constructor
3596 21 Jan 17 peter 69
3596 21 Jan 17 peter 70        \since new in yat 0.15
3583 19 Jan 17 peter 71      */
3691 14 Sep 17 peter 72     BamRead(BamRead&& other) noexcept;
3583 19 Jan 17 peter 73
3583 19 Jan 17 peter 74     /**
2883 03 Dec 12 peter 75        \brief Destructor
2883 03 Dec 12 peter 76      */
2883 03 Dec 12 peter 77     virtual ~BamRead(void);
2883 03 Dec 12 peter 78
2883 03 Dec 12 peter 79     /**
2883 03 Dec 12 peter 80        \brief assignment operator
2883 03 Dec 12 peter 81      */
2883 03 Dec 12 peter 82     const BamRead& operator=(const BamRead& rhs);
2883 03 Dec 12 peter 83
2883 03 Dec 12 peter 84     /**
3583 19 Jan 17 peter 85        \brief move assignment operator
3596 21 Jan 17 peter 86
3596 21 Jan 17 peter 87        \since new in yat 0.15
3583 19 Jan 17 peter 88      */
3583 19 Jan 17 peter 89     BamRead& operator=(BamRead&& rhs);
3583 19 Jan 17 peter 90
3583 19 Jan 17 peter 91     /**
2883 03 Dec 12 peter 92        \return pointer to auxiliary data
3029 21 Apr 13 peter 93
3029 21 Apr 13 peter 94        \see aux_size(void)
2883 03 Dec 12 peter 95      */
2883 03 Dec 12 peter 96     const uint8_t* aux(void) const;
2883 03 Dec 12 peter 97
2883 03 Dec 12 peter 98     /**
3029 21 Apr 13 peter 99        Use bam_aux2? functions (in samtools C api) to convert returned
3029 21 Apr 13 peter 100        pointer to corresponding type.
3029 21 Apr 13 peter 101
3029 21 Apr 13 peter 102        \return pointer to field associated with \a tag, NULL if \a tag
3029 21 Apr 13 peter 103        doesn't exist.
3029 21 Apr 13 peter 104
3029 21 Apr 13 peter 105        \see bam_aux_get
3029 21 Apr 13 peter 106
3028 21 Apr 13 peter 107        \since New in yat 0.11
3028 21 Apr 13 peter 108      */
3028 21 Apr 13 peter 109     const uint8_t* aux(const char tag[2]) const;
3028 21 Apr 13 peter 110
3028 21 Apr 13 peter 111     /**
3029 21 Apr 13 peter 112        \brief append a new tag to aux field
3029 21 Apr 13 peter 113
3029 21 Apr 13 peter 114        \param tag two-charcter tag to append
3029 21 Apr 13 peter 115        \param type describes which type and can be 'iIsScCdfAZH'
3029 21 Apr 13 peter 116        \param len length of data
3029 21 Apr 13 peter 117        \param data pointer to data
3029 21 Apr 13 peter 118
3028 21 Apr 13 peter 119        \since New in yat 0.11
3028 21 Apr 13 peter 120
3029 21 Apr 13 peter 121        \see SAM specification
3028 21 Apr 13 peter 122      */
3028 21 Apr 13 peter 123     void aux_append(const char tag[2], char type, int len, uint8_t* data);
3028 21 Apr 13 peter 124
3028 21 Apr 13 peter 125     /**
3029 21 Apr 13 peter 126        \brief remove a tag in aux field
3029 21 Apr 13 peter 127
3029 21 Apr 13 peter 128        \throw utility::runtime_error if \a tag is not present in read
3029 21 Apr 13 peter 129
3028 21 Apr 13 peter 130        \since New in yat 0.11
3028 21 Apr 13 peter 131      */
3028 21 Apr 13 peter 132     void aux_del(const char tag[2]);
3028 21 Apr 13 peter 133
3028 21 Apr 13 peter 134     /**
3029 21 Apr 13 peter 135        \return length of aux field
3029 21 Apr 13 peter 136
3028 21 Apr 13 peter 137        \since New in yat 0.11
3028 21 Apr 13 peter 138      */
3028 21 Apr 13 peter 139     int aux_size(void) const;
3028 21 Apr 13 peter 140
3028 21 Apr 13 peter 141     /**
2883 03 Dec 12 peter 142        \brief access core data struct
2883 03 Dec 12 peter 143
2883 03 Dec 12 peter 144        \see samtools C api documentaion
2883 03 Dec 12 peter 145      */
2883 03 Dec 12 peter 146     const bam1_core_t& core(void) const;
2883 03 Dec 12 peter 147
2883 03 Dec 12 peter 148     /**
2883 03 Dec 12 peter 149        \brief access core data struct
2883 03 Dec 12 peter 150
2883 03 Dec 12 peter 151        \see samtools C api documentaion
2883 03 Dec 12 peter 152      */
2883 03 Dec 12 peter 153     bam1_core_t& core(void);
2883 03 Dec 12 peter 154
2883 03 Dec 12 peter 155     /**
2883 03 Dec 12 peter 156        \brief access CIGAR array
2883 03 Dec 12 peter 157
2883 03 Dec 12 peter 158        In each element the lower 4 bits gives a CIGAR operation and
2883 03 Dec 12 peter 159        the upper 28 bits keep the length.
3380 25 Feb 15 peter 160
3380 25 Feb 15 peter 161        The size of CIGAR array is accessed via core().n_cigar.
2883 03 Dec 12 peter 162      */
2883 03 Dec 12 peter 163     const uint32_t* cigar(void) const;
2883 03 Dec 12 peter 164
2883 03 Dec 12 peter 165     /**
2883 03 Dec 12 peter 166        \return \a i th element of CIGAR array
2883 03 Dec 12 peter 167      */
2883 03 Dec 12 peter 168     uint32_t cigar(size_t i) const;
2883 03 Dec 12 peter 169
2883 03 Dec 12 peter 170     /**
2883 03 Dec 12 peter 171        \return \a i th CIGAR operation
2883 03 Dec 12 peter 172      */
2883 03 Dec 12 peter 173     uint32_t cigar_op(size_t i) const;
2883 03 Dec 12 peter 174
2883 03 Dec 12 peter 175     /**
2883 03 Dec 12 peter 176        \return length of \a i th CIGAR element
2883 03 Dec 12 peter 177      */
2883 03 Dec 12 peter 178     uint32_t cigar_oplen(size_t i) const;
2883 03 Dec 12 peter 179
2883 03 Dec 12 peter 180     /**
2883 03 Dec 12 peter 181        Translate CIGAR array to a string such as '72M3S'
2883 03 Dec 12 peter 182      */
2883 03 Dec 12 peter 183     std::string cigar_str(void) const;
2883 03 Dec 12 peter 184
2883 03 Dec 12 peter 185     /**
2883 03 Dec 12 peter 186        \brief set CIGAR
2883 03 Dec 12 peter 187
2883 03 Dec 12 peter 188        \param c new cigar
3295 25 Jul 14 peter 189
3295 25 Jul 14 peter 190        \deprecated Provided for backward compatibility with 0.11
3295 25 Jul 14 peter 191        API. Use cigar(const utility::Aligner::Cigar&) instead.
2883 03 Dec 12 peter 192     */
3295 25 Jul 14 peter 193     void cigar(const std::vector<uint32_t>& c) YAT_DEPRECATE;
2883 03 Dec 12 peter 194
2883 03 Dec 12 peter 195     /**
3295 25 Jul 14 peter 196        \brief set CIGAR
3295 25 Jul 14 peter 197
3295 25 Jul 14 peter 198        \param cigar new cigar
3295 25 Jul 14 peter 199
3295 25 Jul 14 peter 200        \since new in yat 0.12
3295 25 Jul 14 peter 201     */
3295 25 Jul 14 peter 202     void cigar(const utility::Aligner::Cigar& cigar);
3295 25 Jul 14 peter 203
3295 25 Jul 14 peter 204     /**
3295 25 Jul 14 peter 205        \brief rightmost coordinate
3295 25 Jul 14 peter 206
3295 25 Jul 14 peter 207        Coordinate is 0-based, i.e., the end is one passed the last
3295 25 Jul 14 peter 208        matching position.
3295 25 Jul 14 peter 209
3295 25 Jul 14 peter 210        \see bam_calend
2883 03 Dec 12 peter 211      */
3295 25 Jul 14 peter 212     int32_t end(void) const;
2883 03 Dec 12 peter 213
2883 03 Dec 12 peter 214     /**
3295 25 Jul 14 peter 215        \brief bitwise flag
3026 06 Apr 13 peter 216
3295 25 Jul 14 peter 217        \see Preprocessor defines BAM_F*
3026 06 Apr 13 peter 218
3295 25 Jul 14 peter 219        \since implemented since yat 0.12
3026 06 Apr 13 peter 220      */
3295 25 Jul 14 peter 221     uint16_t flag(void) const;
3026 06 Apr 13 peter 222
3026 06 Apr 13 peter 223     /**
3732 11 Apr 18 peter 224        \return \c true if each bit set in x is also set in flag()
3732 11 Apr 18 peter 225
3732 11 Apr 18 peter 226        \since new in yat 0.16
3732 11 Apr 18 peter 227      */
3732 11 Apr 18 peter 228     bool flag_bit_is_set(uint16_t x) const;
3732 11 Apr 18 peter 229
3732 11 Apr 18 peter 230     /**
3731 11 Apr 18 peter 231        \return bits set in both \a x and flag()
3731 11 Apr 18 peter 232
3731 11 Apr 18 peter 233        \since new in yat 0.16
3731 11 Apr 18 peter 234      */
3731 11 Apr 18 peter 235     uint16_t get_flag_bit(uint16_t x) const;
3731 11 Apr 18 peter 236
3731 11 Apr 18 peter 237     /**
3731 11 Apr 18 peter 238        For each bit set in x, set bit in flag.
3731 11 Apr 18 peter 239
3731 11 Apr 18 peter 240        flag |= x
3731 11 Apr 18 peter 241
3731 11 Apr 18 peter 242        \since new in yat 0.16
3731 11 Apr 18 peter 243      */
3731 11 Apr 18 peter 244     void set_flag_bit(uint16_t x);
3731 11 Apr 18 peter 245
3731 11 Apr 18 peter 246     /**
3731 11 Apr 18 peter 247        For each bit set in x, unset bit in flag.
3731 11 Apr 18 peter 248
3731 11 Apr 18 peter 249        flag &= ~x
3731 11 Apr 18 peter 250
3731 11 Apr 18 peter 251        \since new in yat 0.16
3731 11 Apr 18 peter 252      */
3731 11 Apr 18 peter 253     void unset_flag_bit(uint16_t x);
3731 11 Apr 18 peter 254
3731 11 Apr 18 peter 255     /**
3731 11 Apr 18 peter 256        For each bit set in x, flip bit in flag.
3731 11 Apr 18 peter 257
3731 11 Apr 18 peter 258        flag ^= x
3731 11 Apr 18 peter 259
3731 11 Apr 18 peter 260        \since new in yat 0.16
3731 11 Apr 18 peter 261      */
3731 11 Apr 18 peter 262     void flip_flag_bit(uint16_t x);
3731 11 Apr 18 peter 263
3731 11 Apr 18 peter 264     /**
3295 25 Jul 14 peter 265        \brief leftmost position for mate
3295 25 Jul 14 peter 266      */
3295 25 Jul 14 peter 267     int32_t mpos(void) const;
3295 25 Jul 14 peter 268
3295 25 Jul 14 peter 269     /**
3295 25 Jul 14 peter 270        \brief Chromosome ID for mate
3295 25 Jul 14 peter 271      */
3295 25 Jul 14 peter 272     int32_t mtid(void) const;
3295 25 Jul 14 peter 273
3295 25 Jul 14 peter 274     /**
2883 03 Dec 12 peter 275        \return query name
2883 03 Dec 12 peter 276
2883 03 Dec 12 peter 277        Length of array is described by core().l_qname
2883 03 Dec 12 peter 278      */
2883 03 Dec 12 peter 279     const char* name(void) const;
2883 03 Dec 12 peter 280
2883 03 Dec 12 peter 281     /**
3055 16 Jun 13 peter 282        \brief modify name
3055 16 Jun 13 peter 283
3055 16 Jun 13 peter 284        \since New in yat 0.11
3055 16 Jun 13 peter 285     */
3055 16 Jun 13 peter 286     void name(const std::string& n);
3055 16 Jun 13 peter 287
3055 16 Jun 13 peter 288     /**
3969 13 Aug 20 peter 289        \brief 0-based leftmost coordinate
2883 03 Dec 12 peter 290      */
2883 03 Dec 12 peter 291     int32_t pos(void) const;
2883 03 Dec 12 peter 292
2883 03 Dec 12 peter 293     /**
3295 25 Jul 14 peter 294        \return Quality of base \a i
2883 03 Dec 12 peter 295      */
3295 25 Jul 14 peter 296     uint8_t qual(size_t i) const;
2883 03 Dec 12 peter 297
2883 03 Dec 12 peter 298     /**
3295 25 Jul 14 peter 299        \brief set quality of a base
2883 03 Dec 12 peter 300
3295 25 Jul 14 peter 301        \param i base to modify
3295 25 Jul 14 peter 302        \param q new quality
3295 25 Jul 14 peter 303
3295 25 Jul 14 peter 304        \since New in yat 0.11
2883 03 Dec 12 peter 305      */
3295 25 Jul 14 peter 306     void qual(size_t i, uint8_t q);
2883 03 Dec 12 peter 307
2883 03 Dec 12 peter 308     /**
2883 03 Dec 12 peter 309        Each character in returned sequence is either A, C, G, T, or N.
2883 03 Dec 12 peter 310
2883 03 Dec 12 peter 311        \return sequence
2883 03 Dec 12 peter 312      */
2883 03 Dec 12 peter 313     std::string sequence(void) const;
2883 03 Dec 12 peter 314
2883 03 Dec 12 peter 315     /**
2883 03 Dec 12 peter 316        4-bit integer describing base \a index
2883 03 Dec 12 peter 317
3843 13 Sep 19 peter 318        \see nt16_str
2883 03 Dec 12 peter 319      */
2883 03 Dec 12 peter 320     uint8_t sequence(size_t index) const;
2883 03 Dec 12 peter 321
2883 03 Dec 12 peter 322     /**
3395 23 Mar 15 peter 323        \return A, C, G, T, or N.
3395 23 Mar 15 peter 324
3395 23 Mar 15 peter 325        \since New in yat 0.13
3395 23 Mar 15 peter 326      */
3395 23 Mar 15 peter 327     char sequence_str(size_t index) const;
3395 23 Mar 15 peter 328
3395 23 Mar 15 peter 329     /**
2985 18 Feb 13 peter 330        \brief modify a base in sequence
2985 18 Feb 13 peter 331
2985 18 Feb 13 peter 332        Set i-th base in sequence to \a x, where seq is a 4-bit integer.
2985 18 Feb 13 peter 333
3843 13 Sep 19 peter 334        \see nt16_table(char)
3026 06 Apr 13 peter 335
3026 06 Apr 13 peter 336        \since New in yat 0.11
2985 18 Feb 13 peter 337      */
2985 18 Feb 13 peter 338     void sequence(size_t i, uint8_t x);
2985 18 Feb 13 peter 339
2985 18 Feb 13 peter 340     /**
3031 27 Apr 13 peter 341        \brief set sequence and quality
3031 27 Apr 13 peter 342
3031 27 Apr 13 peter 343        \since New in yat 0.11
3031 27 Apr 13 peter 344      */
3031 27 Apr 13 peter 345     void sequence(const std::string& seq, const std::vector<uint8_t>& qual);
3031 27 Apr 13 peter 346
3031 27 Apr 13 peter 347     /**
2883 03 Dec 12 peter 348        \see core().l_qseq
2883 03 Dec 12 peter 349      */
2883 03 Dec 12 peter 350     uint32_t sequence_length(void) const;
2883 03 Dec 12 peter 351
2883 03 Dec 12 peter 352     /**
2886 05 Dec 12 peter 353        Exchanging this read with \a other.
3168 21 Jan 14 peter 354
3168 21 Jan 14 peter 355        \see swap(BamRead&, BamRead&)
2886 05 Dec 12 peter 356      */
2886 05 Dec 12 peter 357     void swap(BamRead& other);
2886 05 Dec 12 peter 358
3295 25 Jul 14 peter 359     /**
3733 11 Apr 18 peter 360        \return \c true if mapped to to positive strand, i.e.,
3733 11 Apr 18 peter 361        \c BAM_FREVERSE is not set.
3733 11 Apr 18 peter 362
3733 11 Apr 18 peter 363        \since new in yat 0.16
3733 11 Apr 18 peter 364     */
3733 11 Apr 18 peter 365     bool strand(void) const;
3733 11 Apr 18 peter 366
3733 11 Apr 18 peter 367     /**
3733 11 Apr 18 peter 368        \return \c true if mate is mapped to to positive strand, i.e.,
3733 11 Apr 18 peter 369        \c BAM_FMREVERSE is not set.
3733 11 Apr 18 peter 370
3733 11 Apr 18 peter 371        \since new in yat 0.16
3733 11 Apr 18 peter 372     */
3733 11 Apr 18 peter 373     bool mstrand(void) const;
3733 11 Apr 18 peter 374
3733 11 Apr 18 peter 375     /**
3295 25 Jul 14 peter 376        \brief chromosome ID
3295 25 Jul 14 peter 377      */
3295 25 Jul 14 peter 378     int32_t tid(void) const;
3295 25 Jul 14 peter 379
2883 03 Dec 12 peter 380   private:
3031 27 Apr 13 peter 381     // ensure capacity of data pointer is (at least) n
3630 14 Mar 17 peter 382     void reserve(uint32_t n);
3031 27 Apr 13 peter 383
3889 25 Mar 20 peter 384     /**
3889 25 Mar 20 peter 385        \param k which byte in underlying sequence container to set
3889 25 Mar 20 peter 386        \param x 4-bit integer
3889 25 Mar 20 peter 387        \param y 4-bit integer
3889 25 Mar 20 peter 388
3889 25 Mar 20 peter 389        equivalent to calling
3889 25 Mar 20 peter 390        sequence(2*k, x);
3889 25 Mar 20 peter 391        sequence(2*k+1, y);
3889 25 Mar 20 peter 392     */
3889 25 Mar 20 peter 393     void sequence(size_t k, uint8_t x, uint8_t y);
3889 25 Mar 20 peter 394
2883 03 Dec 12 peter 395     bam1_t* bam_;
2883 03 Dec 12 peter 396
2883 03 Dec 12 peter 397     friend class InBamFile;
2883 03 Dec 12 peter 398     friend class OutBamFile;
2883 03 Dec 12 peter 399     friend class BamReadIterator;
3351 21 Nov 14 peter 400     // access data length, current length of data
3351 21 Nov 14 peter 401     int& data_size(void);
3351 21 Nov 14 peter 402     const int& data_size(void) const;
3630 14 Mar 17 peter 403     uint32_t data_capacity(void) const;
2883 03 Dec 12 peter 404   };
2883 03 Dec 12 peter 405
2883 03 Dec 12 peter 406   /**
3398 25 Mar 15 peter 407      Convert nucleotide character to 4-bit encoding.  The result is
3398 25 Mar 15 peter 408      encoded as 1/2/4/8 for A/C/G/T or combinations of these bits.
3398 25 Mar 15 peter 409
3398 25 Mar 15 peter 410      This is the same functionality that is provided in HTSLIB's
4169 25 Mar 22 peter 411      global array \c seq_nt16_table.
3398 25 Mar 15 peter 412
3398 25 Mar 15 peter 413      This is the inverse of function nt16_str.
3398 25 Mar 15 peter 414
3398 25 Mar 15 peter 415      \since New in yat 0.13
3398 25 Mar 15 peter 416   */
3398 25 Mar 15 peter 417   const unsigned char nt16_table(char);
3398 25 Mar 15 peter 418
3398 25 Mar 15 peter 419   /**
3398 25 Mar 15 peter 420      Convert a 4-bit encoded nucleotide to a character.
3398 25 Mar 15 peter 421
3398 25 Mar 15 peter 422      This is the same functionality that is provided in HTSLIB's
4169 25 Mar 22 peter 423      global array seq_nt16_str.
3398 25 Mar 15 peter 424
3398 25 Mar 15 peter 425      This is the inverse of function nt16_table.
3398 25 Mar 15 peter 426
3398 25 Mar 15 peter 427      \since New in yat 0.13
3398 25 Mar 15 peter 428   */
3398 25 Mar 15 peter 429   const char nt16_str(uint8_t);
3398 25 Mar 15 peter 430
3398 25 Mar 15 peter 431   /**
2886 05 Dec 12 peter 432      Swap specialization for BamRead that is faster than generic
3168 21 Jan 14 peter 433      std::swap as it just swap a pair of pointers.
2886 05 Dec 12 peter 434
2886 05 Dec 12 peter 435      \since New in yat 0.10
2986 18 Feb 13 peter 436
2986 18 Feb 13 peter 437      \relates BamRead
2886 05 Dec 12 peter 438    */
2886 05 Dec 12 peter 439   void swap(BamRead& lhs, BamRead& rhs);
2886 05 Dec 12 peter 440
2886 05 Dec 12 peter 441   /**
2883 03 Dec 12 peter 442      \return \c true if read is soft clipped, either left_soft_clipped
2883 03 Dec 12 peter 443      or right_soft_clipped.
2884 04 Dec 12 peter 444
2884 04 Dec 12 peter 445      \since New in yat 0.10
2986 18 Feb 13 peter 446
2986 18 Feb 13 peter 447      \relates BamRead
2883 03 Dec 12 peter 448    */
2883 03 Dec 12 peter 449   bool soft_clipped(const BamRead& bam);
2883 03 Dec 12 peter 450
2883 03 Dec 12 peter 451   /**
2883 03 Dec 12 peter 452      If read is soft clipped on left side, return how many bases are
2883 03 Dec 12 peter 453      clipped, otherwise return 0.
2884 04 Dec 12 peter 454
2884 04 Dec 12 peter 455      \since New in yat 0.10
2986 18 Feb 13 peter 456
2986 18 Feb 13 peter 457      \relates BamRead
2883 03 Dec 12 peter 458    */
2883 03 Dec 12 peter 459   uint32_t left_soft_clipped(const BamRead& bam);
2883 03 Dec 12 peter 460
2883 03 Dec 12 peter 461   /**
2883 03 Dec 12 peter 462      If read is soft clipped on right side, return how many bases are
2883 03 Dec 12 peter 463      clipped, otherwise return 0.
2884 04 Dec 12 peter 464
2884 04 Dec 12 peter 465      \since New in yat 0.10
2986 18 Feb 13 peter 466
2986 18 Feb 13 peter 467      \relates BamRead
2883 03 Dec 12 peter 468    */
2883 03 Dec 12 peter 469   uint32_t right_soft_clipped(const BamRead& bam);
2883 03 Dec 12 peter 470
2883 03 Dec 12 peter 471   /**
2883 03 Dec 12 peter 472      return \c true if query names are equal
2883 03 Dec 12 peter 473
2883 03 Dec 12 peter 474      \see BamRead::name()
2884 04 Dec 12 peter 475
2884 04 Dec 12 peter 476      \since New in yat 0.10
2986 18 Feb 13 peter 477
2986 18 Feb 13 peter 478      \relates BamRead
2883 03 Dec 12 peter 479    */
2883 03 Dec 12 peter 480   bool same_query_name(const BamRead& lhs, const BamRead& rhs);
2883 03 Dec 12 peter 481
2883 03 Dec 12 peter 482   /**
3384 11 Mar 15 peter 483      return \c true if lhs query name is less than rhs query name
3384 11 Mar 15 peter 484
3384 11 Mar 15 peter 485      \see BamRead::name()
3384 11 Mar 15 peter 486
3384 11 Mar 15 peter 487      \since New in yat 0.13
3384 11 Mar 15 peter 488
3384 11 Mar 15 peter 489      \relates BamRead
3384 11 Mar 15 peter 490    */
3384 11 Mar 15 peter 491   bool less_query_name(const BamRead& lhs, const BamRead& rhs);
3384 11 Mar 15 peter 492
3384 11 Mar 15 peter 493   /**
2883 03 Dec 12 peter 494      Functor to compare two reads with respect to their leftmost
2883 03 Dec 12 peter 495      coordinate.
2883 03 Dec 12 peter 496
2883 03 Dec 12 peter 497      \see BamRead
2884 04 Dec 12 peter 498
2884 04 Dec 12 peter 499      \since New in yat 0.10
2883 03 Dec 12 peter 500    */
2883 03 Dec 12 peter 501   struct BamLessPos
2883 03 Dec 12 peter 502   {
4339 15 Apr 23 peter 503     /// first argument is a const BamRead&
4339 15 Apr 23 peter 504     typedef const BamRead& first_argument_type;
4339 15 Apr 23 peter 505     /// second argument is a const BamRead&
4339 15 Apr 23 peter 506     typedef const BamRead& second_argument_type;
4339 15 Apr 23 peter 507     /// Functor returns a \c bool
4339 15 Apr 23 peter 508     typedef bool result_type;
4339 15 Apr 23 peter 509
2883 03 Dec 12 peter 510     /**
2910 17 Dec 12 peter 511        \return true if lhs tid is less than rhs tid; or tid
2883 03 Dec 12 peter 512        are equal and lhs pos is smaller than rhs pos.
2883 03 Dec 12 peter 513
2910 17 Dec 12 peter 514        \see BamRead::tid() and BamRead::pos()
2883 03 Dec 12 peter 515      */
2883 03 Dec 12 peter 516     bool operator()(const BamRead& lhs, const BamRead& rhs) const;
2883 03 Dec 12 peter 517   };
2883 03 Dec 12 peter 518
2883 03 Dec 12 peter 519
2883 03 Dec 12 peter 520   /**
2883 03 Dec 12 peter 521      Functor to compare two reads with respect to their rightmost
2883 03 Dec 12 peter 522      coordinate.
2883 03 Dec 12 peter 523
2883 03 Dec 12 peter 524      \see BamRead
2884 04 Dec 12 peter 525
2884 04 Dec 12 peter 526      \since New in yat 0.10
2883 03 Dec 12 peter 527    */
2883 03 Dec 12 peter 528   struct BamLessEnd
2883 03 Dec 12 peter 529   {
4339 15 Apr 23 peter 530     /// first argument is a const BamRead&
4339 15 Apr 23 peter 531     typedef const BamRead& first_argument_type;
4339 15 Apr 23 peter 532     /// second argument is a const BamRead&
4339 15 Apr 23 peter 533     typedef const BamRead& second_argument_type;
4339 15 Apr 23 peter 534     /// Functor returns a \c bool
4339 15 Apr 23 peter 535     typedef bool result_type;
4339 15 Apr 23 peter 536
2883 03 Dec 12 peter 537     /**
2910 17 Dec 12 peter 538        \return true if lhs tid is less than rhs tid; or tid
2883 03 Dec 12 peter 539        are equal and lhs end is smaller than rhs end.
2883 03 Dec 12 peter 540
2910 17 Dec 12 peter 541        \see BamRead::tid() and BamRead::end()
2883 03 Dec 12 peter 542      */
2883 03 Dec 12 peter 543     bool operator()(const BamRead& lhs, const BamRead& rhs) const;
2883 03 Dec 12 peter 544   };
2883 03 Dec 12 peter 545
3384 11 Mar 15 peter 546   /**
3384 11 Mar 15 peter 547      Functor to compare two reads with respect to their query names.
3384 11 Mar 15 peter 548
3384 11 Mar 15 peter 549      \see BamRead
3384 11 Mar 15 peter 550      \see less_query_name(const BamRead&, const BamRead&)
3384 11 Mar 15 peter 551
3384 11 Mar 15 peter 552      \since New in yat 0.13
3384 11 Mar 15 peter 553    */
3384 11 Mar 15 peter 554   struct BamLessName
3384 11 Mar 15 peter 555   {
4339 15 Apr 23 peter 556     /// first argument is a const BamRead&
4339 15 Apr 23 peter 557     typedef const BamRead& first_argument_type;
4339 15 Apr 23 peter 558     /// second argument is a const BamRead&
4339 15 Apr 23 peter 559     typedef const BamRead& second_argument_type;
4339 15 Apr 23 peter 560     /// Functor returns a bool
4339 15 Apr 23 peter 561     typedef bool result_type;
4339 15 Apr 23 peter 562
3384 11 Mar 15 peter 563     /**
3384 11 Mar 15 peter 564        \return \c true if lhs query name is less than rhs query name
3384 11 Mar 15 peter 565      */
3384 11 Mar 15 peter 566     bool operator()(const BamRead& lhs, const BamRead& rhs) const;
3384 11 Mar 15 peter 567   };
3384 11 Mar 15 peter 568
3395 23 Mar 15 peter 569
3395 23 Mar 15 peter 570   /**
3395 23 Mar 15 peter 571      Calculate alignment score based on \a bam's cigar and sequence
3395 23 Mar 15 peter 572      and the reference as in [\a ref, ref + bam.end()-bam.pos()).
3395 23 Mar 15 peter 573
3395 23 Mar 15 peter 574      For indels the score is calculated as open_gap + len * gap.  For
3395 23 Mar 15 peter 575      mismatches the score is calculated -mismatch. For matches the
3395 23 Mar 15 peter 576      score is unity.
3395 23 Mar 15 peter 577
3395 23 Mar 15 peter 578      For cigar element BAM_CMATCH the sequence and reference is used
3395 23 Mar 15 peter 579      to differentiate match from mismatch.
3395 23 Mar 15 peter 580
3517 16 Aug 16 peter 581      Type Requirments:
3517 16 Aug 16 peter 582      - \c RandomAccessIterator is a \readable_iterator
3517 16 Aug 16 peter 583      - \c RandomAccessIterator is a \random_access_traversal_iterator
3517 16 Aug 16 peter 584      - \c value_type is equality compareable with \c char
3517 16 Aug 16 peter 585
3395 23 Mar 15 peter 586      \since New in yat 0.13
3395 23 Mar 15 peter 587    */
3395 23 Mar 15 peter 588   template<typename RandomAccessIterator>
3395 23 Mar 15 peter 589   double alignment_score(const BamRead& bam, RandomAccessIterator ref,
3395 23 Mar 15 peter 590                          double mismatch, double gap=0.0, double open_gap=0.0)
3395 23 Mar 15 peter 591   {
3517 16 Aug 16 peter 592     BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<RandomAccessIterator>));
3517 16 Aug 16 peter 593     BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<RandomAccessIterator>));
3395 23 Mar 15 peter 594     double res = 0;
3395 23 Mar 15 peter 595     size_t q = 0;
3395 23 Mar 15 peter 596
3395 23 Mar 15 peter 597     for (uint32_t i = 0; i<bam.core().n_cigar; ++i) {
3395 23 Mar 15 peter 598       switch (bam.cigar_op(i)) {
3395 23 Mar 15 peter 599       case (BAM_CMATCH): {
3395 23 Mar 15 peter 600         RandomAccessIterator end = ref + bam.cigar_oplen(i);
3395 23 Mar 15 peter 601         for (; ref!=end; ++ref) {
3395 23 Mar 15 peter 602           if (bam.sequence_str(q) == *ref)
3395 23 Mar 15 peter 603             res += 1.0;
3395 23 Mar 15 peter 604           else
3395 23 Mar 15 peter 605             res -= mismatch;
3395 23 Mar 15 peter 606           ++q;
3395 23 Mar 15 peter 607         }
3395 23 Mar 15 peter 608         break;
3395 23 Mar 15 peter 609       }
3395 23 Mar 15 peter 610       case (BAM_CEQUAL): {
3395 23 Mar 15 peter 611         q += bam.cigar_oplen(i);
3395 23 Mar 15 peter 612         ref += bam.cigar_oplen(i);
3395 23 Mar 15 peter 613         res += bam.cigar_oplen(i);
3395 23 Mar 15 peter 614         break;
3395 23 Mar 15 peter 615       }
3395 23 Mar 15 peter 616       case (BAM_CDIFF): {
3395 23 Mar 15 peter 617         q += bam.cigar_oplen(i);
3395 23 Mar 15 peter 618         ref += bam.cigar_oplen(i);
3395 23 Mar 15 peter 619         res -= bam.cigar_oplen(i)*mismatch;
3395 23 Mar 15 peter 620         break;
3395 23 Mar 15 peter 621       }
3395 23 Mar 15 peter 622       case (BAM_CSOFT_CLIP): {
3395 23 Mar 15 peter 623         q += bam.cigar_oplen(i);
3395 23 Mar 15 peter 624         break;
3395 23 Mar 15 peter 625       }
3395 23 Mar 15 peter 626       case (BAM_CINS): {
3395 23 Mar 15 peter 627         q += bam.cigar_oplen(i);
3395 23 Mar 15 peter 628         res -= open_gap + bam.cigar_oplen(i) * gap;
3395 23 Mar 15 peter 629         break;
3395 23 Mar 15 peter 630       }
3395 23 Mar 15 peter 631       case (BAM_CDEL): {
3395 23 Mar 15 peter 632         ref += bam.cigar_oplen(i);
3395 23 Mar 15 peter 633         res -= open_gap + bam.cigar_oplen(i) * gap;
3395 23 Mar 15 peter 634         break;
3395 23 Mar 15 peter 635       }
3395 23 Mar 15 peter 636       } // end switch
3395 23 Mar 15 peter 637     }
3395 23 Mar 15 peter 638     return res;
3395 23 Mar 15 peter 639   }
3395 23 Mar 15 peter 640
2883 03 Dec 12 peter 641 }}}
2883 03 Dec 12 peter 642 #endif

Generated on Sat Oct 14 02:22:42 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)