simple_pileup_reader.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_FORMAT_SIMPLE_PILEUP_READER_H_
#define GENESIS_POPULATION_FORMAT_SIMPLE_PILEUP_READER_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
*/
 
#include "genesis/population/variant.hpp"
#include "genesis/sequence/functions/quality.hpp"
#include "genesis/utils/io/input_source.hpp"
#include "genesis/utils/io/input_stream.hpp"
 
#include <array>
#include <string>
#include <vector>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Simple (m)pileup Reader
// =================================================================================================
 
class SimplePileupReader
{
public:
 
    // -------------------------------------------------------------------------
    //     Typedefs and Enums
    // -------------------------------------------------------------------------
 
    struct Sample
    {
        size_t read_depth = 0;
 
        std::string read_bases;
 
        std::vector<unsigned char> phred_scores;
 
        char ancestral_base = '\0';
    };
 
    struct Record
    {
        std::string         chromosome = "";
        size_t              position = 0;
        char                reference_base = 'N';
        std::vector<Sample> samples;
    };
 
    using self_type = SimplePileupReader;
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    SimplePileupReader() = default;
    ~SimplePileupReader() = default;
 
    SimplePileupReader( self_type const& ) = default;
    SimplePileupReader( self_type&& )      = default;
 
    self_type& operator= ( self_type const& ) = default;
    self_type& operator= ( self_type&& )      = default;
 
    // ---------------------------------------------------------------------
    //     Reading Records
    // ---------------------------------------------------------------------
 
    std::vector<Record> read_records( std::shared_ptr< utils::BaseInputSource > source ) const;
 
    // /* *
    //  * @brief Read an (m)pileup file line by line, but only the samples at the given indices,
    //  * as pileup Record%s.
    //  */
    // std::vector<Record> read_records(
    //     std::shared_ptr< utils::BaseInputSource > source,
    //     std::vector<size_t> const&                sample_indices
    // ) const;
 
    std::vector<Record> read_records(
        std::shared_ptr< utils::BaseInputSource > source,
        std::vector<bool> const&                  sample_filter
    ) const;
 
    // ---------------------------------------------------------------------
    //     Reading Variants
    // ---------------------------------------------------------------------
 
    std::vector<Variant> read_variants( std::shared_ptr< utils::BaseInputSource > source ) const;
 
    // /* *
    //  * @brief Read an (m)pileup file line by line, but only the samples at the given indices,
    //  * as Variant%s.
    //  */
    // std::vector<Variant> read_variants(
    //     std::shared_ptr< utils::BaseInputSource > source,
    //     std::vector<size_t> const&                sample_indices
    // ) const;
 
    std::vector<Variant> read_variants(
        std::shared_ptr< utils::BaseInputSource > source,
        std::vector<bool> const&                  sample_filter
    ) const;
 
    // -------------------------------------------------------------------------
    //     Parsing Records
    // -------------------------------------------------------------------------
 
    bool parse_line_record(
        utils::InputStream& input_stream,
        Record&             record
    ) const;
 
    bool parse_line_record(
        utils::InputStream&      input_stream,
        Record&                  record,
        std::vector<bool> const& sample_filter
    ) const;
 
    // -------------------------------------------------------------------------
    //     Parsing Variants
    // -------------------------------------------------------------------------
 
    bool parse_line_variant(
        utils::InputStream& input_stream,
        Variant&            variant
    ) const;
 
    bool parse_line_variant(
        utils::InputStream&      input_stream,
        Variant&                 variant,
        std::vector<bool> const& sample_filter
    ) const;
 
    // -------------------------------------------------------------------------
    //     General Settings
    // -------------------------------------------------------------------------
 
    bool strict_bases() const
    {
        return strict_bases_;
    }
 
    self_type& strict_bases( bool value )
    {
        strict_bases_ = value;
        return *this;
    }
 
    bool with_quality_string() const
    {
        return with_quality_string_;
    }
 
    self_type& with_quality_string( bool value )
    {
        with_quality_string_ = value;
        return *this;
    }
 
    sequence::QualityEncoding quality_encoding() const
    {
        return quality_encoding_;
    }
 
    self_type& quality_encoding( sequence::QualityEncoding value )
    {
        quality_encoding_ = value;
        return *this;
    }
 
    std::array<size_t, 128> const& quality_code_counts() const
    {
        return quality_code_counts_;
    }
 
    bool with_ancestral_base() const
    {
        return with_ancestral_base_;
    }
 
    self_type& with_ancestral_base( bool value )
    {
        with_ancestral_base_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Variant Settings
    // -------------------------------------------------------------------------
 
    size_t min_base_quality() const
    {
        return min_base_quality_;
    }
 
    self_type& min_base_quality( size_t value )
    {
        min_base_quality_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Internal Members
    // -------------------------------------------------------------------------
 
private:
 
    void reset_status_( Variant& variant ) const;
 
    template<class T>
    bool parse_line_(
        utils::InputStream&      input_stream,
        T&                       target,
        std::vector<bool> const& sample_filter,
        bool                     use_sample_filter
    ) const;
 
    template<class T, class S>
    void process_sample_(
        utils::InputStream& input_stream,
        T const&            target,
        S&                  sample
    ) const;
 
    template<class T>
    void set_target_alternative_base_(
        T& target
    ) const;
 
    template<class S>
    void set_sample_read_depth_(
        size_t read_depth,
        S& sample
    ) const;
 
    bool process_sample_read_bases_buffer_(
        utils::InputStream& input_stream,
        char reference_base
    ) const;
 
    void process_sample_read_bases_stream_(
        utils::InputStream& input_stream,
        char reference_base
    ) const;
 
    template<class S>
    void set_sample_read_bases_(
        std::string const& read_bases,
        S& sample
    ) const;
 
    template<class S>
    void process_quality_string_(
        utils::InputStream& input_stream,
        S& sample
    ) const;
 
    void tally_base_(
        utils::InputStream& input_stream,
        SampleCounts& sample,
        char b
    ) const;
 
    template<class S>
    void process_ancestral_base_(
        utils::InputStream& input_stream,
        S& sample
    ) const;
 
    void skip_sample_(
        utils::InputStream& input_stream
    ) const;
 
    void next_field_(
        utils::InputStream& input_stream
    ) const;
 
    // -------------------------------------------------------------------------
    //     Data Members
    // -------------------------------------------------------------------------
 
private:
 
    // If set, we expect bases to be ACGTN. If not set, we will fix any that are not to N.
    bool strict_bases_ = false;
 
    // Set whether the file contains the base quality score column, and if so, how its encoded
    // (we default to Sanger with offset 33), and if we want to skip low quality bases.
    bool with_quality_string_ = true;
    sequence::QualityEncoding quality_encoding_ = sequence::QualityEncoding::kSanger;
    size_t min_base_quality_ = 0;
 
    // We also keep track of the base codes found, to check that we have the right encoding.
    mutable std::array<size_t, 128> quality_code_counts_;
 
    // Set whether the last part of the sample line contains the base of the ancestral allele.
    bool with_ancestral_base_ = false;
 
    // Internal buffer to read bases into. Used for speedup to avoid reallocations.
    mutable std::string base_buffer_;
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard