genome_region_list.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_GENOME_REGION_LIST_H_
#define GENESIS_POPULATION_GENOME_REGION_LIST_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 <cassert>
#include <cstdint>
#include <map>
#include <stdexcept>
#include <string>
#include <vector>
 
#include "genesis/population/genome_locus.hpp"
#include "genesis/population/genome_region.hpp"
#include "genesis/utils/containers/interval_tree.hpp"
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Genome Data
// =================================================================================================
 
struct EmptyGenomeData
{};
 
// =================================================================================================
//     Genome Region List
// =================================================================================================
 
// template<class DataType = EmptyGenomeData>
class GenomeRegionList
{
public:
 
    // -------------------------------------------------------------------------
    //     Typedefs and Enums
    // -------------------------------------------------------------------------
 
    // using data_type = DataType;
    using data_type = EmptyGenomeData;
    using numerical_type = size_t;
    using tree_type = genesis::utils::IntervalTree<
        data_type, numerical_type, genesis::utils::IntervalClosed
    >;
    using self_type = GenomeRegionList;
 
    using iterator               = typename tree_type::iterator;
    using const_iterator         = typename tree_type::const_iterator;
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    GenomeRegionList() = default;
    ~GenomeRegionList() = default;
 
    GenomeRegionList( GenomeRegionList const& ) = default;
    GenomeRegionList( GenomeRegionList&& )      = default;
 
    GenomeRegionList& operator= ( GenomeRegionList const& ) = default;
    GenomeRegionList& operator= ( GenomeRegionList&& )      = default;
 
    // -------------------------------------------------------------------------
    //     Modifiers
    // -------------------------------------------------------------------------
 
    // -------------------------------------------
    //         Add Manually
    // -------------------------------------------
 
    void add( std::string const& chromosome )
    {
        // We use the special value 0 to denote that we want the whole chromosome.
        add( chromosome, 0, 0 );
    }
 
    void add(
        std::string const& chromosome,
        numerical_type start,
        numerical_type end,
        bool overlap = false
    ) {
        // Check chromosome.
        if( chromosome.empty() ) {
            throw std::invalid_argument(
                "Cannot add region to GenomeRegionList with empty chromosome name, "
                "as this denotes an invalid chromosome."
            );
        }
 
        // Check positions.
        // The start and end are also checked in the interval tree, but let's do it here
        // so that the error message is nicer in case they are wrong.
        if( start > end ) {
            throw std::invalid_argument(
                "Cannot add region to GenomeRegionList with start == " +
                std::to_string( start ) + " > end == " + std::to_string( end )
            );
        }
        if(( start == 0 ) ^ ( end == 0 )) {
            throw std::invalid_argument(
                "Cannot add region to GenomeRegionList with either start == 0 or end == 0, "
                "but not both, as we use 1-base indexing, with both being 0 being interpreted "
                "as the special case of denoting the whole chromosome. "
                "Hence either both start and end have to be 0, or neither."
            );
        }
 
        // Insert, either by merging with an existing, or as a new region.
        // We just get the chromosome from the map via array access, which will create it
        // if it is not yet present.
        if( overlap ) {
            regions_[ chromosome ].insert_overlap({ start, end });
        } else {
            regions_[ chromosome ].insert({ start, end });
        }
    }
 
    // -------------------------------------------
    //         Add Locus
    // -------------------------------------------
 
    void add( GenomeLocus const& locus, bool overlap = false )
    {
        add( locus.chromosome, locus.position, locus.position, overlap );
    }
 
    void add( GenomeLocus const& start, GenomeLocus const& end, bool overlap = false )
    {
        if( start.chromosome != end.chromosome ) {
            throw std::invalid_argument(
                "Cannot use two GenomeLocus instances with different chromosomes ( start == \"" +
                start.chromosome + "\", end == \"" + end.chromosome + "\") as an entry in a "
                "GenomeRegionList."
            );
        }
        add( start.chromosome, start.position, end.position, overlap );
    }
 
    // -------------------------------------------
    //         Add Region
    // -------------------------------------------
 
    void add( GenomeRegion const& region, bool overlap = false )
    {
        add( region.chromosome, region.start, region.end, overlap );
    }
 
    // -------------------------------------------
    //         Add Region List
    // -------------------------------------------
 
    void add( GenomeRegionList const& other, bool overlap = false )
    {
        for( auto const& chr : other.regions_ ) {
            for( auto const& interval : chr.second ) {
                add( chr.first, interval.low(), interval.high(), overlap );
            }
        }
    }
 
    // -------------------------------------------
    //         Clear
    // -------------------------------------------
 
    void clear()
    {
        regions_.clear();
    }
 
    void clear( std::string const& chromosome )
    {
        if( regions_.count( chromosome ) == 0 ) {
            throw std::invalid_argument(
                "Chromosome name \"" + chromosome + "\" not found in GenomeRegionList"
            );
        }
        regions_.erase( chromosome );
    }
 
    // -------------------------------------------------------------------------
    //     Locus Queries
    // -------------------------------------------------------------------------
 
    bool is_covered( std::string const& chromosome, numerical_type position ) const
    {
        // Using find(), so we only have to search in the map once, for speed.
        auto const it = regions_.find( chromosome );
        if( it == regions_.end() ) {
            return false;
        }
        auto const& chrom_tree = it->second;
 
        // If the chromosome in our interval tree contains the 0 interval, we consider that
        // as having the whole chromosome covered.
        if( chrom_tree.overlap_find( 0 ) != chrom_tree.end() ) {
            return true;
        }
 
        // If the above is not the case, check the actual position.
        return chrom_tree.overlap_find( position ) != chrom_tree.end();
    }
 
    bool is_covered( std::string const& chromosome ) const
    {
        auto const it = regions_.find( chromosome );
        if( it == regions_.end() ) {
            return false;
        }
        auto const& chrom_tree = it->second;
        return chrom_tree.overlap_find( 0 ) != chrom_tree.end();
    }
 
    size_t cover_count(
        std::string const& chromosome,
        numerical_type position,
        bool whole_chromosome = false
    ) const {
        // Using find(), so we only have to search in the map once, for speed.
        auto const it = regions_.find( chromosome );
        if( it == regions_.end() ) {
            return 0;
        }
        auto const& chrom_tree = it->second;
 
        // Count overlaps.
        size_t result = 0;
        chrom_tree.overlap_find_all( position, [&]( const_iterator it ){
            (void) it;
            assert( it->interval().within( position ));
            ++result;
            return true;
        });
 
        // Add whole chromosome if needed.
        if( whole_chromosome && ( chrom_tree.overlap_find( 0 ) != chrom_tree.end() )) {
            ++result;
        }
 
        return result;
    }
 
    // Not used at the moment, as we have no access to the end iterator to check for a valid find.
 
    // /* *
    //  * @brief Return an iterator to the GenomeRegion that covers the given position on the
    //  * chromosome, or an end iterator if there is no region that covers the position.
    //  *
    //  * Throws an exception if the chromosome is not in the list at all.
    //  */
    // const_iterator find( std::string const& chromosome, size_t position ) const
    // {
    //     // Using find(), so we only have to search in the map once, for speed.
    //     auto const it = regions_.find( chromosome );
    //     if( it == regions_.end() ) {
    //         throw std::invalid_argument(
    //             "GenomeRegionList does not contain chromosome \"" + chromosome + "\""
    //         );
    //     }
    //     auto const& chrom_tree = it->second;
    //     return chrom_tree.overlap_find( position );
    // }
 
    // -------------------------------------------------------------------------
    //     Chromosome Accessors
    // -------------------------------------------------------------------------
 
    bool empty() const
    {
        return regions_.empty();
    }
 
    size_t chromosome_count() const
    {
        return regions_.size();
    }
 
    std::vector<std::string> chromosome_names() const
    {
        std::vector<std::string> result;
        for( auto const& p : regions_ ) {
            result.push_back( p.first );
        }
        return result;
    }
 
    bool has_chromosome( std::string const& chromosome ) const
    {
        return regions_.count( chromosome ) > 0;
    }
 
    tree_type const& chromosome_regions( std::string const& chromosome ) const
    {
        return regions_.at( chromosome );
    }
 
    tree_type& chromosome_regions( std::string const& chromosome )
    {
        return regions_.at( chromosome );
    }
 
    size_t region_count( std::string const& chromosome ) const
    {
        if( regions_.count( chromosome ) == 0 ) {
            throw std::invalid_argument(
                "Chromosome name \"" + chromosome + "\" not found in GenomeRegionList"
            );
        }
        return regions_.at( chromosome ).size();
    }
 
    size_t total_region_count() const
    {
        size_t cnt = 0;
        for( auto const& reg : regions_ ) {
            cnt += reg.second.size();
        }
        return cnt;
    }
 
    std::map<std::string, tree_type> const& chromosome_map() const
    {
        return regions_;
    }
 
    std::map<std::string, tree_type>& chromosome_map()
    {
        return regions_;
    }
 
    // -------------------------------------------------------------------------
    //     Iterator
    // -------------------------------------------------------------------------
 
    // iterator begin()
    // {
    //     return regions_.begin();
    // }
 
    // const_iterator begin() const
    // {
    //     return regions_.begin();
    // }
 
    // iterator end()
    // {
    //     return regions_.end();
    // }
 
    // const_iterator end() const
    // {
    //     return regions_.end();
    // }
 
    // -------------------------------------------------------------------------
    //     Settings
    // -------------------------------------------------------------------------
 
    // bool allow_overlap() const
    // {
    //     return allow_overlap_;
    // }
    //
    // self_type& allow_overlap( bool value )
    // {
    //     when set this way, it will only affect newly added regions...
    //     better in constructor and const?
    //     allow_overlap_ = value;
    //     return *this;
    // }
 
    // -------------------------------------------------------------------------
    //     Data Members
    // -------------------------------------------------------------------------
 
private:
 
    std::map<std::string, tree_type> regions_;
 
    // bool allow_overlap_ = true;
 
};
 
} // namespace population
} // namespace genesis
 
#endif // include guard