yat/omic/BamFile.h

Code
Comments
Other
Rev Date Author Line
2883 03 Dec 12 peter 1 #ifndef theplu_yat_omic_bam_file
2883 03 Dec 12 peter 2 #define theplu_yat_omic_bam_file
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 /*
3999 08 Oct 20 peter 7   Copyright (C) 2012, 2013, 2014, 2016, 2017, 2018, 2020 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
2883 03 Dec 12 peter 25 #include "BamHeader.h"
2883 03 Dec 12 peter 26 #include "BamRead.h"
2883 03 Dec 12 peter 27
3883 24 Mar 20 peter 28 #include "yat/utility/deprecate.h"
2897 11 Dec 12 peter 29 #include "yat/utility/Exception.h"
3674 01 Aug 17 peter 30 #include "yat/utility/FileUtil.h"
3048 12 Jun 13 peter 31 #include "yat/utility/yat_assert.h"
2897 11 Dec 12 peter 32
3883 24 Mar 20 peter 33 #include <htslib/sam.h>
2883 03 Dec 12 peter 34
2883 03 Dec 12 peter 35 #include <boost/utility.hpp>
2883 03 Dec 12 peter 36
2883 03 Dec 12 peter 37 #include <cstdio>
2883 03 Dec 12 peter 38 #include <sstream>
2883 03 Dec 12 peter 39 #include <stdexcept>
2883 03 Dec 12 peter 40 #include <string>
2883 03 Dec 12 peter 41
2883 03 Dec 12 peter 42 namespace theplu {
2883 03 Dec 12 peter 43 namespace yat {
2883 03 Dec 12 peter 44 namespace omic {
2883 03 Dec 12 peter 45
2883 03 Dec 12 peter 46   /**
2883 03 Dec 12 peter 47      Base class for bam files
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   template<typename Derived>
2883 03 Dec 12 peter 52   class BamFile : boost::noncopyable
2883 03 Dec 12 peter 53   {
2883 03 Dec 12 peter 54     typedef Derived derived_type;
2883 03 Dec 12 peter 55   public:
2883 03 Dec 12 peter 56     /**
2883 03 Dec 12 peter 57        Default Constructor
2883 03 Dec 12 peter 58      */
2883 03 Dec 12 peter 59     BamFile(void);
2883 03 Dec 12 peter 60
2883 03 Dec 12 peter 61     /**
2883 03 Dec 12 peter 62        \brief Destructor
2883 03 Dec 12 peter 63
2883 03 Dec 12 peter 64        Closes file if it is open.
2883 03 Dec 12 peter 65      */
2883 03 Dec 12 peter 66     virtual ~BamFile(void);
2883 03 Dec 12 peter 67
2883 03 Dec 12 peter 68     /**
2883 03 Dec 12 peter 69        \brief close file
2883 03 Dec 12 peter 70      */
2883 03 Dec 12 peter 71     void close(void);
2883 03 Dec 12 peter 72
2883 03 Dec 12 peter 73     /**
3674 01 Aug 17 peter 74        \return true if file has an index file
3674 01 Aug 17 peter 75
3674 01 Aug 17 peter 76        \since New in yat 0.15
3674 01 Aug 17 peter 77      */
3674 01 Aug 17 peter 78     bool have_index(void) const;
3674 01 Aug 17 peter 79
3674 01 Aug 17 peter 80     /**
2883 03 Dec 12 peter 81        \return \c true iff open
2883 03 Dec 12 peter 82      */
2883 03 Dec 12 peter 83     bool is_open(void) const;
2883 03 Dec 12 peter 84   protected:
2883 03 Dec 12 peter 85     /**
3674 01 Aug 17 peter 86        \since New in yat 0.15
3674 01 Aug 17 peter 87      */
3674 01 Aug 17 peter 88     void build_index_base(void) const;
3674 01 Aug 17 peter 89
3674 01 Aug 17 peter 90     /**
2883 03 Dec 12 peter 91        open a bam file named \a fn with mode \a mode
2883 03 Dec 12 peter 92
3883 24 Mar 20 peter 93        \see sam_open in htslib
3357 23 Nov 14 peter 94
3357 23 Nov 14 peter 95        \note \a aux is ignored when compiling against htslib
3883 24 Mar 20 peter 96
3883 24 Mar 20 peter 97        \deprecated Provided for backward compatibility with the yat
3883 24 Mar 20 peter 98        0.17 API. Use open_base(2) instead.
2883 03 Dec 12 peter 99      */
2883 03 Dec 12 peter 100     void open_base(const std::string& fn, const std::string& mode,
3883 24 Mar 20 peter 101                    const void* aux) YAT_DEPRECATE;
3353 22 Nov 14 peter 102
3363 25 Nov 14 peter 103     /**
3883 24 Mar 20 peter 104        open a bam file named \a fn with mode \a mode
3363 25 Nov 14 peter 105
3883 24 Mar 20 peter 106        \see sam_open in htslib
3883 24 Mar 20 peter 107
3883 24 Mar 20 peter 108        \since New in yat 0.18
3363 25 Nov 14 peter 109     */
3883 24 Mar 20 peter 110     void open_base(const std::string& fn, const std::string& mode);
3353 22 Nov 14 peter 111
2883 03 Dec 12 peter 112     /**
2883 03 Dec 12 peter 113        bam file handler
2883 03 Dec 12 peter 114      */
3353 22 Nov 14 peter 115     samFile* sf_;
2883 03 Dec 12 peter 116
2883 03 Dec 12 peter 117     /**
2883 03 Dec 12 peter 118        \brief filename of bam file
2883 03 Dec 12 peter 119
2883 03 Dec 12 peter 120        The filename is set in open_base().
2883 03 Dec 12 peter 121      */
2883 03 Dec 12 peter 122     const std::string& filename(void) const { return filename_; }
3674 01 Aug 17 peter 123
2883 03 Dec 12 peter 124   private:
2883 03 Dec 12 peter 125     std::string filename_;
2883 03 Dec 12 peter 126   };
2883 03 Dec 12 peter 127
2883 03 Dec 12 peter 128
2883 03 Dec 12 peter 129   /**
2883 03 Dec 12 peter 130      This class supports reading from a bam file.
2884 04 Dec 12 peter 131
2884 04 Dec 12 peter 132      \since New in yat 0.10
2883 03 Dec 12 peter 133    */
2883 03 Dec 12 peter 134   class InBamFile : public BamFile<InBamFile>
2883 03 Dec 12 peter 135   {
2883 03 Dec 12 peter 136     typedef BamFile<InBamFile> super_t;
2883 03 Dec 12 peter 137   public:
2883 03 Dec 12 peter 138     /**
3353 22 Nov 14 peter 139        Struct describing bam index.
3353 22 Nov 14 peter 140
3883 24 Mar 20 peter 141        \see hts_idx_t in htslib
3353 22 Nov 14 peter 142
3363 25 Nov 14 peter 143        \since New in yat 0.13
3363 25 Nov 14 peter 144
3363 25 Nov 14 peter 145        \see function index(void) const
3353 22 Nov 14 peter 146      */
3353 22 Nov 14 peter 147     typedef hts_idx_t index_type;
3353 22 Nov 14 peter 148
3353 22 Nov 14 peter 149     /**
2883 03 Dec 12 peter 150        \brief Default constructor
2883 03 Dec 12 peter 151      */
2883 03 Dec 12 peter 152     InBamFile(void);
2883 03 Dec 12 peter 153
2883 03 Dec 12 peter 154     /**
2883 03 Dec 12 peter 155        Create an input file
2883 03 Dec 12 peter 156
2883 03 Dec 12 peter 157        \param fn string specifying the filename
2883 03 Dec 12 peter 158
2883 03 Dec 12 peter 159        \see open(const std::string&)
2883 03 Dec 12 peter 160      */
2883 03 Dec 12 peter 161     explicit InBamFile(const std::string& fn);
2883 03 Dec 12 peter 162
2883 03 Dec 12 peter 163     /**
2883 03 Dec 12 peter 164        \brief destructor
2883 03 Dec 12 peter 165      */
2883 03 Dec 12 peter 166     virtual ~InBamFile(void);
2883 03 Dec 12 peter 167
2883 03 Dec 12 peter 168     /**
3674 01 Aug 17 peter 169        \since New in yat 0.15
3674 01 Aug 17 peter 170      */
3674 01 Aug 17 peter 171     void build_index(void) const;
3674 01 Aug 17 peter 172
3674 01 Aug 17 peter 173     /**
2883 03 Dec 12 peter 174        \return header
2883 03 Dec 12 peter 175      */
2883 03 Dec 12 peter 176     const BamHeader& header(void) const;
2883 03 Dec 12 peter 177
2883 03 Dec 12 peter 178     /**
2883 03 Dec 12 peter 179        \return index associated with BamFile
2883 03 Dec 12 peter 180
2883 03 Dec 12 peter 181        First time this function is called an index file is loaded from
2883 03 Dec 12 peter 182        disk. If bam file is named 'foo.bam', the index file should be
2883 03 Dec 12 peter 183        named 'foo.bam.bai'. If no such file exists or this bam file
3077 11 Sep 13 peter 184        reads from stdin, this function throws.
2883 03 Dec 12 peter 185      */
3353 22 Nov 14 peter 186     const index_type* index(void) const;
2883 03 Dec 12 peter 187
2883 03 Dec 12 peter 188     /**
3464 09 Feb 16 peter 189        \return number of mapped reads in segment \a tid
3464 09 Feb 16 peter 190      */
3464 09 Feb 16 peter 191     uint64_t n_mapped(int tid) const;
3464 09 Feb 16 peter 192
3464 09 Feb 16 peter 193     /**
3464 09 Feb 16 peter 194        Typically unmapped reads have same coordinate as its mate.
3464 09 Feb 16 peter 195
3464 09 Feb 16 peter 196        \return number of unmapped reads with segment \a tid
3464 09 Feb 16 peter 197      */
3464 09 Feb 16 peter 198     uint64_t n_unmapped(int tid) const;
3464 09 Feb 16 peter 199
3464 09 Feb 16 peter 200     /**
3464 09 Feb 16 peter 201        \return number of reads with no coordinate
3464 09 Feb 16 peter 202      */
3464 09 Feb 16 peter 203     uint64_t n_no_coordinate(void) const;
3464 09 Feb 16 peter 204
3464 09 Feb 16 peter 205     /**
2883 03 Dec 12 peter 206        \brief Open an input bam file.
2883 03 Dec 12 peter 207
3077 11 Sep 13 peter 208        Open a bam file named \a fn. If \a fn is "-", \c stdin is used.
3077 11 Sep 13 peter 209
2883 03 Dec 12 peter 210        \param fn string specifying the filename
2883 03 Dec 12 peter 211      */
2883 03 Dec 12 peter 212     void open(const std::string& fn);
2883 03 Dec 12 peter 213
2883 03 Dec 12 peter 214     /**
2883 03 Dec 12 peter 215        \brief read the next BamRead
2883 03 Dec 12 peter 216
2883 03 Dec 12 peter 217        \return true on success
2883 03 Dec 12 peter 218      */
2883 03 Dec 12 peter 219     bool read(BamRead& read);
2883 03 Dec 12 peter 220
2883 03 Dec 12 peter 221     /**
2883 03 Dec 12 peter 222        \brief read the next BamRead
2883 03 Dec 12 peter 223
3363 25 Nov 14 peter 224        This function is used in BamReadIterator. Typically it is more
3363 25 Nov 14 peter 225        convenient to read bams via a BamReadIterator than through this
3363 25 Nov 14 peter 226        function directly.
3363 25 Nov 14 peter 227
2883 03 Dec 12 peter 228        \return true on success
2883 03 Dec 12 peter 229      */
3353 22 Nov 14 peter 230     bool read(BamRead& read, hts_itr_t* iter);
2883 03 Dec 12 peter 231   private:
3464 09 Feb 16 peter 232     uint64_t get_idx_stat(int tid, bool return_mapped) const;
2883 03 Dec 12 peter 233     BamHeader header_;
3464 09 Feb 16 peter 234     // always access index_ via function index(), so index is loaded
3464 09 Feb 16 peter 235     // if needed
3353 22 Nov 14 peter 236     mutable index_type* index_;
2883 03 Dec 12 peter 237   };
2883 03 Dec 12 peter 238
2883 03 Dec 12 peter 239
2883 03 Dec 12 peter 240   /**
2883 03 Dec 12 peter 241      This class supports writing to a bam file.
2884 04 Dec 12 peter 242
2884 04 Dec 12 peter 243      \since New in yat 0.10
2883 03 Dec 12 peter 244    */
2883 03 Dec 12 peter 245   class OutBamFile : public BamFile<OutBamFile>
2883 03 Dec 12 peter 246   {
2883 03 Dec 12 peter 247     typedef BamFile<OutBamFile> super_t;
2883 03 Dec 12 peter 248   public:
2883 03 Dec 12 peter 249     /**
2883 03 Dec 12 peter 250        Create an output bam file
2883 03 Dec 12 peter 251      */
2883 03 Dec 12 peter 252     OutBamFile(void);
2883 03 Dec 12 peter 253
2883 03 Dec 12 peter 254     /**
2883 03 Dec 12 peter 255        \brief Create an output bam file.
2883 03 Dec 12 peter 256
2883 03 Dec 12 peter 257        Equivalent to default constructor followed by a call to open(2).
2883 03 Dec 12 peter 258
2883 03 Dec 12 peter 259        \see open(const std::string&, const BamHeader&)
2883 03 Dec 12 peter 260      */
2883 03 Dec 12 peter 261     OutBamFile(const std::string&, const BamHeader& header);
2883 03 Dec 12 peter 262
2883 03 Dec 12 peter 263     /**
3160 13 Jan 14 peter 264        \brief Create an output bam file.
3160 13 Jan 14 peter 265
3160 13 Jan 14 peter 266        Equivalent to default constructor followed by a call to open(3).
3160 13 Jan 14 peter 267
3160 13 Jan 14 peter 268        \see open(const std::string&, const BamHeader&, unsigned int)
3160 13 Jan 14 peter 269
3160 13 Jan 14 peter 270        \since yat 0.12
3160 13 Jan 14 peter 271      */
3160 13 Jan 14 peter 272     OutBamFile(const std::string&, const BamHeader& header,
3160 13 Jan 14 peter 273                unsigned int compression);
3160 13 Jan 14 peter 274
3160 13 Jan 14 peter 275     /**
3674 01 Aug 17 peter 276        \since New in yat 0.15
3674 01 Aug 17 peter 277
3674 01 Aug 17 peter 278        \throw if file is open (\see is_open())
3674 01 Aug 17 peter 279      */
3674 01 Aug 17 peter 280     void build_index(void) const;
3674 01 Aug 17 peter 281
3674 01 Aug 17 peter 282     /**
2883 03 Dec 12 peter 283        \brief Open an output bam file.
2883 03 Dec 12 peter 284
2883 03 Dec 12 peter 285        Opens an output bam file and writes the header contained in \a
3077 11 Sep 13 peter 286        hdr. If \a fn is "-", \c stdout is used.
2883 03 Dec 12 peter 287
2883 03 Dec 12 peter 288        \param fn string specifying the filename
2883 03 Dec 12 peter 289        \param hdr header
2883 03 Dec 12 peter 290      */
2883 03 Dec 12 peter 291     void open(const std::string& fn, const BamHeader& hdr);
2883 03 Dec 12 peter 292
2883 03 Dec 12 peter 293     /**
3160 13 Jan 14 peter 294        \brief Open an output bam file.
3160 13 Jan 14 peter 295
3160 13 Jan 14 peter 296        Opens an output bam file and writes the header contained in \a
3160 13 Jan 14 peter 297        hdr. If \a fn is "-", \c stdout is used.
3160 13 Jan 14 peter 298
3160 13 Jan 14 peter 299        \param fn string specifying the filename
3160 13 Jan 14 peter 300        \param hdr header
3160 13 Jan 14 peter 301        \param compression a number [0,9] indicating level of
3160 13 Jan 14 peter 302        compression. 9 gives highest compression and 0 indicates no
3363 25 Nov 14 peter 303        compression (suitable for piping between applications).
3160 13 Jan 14 peter 304
3160 13 Jan 14 peter 305        \since yat 0.12
3160 13 Jan 14 peter 306      */
3160 13 Jan 14 peter 307     void open(const std::string& fn, const BamHeader& hdr,
3160 13 Jan 14 peter 308               unsigned int compression);
3160 13 Jan 14 peter 309
3160 13 Jan 14 peter 310     /**
2883 03 Dec 12 peter 311        \brief write a read to output file
2983 04 Feb 13 peter 312
3166 17 Jan 14 peter 313        \throw OutBamFile::error if write failed
3166 17 Jan 14 peter 314
3166 17 Jan 14 peter 315        \since yat 0.12
2883 03 Dec 12 peter 316      */
2883 03 Dec 12 peter 317     void write(const BamRead& read);
3166 17 Jan 14 peter 318
3166 17 Jan 14 peter 319     /**
3166 17 Jan 14 peter 320        \brief Error thrown from OutBamFile::write(const BamRead&) at failure
3166 17 Jan 14 peter 321      */
3166 17 Jan 14 peter 322     class error : public utility::IO_error
3166 17 Jan 14 peter 323     {
3166 17 Jan 14 peter 324     public:
3166 17 Jan 14 peter 325       /// \brief Constructor
3166 17 Jan 14 peter 326       error(const BamRead&);
3166 17 Jan 14 peter 327       /// \brief Destructor
3166 17 Jan 14 peter 328       // has to be throw() since base class destructor is
3166 17 Jan 14 peter 329       virtual ~error(void) throw();
3166 17 Jan 14 peter 330       /**
3166 17 Jan 14 peter 331          \return bam read that failed
3166 17 Jan 14 peter 332        */
3166 17 Jan 14 peter 333       const BamRead& read(void) const;
3166 17 Jan 14 peter 334     private:
3166 17 Jan 14 peter 335       BamRead read_;
3166 17 Jan 14 peter 336     }; // end of class error
3166 17 Jan 14 peter 337
2883 03 Dec 12 peter 338   private:
3630 14 Mar 17 peter 339     void open(const std::string& fn, const std::string& mode,
3630 14 Mar 17 peter 340               const BamHeader& h);
3630 14 Mar 17 peter 341
2883 03 Dec 12 peter 342   };
2883 03 Dec 12 peter 343
2883 03 Dec 12 peter 344
2883 03 Dec 12 peter 345   // template implementations
2883 03 Dec 12 peter 346   template<class Derived>
2883 03 Dec 12 peter 347   BamFile<Derived>::BamFile(void)
2883 03 Dec 12 peter 348   : sf_(NULL) {}
2883 03 Dec 12 peter 349
2883 03 Dec 12 peter 350
2883 03 Dec 12 peter 351   template<class Derived>
2883 03 Dec 12 peter 352   BamFile<Derived>::~BamFile(void)
2883 03 Dec 12 peter 353   {
2883 03 Dec 12 peter 354     close();
2883 03 Dec 12 peter 355   }
2883 03 Dec 12 peter 356
2883 03 Dec 12 peter 357
2883 03 Dec 12 peter 358   template<class Derived>
3674 01 Aug 17 peter 359   void BamFile<Derived>::build_index_base(void) const
3674 01 Aug 17 peter 360   {
3674 01 Aug 17 peter 361     int res = bam_index_build(filename_.c_str(), 0);
3674 01 Aug 17 peter 362     if (res) {
3674 01 Aug 17 peter 363       std::ostringstream msg;
3674 01 Aug 17 peter 364       msg << "failed building index file '" << filename_ << ".bai': ";
3674 01 Aug 17 peter 365       if (res == -1)
3674 01 Aug 17 peter 366         msg << "failed to build index";
3674 01 Aug 17 peter 367       else if (res == -2)
3674 01 Aug 17 peter 368         msg << "failed to open file '" << filename_ << "'";
3674 01 Aug 17 peter 369       else if (res == -4)
3674 01 Aug 17 peter 370         msg << "failed to save file";
3674 01 Aug 17 peter 371       throw utility::runtime_error(msg.str());
3674 01 Aug 17 peter 372     }
3674 01 Aug 17 peter 373   }
3674 01 Aug 17 peter 374
3674 01 Aug 17 peter 375
3674 01 Aug 17 peter 376   template<class Derived>
2883 03 Dec 12 peter 377   void BamFile<Derived>::close(void)
2883 03 Dec 12 peter 378   {
3357 23 Nov 14 peter 379     if (sf_==NULL)
3357 23 Nov 14 peter 380       return;
3357 23 Nov 14 peter 381     if (sam_close(sf_))
3357 23 Nov 14 peter 382       throw utility::IO_error("BamFile::close() failed");
2883 03 Dec 12 peter 383     sf_ = NULL;
2883 03 Dec 12 peter 384   }
2883 03 Dec 12 peter 385
2883 03 Dec 12 peter 386
2883 03 Dec 12 peter 387   template<class Derived>
2883 03 Dec 12 peter 388   bool BamFile<Derived>::is_open(void) const
2883 03 Dec 12 peter 389   {
2883 03 Dec 12 peter 390     return sf_;
2883 03 Dec 12 peter 391   }
2883 03 Dec 12 peter 392
2883 03 Dec 12 peter 393
2883 03 Dec 12 peter 394   template<class Derived>
3674 01 Aug 17 peter 395   bool BamFile<Derived>::have_index(void) const
3674 01 Aug 17 peter 396   {
3674 01 Aug 17 peter 397     return utility::FileUtil(filename() + ".bai").exists();
3674 01 Aug 17 peter 398   }
3674 01 Aug 17 peter 399
3674 01 Aug 17 peter 400
3674 01 Aug 17 peter 401   template<class Derived>
2883 03 Dec 12 peter 402   void BamFile<Derived>::open_base(const std::string& fn,
2883 03 Dec 12 peter 403                                    const std::string& mode,
2883 03 Dec 12 peter 404                                    const void* aux)
2883 03 Dec 12 peter 405   {
3883 24 Mar 20 peter 406     YAT_ASSERT(aux==NULL); // aux is ignored in htslib mode
3883 24 Mar 20 peter 407     open_base(fn, mode);
3883 24 Mar 20 peter 408   }
3883 24 Mar 20 peter 409
3883 24 Mar 20 peter 410
3883 24 Mar 20 peter 411   template<class Derived>
3883 24 Mar 20 peter 412   void BamFile<Derived>::open_base(const std::string& fn,
3883 24 Mar 20 peter 413                                    const std::string& mode)
3883 24 Mar 20 peter 414   {
2883 03 Dec 12 peter 415     filename_ = fn;
3048 12 Jun 13 peter 416     YAT_ASSERT(!sf_);
3357 23 Nov 14 peter 417     sf_ = sam_open(fn.c_str(), mode.c_str());
2883 03 Dec 12 peter 418     if (!sf_) {
2883 03 Dec 12 peter 419       std::ostringstream ss;
3078 11 Sep 13 peter 420       ss << "failed open '" << fn << "'";
2897 11 Dec 12 peter 421       throw utility::runtime_error(ss.str());
2883 03 Dec 12 peter 422     }
2883 03 Dec 12 peter 423   }
2883 03 Dec 12 peter 424
2883 03 Dec 12 peter 425 }}}
2883 03 Dec 12 peter 426 #endif