frequency_table_input_stream.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_FORMAT_FREQUENCY_TABLE_INPUT_STREAM_H_
#define GENESIS_POPULATION_FORMAT_FREQUENCY_TABLE_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
*/
 
#include "genesis/population/variant.hpp"
 
#include "genesis/sequence/reference_genome.hpp"
#include "genesis/utils/io/input_source.hpp"
#include "genesis/utils/io/input_stream.hpp"
 
#include <functional>
#include <limits>
#include <memory>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Frequency Table Input Stream
// =================================================================================================
 
class FrequencyTableInputStream
{
public:
 
    // -------------------------------------------------------------------------
    //     Member Types
    // -------------------------------------------------------------------------
 
    using self_type         = FrequencyTableInputStream;
    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;
 
    // ======================================================================================
    //      Internal Iterator
    // ======================================================================================
 
public:
 
    class Iterator
    {
        // -------------------------------------------------------------------------
        //     Typedefs and Enums
        // -------------------------------------------------------------------------
 
    public:
 
        using self_type         = FrequencyTableInputStream::Iterator;
        using value_type        = Variant;
        using pointer           = value_type const*;
        using reference         = value_type const&;
        using iterator_category = std::input_iterator_tag;
 
    private:
 
        // When parsing the header line, we need to keep track of what information is present.
        // For this, we use information per sample, as well as about the constant fields.
 
        // Per sample information, to make sure that we can actually process a sample.
        // Do we have ref and alt counts? Do we have frequency and/or read depth?
        struct SampleInfo
        {
            size_t index = std::numeric_limits<size_t>::max();
            bool has_ref = false;
            bool has_alt = false;
            bool has_frq = false;
            bool has_cov = false;
        };
 
        // All column info, to make sure that we have all the data needed before processing.
        // Do we have the chromosome and position column? Do we have the ref and alt base column?
        // Which samples are present, and which data do they contain?
        struct HeaderInfo
        {
            // We need to make sure that we only see each data type once.
            bool has_chr = false;
            bool has_pos = false;
            bool has_ref = false;
            bool has_alt = false;
 
            // We make a list of all sample names that we have seen, noting which type
            // of field we found, to make sure that we have all relevant ones. The order of the
            // indices is the order in which the SampleCounts in the final Variant are produced.
            // Shall be accessed through get_sample_info_() to ensure correct setup of the index.
            std::unordered_map<std::string, SampleInfo> sample_infos;
        };
 
        // Store the data that is parsed per sample, before processing it into our final Variant.
        struct SampleData
        {
            bool   is_missing = false;
            size_t ref_cnt    = 0;
            size_t alt_cnt    = 0;
            size_t cov        = 0;
            double frq        = 0.0;
        };
 
        // -------------------------------------------------------------------------
        //     Constructors and Rule of Five
        // -------------------------------------------------------------------------
 
    private:
 
        Iterator() = default;
        Iterator( FrequencyTableInputStream const* parent )
            : parent_( parent )
            , input_stream_(    std::make_shared<utils::InputStream>( parent_->input_source_ ))
            , sample_data_(     std::make_shared<std::vector<SampleData>>() )
            , current_variant_( std::make_shared<Variant>() )
        {
            // Safeguard
            if( ! parent_ ) {
                return;
            }
 
            // We need a manually constructored shared ptr, as this is private,
            // and make_shared does not have access to this. Hm nope, seems to work?!
            // sample_data_ = std::shared_ptr<std::vector<SampleData>>( new std::vector<SampleData>() );
 
            // Start streaming the data
            parse_header_();
            increment_();
        }
 
    public:
 
        ~Iterator() = default;
 
        Iterator( self_type const& ) = default;
        Iterator( self_type&& )      = default;
 
        Iterator& operator= ( self_type const& ) = default;
        Iterator& operator= ( self_type&& )      = default;
 
        friend FrequencyTableInputStream;
 
        // -------------------------------------------------------------------------
        //     Accessors
        // -------------------------------------------------------------------------
 
        value_type const* operator->() const
        {
            assert( current_variant_ );
            return &*current_variant_;
        }
 
        // value_type* operator->()
        // {
        //     return &current_variant_;
        // }
 
        value_type const& operator*() const
        {
            assert( current_variant_ );
            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> sample_names() const;
 
        // -------------------------------------------------------------------------
        //     Internal Members
        // -------------------------------------------------------------------------
 
    private:
 
        // ---------------------------------------------
        //     Init and Header Parsing
        // ---------------------------------------------
 
        void parse_header_();
 
        void check_header_fields_(
            std::unordered_set<std::string> const& all_samplenames
        ) const;
 
        // ---------------------------------------------
        //     Parse Header Fields
        // ---------------------------------------------
 
        void parse_header_field_(
            std::string const& field,
            std::unordered_set<std::string>& all_samplenames
        );
 
        int evaluate_if_field_is_chr_( std::string const& field );
        int evaluate_if_field_is_pos_( std::string const& field );
        int evaluate_if_field_is_ref_( std::string const& field );
        int evaluate_if_field_is_alt_( std::string const& field );
 
        int evaluate_if_field_is_sample_ref_(
            std::string const& field,
            std::unordered_set<std::string>& all_samplenames
        );
        int evaluate_if_field_is_sample_alt_(
            std::string const& field,
            std::unordered_set<std::string>& all_samplenames
        );
        int evaluate_if_field_is_sample_frq_(
            std::string const& field,
            std::unordered_set<std::string>& all_samplenames
        );
        int evaluate_if_field_is_sample_cov_(
            std::string const& field,
            std::unordered_set<std::string>& all_samplenames
        );
 
        // ---------------------------------------------
        //     Sample Helpers
        // ---------------------------------------------
 
        SampleInfo& get_sample_info_( std::string const& samplename );
        bool is_ignored_sample_( std::string const& samplename ) const;
        static bool parse_if_missing_(
            FrequencyTableInputStream const* parent,
            genesis::utils::InputStream& input_stream
        );
 
        // ---------------------------------------------
        //     String Matching Helpers
        // ---------------------------------------------
 
        bool match_header_field_(
            std::string const& field,
            std::string const& user_string,
            std::vector<std::string> const& predefined_list
        ) const;
 
        bool match_header_sample_(
            std::string const& field,
            std::string const& user_substring,
            std::vector<std::string> const& predefined_list,
            std::string& samplename
        ) const;
 
        bool match_header_sample_(
            std::string const& field,
            std::string const& user_substring,
            std::vector<std::string> const& predefined_list1,
            std::vector<std::string> const& predefined_list2,
            std::string& samplename
        ) const;
 
        bool match_header_sample_user_partial_(
            std::string const& field,
            std::string const& substring,
            std::string& samplename
        ) const;
 
        bool match_header_sample_predefined_partial_(
            std::string const& field,
            std::string const& substring,
            std::string& samplename
        ) const;
 
        // ---------------------------------------------
        //     Increment and Processing Samples
        // ---------------------------------------------
 
        void increment_();
 
        void process_sample_data_(
            SampleInfo const& sample_info,
            SampleData const& sample_data,
            Variant& variant,
            size_t sample_index
        );
 
        // -------------------------------------------------------------------------
        //     Data Members
        // -------------------------------------------------------------------------
 
    private:
 
        // Parent. If null, this indicates the end of the input and that we are done iterating.
        FrequencyTableInputStream const* parent_ = nullptr;
 
        // Data stream to read from.
        std::shared_ptr<utils::InputStream> input_stream_;
 
        // We keep information about the header, and which samples there are.
        // This also stores which sample has which index, so that we can access them in the Variant.
        HeaderInfo header_info_;
 
        // At the beginning, we guess which columns contain which information, based on their
        // headers. In order to keep the actual parsing as fast as possible, we use a lambda
        // for each column that is set up in the beginning. This avoids checking types
        // of columns for each row over and over again - we simply need to call all processors
        // in order, across the input file line.
        std::vector<std::function<void(genesis::utils::InputStream&)>> column_processors_;
 
        // We use a set of sample data objects to buffer values in, using a pointer to keep its
        // address stable (as iterators might be copied or moved implicitly while setting up
        // a for loop for instance). This is accessed via lambdas that capture the shared pointer.
        // This data is used to process them into the actual variant later.
        std::shared_ptr<std::vector<SampleData>> sample_data_;
 
        // Same reasoning as above: Using a pointer to keep the address of the Variant stable.
        // This is the final object that is exposed to the user.
        std::shared_ptr<Variant> current_variant_;
    };
 
    // ======================================================================================
    //      Main Class
    // ======================================================================================
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    FrequencyTableInputStream() = default;
 
    explicit FrequencyTableInputStream(
        std::shared_ptr<utils::BaseInputSource> input_source
    )
        // Call the other constuctor, to avoid code duplication.
        : FrequencyTableInputStream( input_source, std::unordered_set<std::string>{}, false )
    {}
 
    FrequencyTableInputStream(
        std::shared_ptr<utils::BaseInputSource> input_source,
        std::unordered_set<std::string> const& sample_names_filter,
        bool inverse_sample_names_filter = false
    )
        : input_source_( input_source )
        , sample_names_filter_( sample_names_filter )
        , inverse_sample_names_filter_( inverse_sample_names_filter )
    {}
 
    ~FrequencyTableInputStream() = default;
 
    FrequencyTableInputStream( self_type const& ) = default;
    FrequencyTableInputStream( self_type&& )      = default;
 
    self_type& operator= ( self_type const& ) = default;
    self_type& operator= ( self_type&& )      = default;
 
    // -------------------------------------------------------------------------
    //     Iteration
    // -------------------------------------------------------------------------
 
    Iterator begin() const
    {
        return Iterator( this );
    }
 
    Iterator end() const
    {
        return Iterator();
    }
 
    // -------------------------------------------------------------------------
    //     Basic Input Settings
    // -------------------------------------------------------------------------
 
    std::shared_ptr<utils::BaseInputSource> input_source() const
    {
        return input_source_;
    }
 
    self_type& input_source( std::shared_ptr<utils::BaseInputSource> value )
    {
        input_source_ = value;
        return *this;
    }
 
    std::unordered_set<std::string> const& sample_names_filter() const
    {
        return sample_names_filter_;
    }
 
    self_type& sample_names_filter( std::unordered_set<std::string> const& value )
    {
        sample_names_filter_ = value;
        return *this;
    }
 
    bool inverse_sample_names_filter() const
    {
        return inverse_sample_names_filter_;
    }
 
    self_type& inverse_sample_names_filter( bool value )
    {
        inverse_sample_names_filter_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Header Parsing Settings
    // -------------------------------------------------------------------------
 
    self_type& header_chromosome_string( std::string const& str )
    {
        usr_chr_name_ = str;
        return *this;
    }
 
    std::string const& header_chromosome_string() const
    {
        return usr_chr_name_;
    }
 
    self_type& header_position_string( std::string const& str )
    {
        usr_pos_name_ = str;
        return *this;
    }
 
    std::string const& header_position_string() const
    {
        return usr_pos_name_;
    }
 
    self_type& header_reference_base_string( std::string const& str )
    {
        usr_ref_name_ = str;
        return *this;
    }
 
    std::string const& header_reference_base_string() const
    {
        return usr_ref_name_;
    }
 
    self_type& header_alternative_base_string( std::string const& str )
    {
        usr_alt_name_ = str;
        return *this;
    }
 
    std::string const& header_alternative_base_string() const
    {
        return usr_alt_name_;
    }
 
    self_type& header_sample_reference_count_substring( std::string const& str )
    {
        usr_smp_ref_name_ = str;
        return *this;
    }
 
    std::string const& header_sample_reference_count_substring() const
    {
        return usr_smp_ref_name_;
    }
 
    self_type& header_sample_alternative_count_substring( std::string const& str )
    {
        usr_smp_alt_name_ = str;
        return *this;
    }
 
    std::string const& header_sample_alternative_count_substring() const
    {
        return usr_smp_alt_name_;
    }
 
    self_type& header_sample_frequency_substring( std::string const& str )
    {
        usr_smp_frq_name_ = str;
        return *this;
    }
 
    std::string const& header_sample_frequency_substring() const
    {
        return usr_smp_frq_name_;
    }
 
    self_type& header_sample_read_depth_substring( std::string const& str )
    {
        usr_smp_cov_name_ = str;
        return *this;
    }
 
    std::string const& header_sample_read_depth_substring() const
    {
        return usr_smp_cov_name_;
    }
 
    // -------------------------------------------------------------------------
    //     Settings
    // -------------------------------------------------------------------------
 
    std::shared_ptr<::genesis::sequence::ReferenceGenome> reference_genome() const
    {
        return ref_genome_;
    }
 
    self_type& reference_genome( std::shared_ptr<::genesis::sequence::ReferenceGenome> value )
    {
        ref_genome_ = value;
        return *this;
    }
 
    char separator_char() const
    {
        return separator_char_;
    }
 
    self_type& separator_char( char value )
    {
        separator_char_ = value;
        return *this;
    }
 
    std::string const& missing_value() const
    {
        return usr_missing_;
    }
 
    self_type& missing_value( std::string const& value )
    {
        usr_missing_ = value;
        return *this;
    }
 
    double int_factor() const
    {
        return int_factor_;
    }
 
    self_type& int_factor( double value )
    {
        if(
            value > max_int_factor_ ||
            static_cast<double>( static_cast<SampleCounts::size_type>( value )) != value
        ) {
            throw std::runtime_error(
                "Cannot set int_factor to " + std::to_string( value ) + " as this is out of range "
                "of the int type used for storing base counts."
            );
        }
        int_factor_ = value;
        return *this;
    }
 
    double allowed_relative_frequency_error() const
    {
        return allowed_rel_freq_error_;
    }
 
    self_type& allowed_relative_frequency_error( double value )
    {
        allowed_rel_freq_error_ = value;
        return *this;
    }
 
    bool frequency_is_ref() const
    {
        return frequency_is_ref_;
    }
 
    self_type& frequency_is_ref( bool value )
    {
        frequency_is_ref_ = value;
        return *this;
    }
 
    // -------------------------------------------------------------------------
    //     Data Members
    // -------------------------------------------------------------------------
 
private:
 
    // Input data.
    std::shared_ptr<utils::BaseInputSource> input_source_;
    std::unordered_set<std::string> sample_names_filter_;
    bool inverse_sample_names_filter_ = false;
 
    // Input settings.
    std::shared_ptr<::genesis::sequence::ReferenceGenome> ref_genome_;
    char separator_char_ = '\t';
 
    // Missing data indicators. We have a default set that we use, and a user-provided one.
    // Both are used case-insensitive. Need to be sorted so that `nan` comes before `na`,
    // as otherwise, `na` would match already when its actually `nan`, leading to an error.
    std::vector<std::string> missing_ = { ".", "nan", "na" };
    std::string usr_missing_;
 
    // When reading frequencies, for now, we want to turn them into counts, as this is what
    // our data infrastructure of Variant expects. To lose as little precision as possible, we
    // multiply the frequency [ 0.0, 1.0 ] by the largest integer for which itself and all
    // smaller integers can be stored in a double exactly. This guarantees that all frequencies
    // are mapped into the largest double range that is as exact as we can manage here.
    // Using these fake large counts also basically eliminates Bessel's correction for read depth
    // that we have in some of our pop gen equations, as we do not want that when working with
    // frequencies anyway, e.g., when obtaining them from HAF-pipe.
    // See https://stackoverflow.com/q/1848700/4184258 for the exact double value used here.
    static constexpr double max_int_factor_ = 9007199254740992.0;
 
    // We use a smaller factor by default, to make sure that we can add numbers without reaching
    // the max int precesion point of double, just in case. We are internally using size_t for
    // counts, so that should not happen, but prevents accidents. Also, using a power of ten here
    // it makes numbers converted from data nicer for users and potentially easier to understand,
    // as they more clearly correspond to the frequencies actually specified in the input data.
    static constexpr double default_int_factor_ = 1000000.0;
 
    // Make sure that this actually fits into the SampleCounts values.
    static_assert(
        static_cast<double>( static_cast<SampleCounts::size_type>( max_int_factor_ )) == max_int_factor_,
        "Numeric type for SampleCounts does not fit for FrequencyTableInputStream::max_int_factor_"
    );
 
    // The above is the max that we can use, but we allow users to set the used int factor that is
    // used for frequency-based computations, so that they can for example use the expected
    // read depth information from HAF-pipe, or similar values instead.
    // double int_factor_ = max_int_factor_;
    double int_factor_ = default_int_factor_;
 
    // When we have multiple pieces of information for a sample, we do cross checks, to make
    // sure that everything is in order. This relative value here is used for frequencies.
    double allowed_rel_freq_error_ = 0.001;
 
    // What does the frequency mean? We use: true = ref, false = alt frequency.
    // Probably would be more neat with an enum, but it's only used internally, so that's okay.
    bool frequency_is_ref_ = true;
 
    // Default names for header fields in a csv file that typically describe
    // the columns and quantities we are interested in. All used case-insensitively.
    std::vector<std::string> chr_names_ = { "chromosome", "chrom", "chr", "contig" };
    std::vector<std::string> pos_names_ = { "position", "pos" };
    std::vector<std::string> ref_names_ = { "reference", "referencebase", "ref", "refbase" };
    std::vector<std::string> alt_names_ = { "alternative", "alternativebase", "alt", "altbase" };
    std::vector<std::string> cnt_names_ = { "counts", "count", "cnt", "ct" };
    std::vector<std::string> frq_names_ = { "frequency", "freq", "maf", "af", "allelefrequency" };
    std::vector<std::string> cov_names_ = { "coverage", "cov", "readdepth", "depth", "ad" };
 
    // User supplied overwrites for the above automatic terms.
    // If either of them is given, we use those instead of the generic word lists.
    std::string usr_chr_name_;
    std::string usr_pos_name_;
    std::string usr_ref_name_;
    std::string usr_alt_name_;
    std::string usr_smp_ref_name_;
    std::string usr_smp_alt_name_;
    std::string usr_smp_frq_name_;
    std::string usr_smp_cov_name_;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard