sam_variant_input_stream.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_FORMAT_SAM_VARIANT_INPUT_STREAM_H_
#define GENESIS_POPULATION_FORMAT_SAM_VARIANT_INPUT_STREAM_H_
 
/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2024 Lucas Czech
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
    Contact:
    Lucas Czech <lucas.czech@sund.ku.dk>
    University of Copenhagen, Globe Institute, Section for GeoGenetics
    Oster Voldgade 5-7, 1350 Copenhagen K, Denmark
*/
 
#ifdef GENESIS_HTSLIB
 
#include "genesis/population/sample_counts.hpp"
#include "genesis/population/function/genome_locus.hpp"
#include "genesis/population/genome_locus_set.hpp"
#include "genesis/population/genome_locus.hpp"
#include "genesis/population/variant.hpp"
#include "genesis/population/variant.hpp"
#include "genesis/utils/core/fs.hpp"
 
#include <cstdint>
#include <memory>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <vector>
 
// Forward declarations of htslib structs, so that we do not have to include their headers
// and completely spam our namespace here.
extern "C" {
 
    struct htsFile;
    struct hts_itr_t;
    struct sam_hdr_t;
    struct bam1_t;
 
    struct bam_plp_s;
    typedef struct bam_plp_s *bam_plp_t;
 
    struct bam_pileup1_t;
    // typedef struct bam_pileup1_t *bam_pileup1_t;
 
}
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     SAM/BAM/CRAM File Input Stream
// =================================================================================================
 
class SamVariantInputStream
{
public:
 
    // -------------------------------------------------------------------------
    //     Member Types
    // -------------------------------------------------------------------------
 
    using self_type         = SamVariantInputStream;
    using value_type        = Variant;
    using pointer           = value_type const*;
    using reference         = value_type const&;
    using difference_type   = std::ptrdiff_t;
    using iterator_category = std::input_iterator_tag;
 
    // ======================================================================================
    //      Sam File Handle
    // ======================================================================================
 
private:
 
    struct SamFileHandle;
 
    // ======================================================================================
    //      Internal Iterator
    // ======================================================================================
 
public:
 
    class Iterator
    {
    public:
 
        // -------------------------------------------------------------------------
        //     Constructors and Rule of Five
        // -------------------------------------------------------------------------
 
        using self_type         = SamVariantInputStream::Iterator;
        using value_type        = Variant;
        using pointer           = value_type const*;
        using reference         = value_type const&;
        using iterator_category = std::input_iterator_tag;
 
    private:
 
        Iterator() = default;
        Iterator( SamVariantInputStream const* parent );
 
    public:
 
        ~Iterator() = default;
 
        Iterator( self_type const& ) = default;
        Iterator( self_type&& )      = default;
 
        Iterator& operator= ( self_type const& ) = default;
        Iterator& operator= ( self_type&& )      = default;
 
        friend SamVariantInputStream;
        friend SamFileHandle;
 
        // -------------------------------------------------------------------------
        //     Accessors
        // -------------------------------------------------------------------------
 
        value_type const * operator->() const
        {
            return &current_variant_;
        }
 
        value_type * operator->()
        {
            return &current_variant_;
        }
 
        value_type const & operator*() const
        {
            return current_variant_;
        }
 
        value_type & operator*()
        {
            return current_variant_;
        }
 
        // -------------------------------------------------------------------------
        //     Iteration
        // -------------------------------------------------------------------------
 
        self_type& operator ++ ()
        {
            increment_();
            return *this;
        }
 
        // self_type operator ++(int)
        // {
        //     auto cpy = *this;
        //     increment_();
        //     return cpy;
        // }
 
        bool operator==( self_type const& it ) const
        {
            return parent_ == it.parent_;
        }
 
        bool operator!=( self_type const& it ) const
        {
            return !(*this == it);
        }
 
        // -------------------------------------------------------------------------
        //     Data Access
        // -------------------------------------------------------------------------
 
        std::vector<std::string> rg_tags( bool all_header_tags = false ) const;
 
        size_t sample_size() const;
 
        // -------------------------------------------------------------------------
        //     Internal Members
        // -------------------------------------------------------------------------
 
    private:
 
        void increment_();
 
        void process_base_( bam_pileup1_t const* p );
 
        size_t get_sample_index_( bam_pileup1_t const* p ) const;
 
    private:
 
        // Parent.
        SamVariantInputStream const* parent_ = nullptr;
 
        // htslib specific file handling pointers during iteration.
        // We put this in a shared ptr so that the iterator can be copied,
        // but the htslib data gets only freed once all instances are done.
        // This also allows PIMPL to avoid including the htslib headers here.
        std::shared_ptr<SamFileHandle> handle_;
 
        // Current variant object, keeping the base tally of the current locus.
        Variant current_variant_;
 
    };
 
    // ======================================================================================
    //      Main Class
    // ======================================================================================
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    SamVariantInputStream()
        : SamVariantInputStream( std::string{} )
    {}
 
    explicit SamVariantInputStream(
        std::string const& input_file
    )
        : SamVariantInputStream( input_file, std::unordered_set<std::string>{}, false )
    {}
 
    SamVariantInputStream(
        std::string const& input_file,
        std::unordered_set<std::string> const& rg_tag_filter,
        bool inverse_rg_tag_filter = false
    );
 
    ~SamVariantInputStream() = default;
 
    SamVariantInputStream( SamVariantInputStream const& ) = default;
    SamVariantInputStream( SamVariantInputStream&& )      = default;
 
    SamVariantInputStream& operator= ( SamVariantInputStream const& ) = default;
    SamVariantInputStream& operator= ( SamVariantInputStream&& )      = default;
 
    // -------------------------------------------------------------------------
    //     Iteration
    // -------------------------------------------------------------------------
 
    Iterator begin() const
    {
        // Better error messaging that what htslib would give us if the file did not exist.
        // We check here, as input_file() allows to change the file after construction,
        // so we only do the check once we know that we are good to go.
        std::string err_str;
        if( ! utils::file_is_readable( input_file_, err_str )) {
            throw std::runtime_error(
                "Cannot open input sam/bam/cram file '" + input_file_ + "': " + err_str
            );
        }
 
        return Iterator( this );
    }
 
    Iterator end() const
    {
        return Iterator();
    }
 
    // -------------------------------------------------------------------------
    //     Basic Input Settings
    // -------------------------------------------------------------------------
 
    std::string const& input_file() const
    {
        return input_file_;
    }
 
    self_type& input_file( std::string const& value )
    {
        input_file_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Flag Settings
    // -------------------------------------------------------------------------
 
    uint32_t flags_include_all() const
    {
        return flags_include_all_;
    }
 
    self_type& flags_include_all( uint32_t value )
    {
        flags_include_all_ = value;
        return *this;
    }
 
    uint32_t flags_include_any() const
    {
        return flags_include_any_;
    }
 
    self_type& flags_include_any( uint32_t value )
    {
        flags_include_any_ = value;
        return *this;
    }
 
    uint32_t flags_exclude_all() const
    {
        return flags_exclude_all_;
    }
 
    self_type& flags_exclude_all( uint32_t value )
    {
 
        flags_exclude_all_ = value;
        return *this;
    }
 
    uint32_t flags_exclude_any() const
    {
        return flags_exclude_any_;
    }
 
    self_type& flags_exclude_any( uint32_t value )
    {
        flags_exclude_any_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Region Filter
    // -------------------------------------------------------------------------
 
    self_type& region_filter( std::shared_ptr<GenomeLocusSet> region_filter )
    {
        region_filter_ = region_filter;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Quality Settings
    // -------------------------------------------------------------------------
 
    int min_map_qual() const
    {
        return min_map_qual_;
    }
 
    self_type& min_map_qual( int value )
    {
        min_map_qual_ = value;
        return *this;
    }
 
    int min_base_qual() const
    {
        return min_base_qual_;
    }
 
    self_type& min_base_qual( int value )
    {
        min_base_qual_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Depth Settings
    // -------------------------------------------------------------------------
 
    int min_depth() const
    {
        return min_depth_;
    }
 
    self_type& min_depth( int value )
    {
        min_depth_ = value;
        return *this;
    }
 
    int max_depth() const
    {
        return max_depth_;
    }
 
    self_type& max_depth( int value )
    {
        max_depth_ = value;
        return *this;
    }
 
    int max_accumulation_depth() const
    {
        return max_acc_depth_;
    }
 
    self_type& max_accumulation_depth( int value )
    {
        max_acc_depth_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Read Group Settings
    // -------------------------------------------------------------------------
 
    bool split_by_rg() const
    {
        return split_by_rg_;
    }
 
    self_type& split_by_rg( bool value )
    {
        split_by_rg_ = value;
        return *this;
    }
 
    bool with_unaccounted_rg() const
    {
        return with_unaccounted_rg_;
    }
 
    self_type& with_unaccounted_rg( bool value )
    {
        with_unaccounted_rg_ = value;
        return *this;
    }
 
    std::unordered_set<std::string> const& rg_tag_filter() const
    {
        return rg_tag_filter_;
    }
 
    self_type& rg_tag_filter( std::unordered_set<std::string> const& value )
    {
        rg_tag_filter_ = value;
        return *this;
    }
 
    bool inverse_rg_tag_filter() const
    {
        return inverse_rg_tag_filter_;
    }
 
    self_type& inverse_rg_tag_filter( bool value )
    {
        inverse_rg_tag_filter_ = value;
        return *this;
    }
 
    // bool skip_empty_positions() const
    // {
    //     return skip_empty_positions_;
    // }
    //
    // self_type& skip_empty_positions( bool value )
    // {
    //     skip_empty_positions_ = value;
    //     return *this;
    // }
 
    // -------------------------------------------------------------------------
    //     Data Members
    // -------------------------------------------------------------------------
 
private:
 
    // Input data
    std::string input_file_;
 
    // Settings
 
    // Read filtering flags, as used by htslib and samtools
    uint32_t flags_include_all_ = 0;
    uint32_t flags_include_any_ = 0;
    uint32_t flags_exclude_all_ = 0;
    uint32_t flags_exclude_any_ = 0;
 
    // Filter for which positions to consider or skip
    std::shared_ptr<GenomeLocusSet> region_filter_;
 
    // Minimum mapping and base qualities
    int min_map_qual_ = 0;
    int min_base_qual_ = 0;
 
    // Read depth / coverage
    int min_depth_ = 0;
    int max_depth_ = 0;
    int max_acc_depth_ = 0;
 
    bool split_by_rg_ = false;
    bool with_unaccounted_rg_ = false;
    std::unordered_set<std::string> rg_tag_filter_;
    bool inverse_rg_tag_filter_ = false;
 
    // Unused / need to investiage if needed
    // int max_indel_depth_ = 250;
    // bool skip_empty_positions_ = true;
    // bool detect_overlapping_read_pairs_ = false;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // htslib guard
#endif // include guard