kmer.hpp

Go to the documentation of this file.

#ifndef GENESIS_SEQUENCE_KMER_KMER_H_
#define GENESIS_SEQUENCE_KMER_KMER_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 <array>
#include <cassert>
#include <climits>
#include <cstdint>
#include <stdexcept>
#include <string>
 
namespace genesis {
namespace sequence {
 
// =================================================================================================
//     Kmer
// =================================================================================================
 
struct KmerTagDefault{};
 
// TODO template params: Tag, Encoding, Storage
// encoding is a policy class that does all the: how many bits, which bits stand for what char, which alphabet, how to convert the alpahbet --- all of that!
// storage defaults to a simple wrapper around 64bit uint that gives access to consecutive bits according to the encoding bit width.
// but also have a std array class that allows to do the same on larger k mers
 
template<typename Tag = KmerTagDefault>
class Kmer
{
public:
 
    // -------------------------------------------------------------------------
    //     Typedefs and Enums
    // -------------------------------------------------------------------------
 
    using WordType = uint64_t;
 
    static const size_t BIT_SIZE = sizeof(WordType) * CHAR_BIT;
    static_assert( CHAR_BIT == 8, "CHAR_BIT != 8" );
 
    static const uint8_t BITS_PER_CHAR = 2;
 
    static const uint8_t MAX_CHARS_PER_KMER = BIT_SIZE / BITS_PER_CHAR;
    static_assert( BIT_SIZE % BITS_PER_CHAR == 0, "BIT_SIZE % BITS_PER_CHAR != 0" );
 
    // -------------------------------------------------------------------------
    //     Constructors and Rule of Five
    // -------------------------------------------------------------------------
 
    // Constructor to initialize k-mer data
    Kmer( WordType data )
        : data_( data )
    {
        // assert( data & dna_ones_mask_ == data );
    }
 
    // -------------------------------------------------------------------------
    //     K
    // -------------------------------------------------------------------------
 
    static uint8_t k()
    {
        return k_;
    }
 
    static void set_k( uint8_t k )
    {
        if( k_ != 0 && k != k_ ) {
            throw std::runtime_error( "Cannot set k for a given Tag multiple times");
        }
        reset_k( k );
    }
 
    static void reset_k( uint8_t k )
    {
        if( k == 0 || k > MAX_CHARS_PER_KMER ) {
            throw std::invalid_argument( "Cannot use k-mer with k==" + std::to_string( k ));
        }
        k_ = k;
    }
 
    // -------------------------------------------------------------------------
    //     Data
    // -------------------------------------------------------------------------
 
    uint8_t operator[] ( size_t position ) const
    {
        assert( k_ > 0 && k_ <= MAX_CHARS_PER_KMER );
        if( position >= k_ ) {
            throw std::runtime_error(
                "Invalid position " + std::to_string( position ) +
                " in k-mer of size " + std::to_string( k_ )
            );
        }
        return ( data_ & dna_char_mask_[position] ) >> ( position * BITS_PER_CHAR );
    }
 
    WordType const& value() const
    {
        return data_;
    }
 
    WordType& value()
    {
        return data_;
    }
 
    // -------------------------------------------------------------------------
    //     Internal Data
    // -------------------------------------------------------------------------
 
private:
 
    static const WordType all_0_ = 0;
    static const WordType all_1_ = (((1ul << 32) - 1) << 32)  +  ((1ul << 32) - 1);
 
    static const std::array<WordType, MAX_CHARS_PER_KMER> dna_char_mask_;
    // static const std::array<WordType, MAX_CHARS_PER_KMER> dna_char_mask_ =
    // {{
    //     1ul << 0  | 1ul << 1,  1ul << 2  | 1ul << 3,  1ul << 4  | 1ul << 5,  1ul << 6  | 1ul << 7,
    //     1ul << 8  | 1ul << 9,  1ul << 10 | 1ul << 11, 1ul << 12 | 1ul << 13, 1ul << 14 | 1ul << 15,
    //     1ul << 16 | 1ul << 17, 1ul << 18 | 1ul << 19, 1ul << 20 | 1ul << 21, 1ul << 22 | 1ul << 23,
    //     1ul << 24 | 1ul << 25, 1ul << 26 | 1ul << 27, 1ul << 28 | 1ul << 29, 1ul << 30 | 1ul << 31,
    //     1ul << 32 | 1ul << 33, 1ul << 34 | 1ul << 35, 1ul << 36 | 1ul << 37, 1ul << 38 | 1ul << 39,
    //     1ul << 40 | 1ul << 41, 1ul << 42 | 1ul << 43, 1ul << 44 | 1ul << 45, 1ul << 46 | 1ul << 47,
    //     1ul << 48 | 1ul << 49, 1ul << 50 | 1ul << 51, 1ul << 52 | 1ul << 53, 1ul << 54 | 1ul << 55,
    //     1ul << 56 | 1ul << 57, 1ul << 58 | 1ul << 59, 1ul << 60 | 1ul << 61, 1ul << 62 | 1ul << 63
    // }};
    // static_assert( dna_char_mask_[MAX_CHARS_PER_KMER-1] != 0, "dna_char_mask_ == 0" );
 
    // static const std::array<WordType, MAX_CHARS_PER_KMER> dna_ones_mask_ =
    // {{
    //     all_0_,       all_1_ >> 63, all_1_ >> 62, all_1_ >> 61,
    //     all_1_ >> 60, all_1_ >> 59, all_1_ >> 58, all_1_ >> 57,
    //     all_1_ >> 56, all_1_ >> 55, all_1_ >> 54, all_1_ >> 53,
    //     all_1_ >> 52, all_1_ >> 51, all_1_ >> 50, all_1_ >> 49,
    //     all_1_ >> 48, all_1_ >> 47, all_1_ >> 46, all_1_ >> 45,
    //     all_1_ >> 44, all_1_ >> 43, all_1_ >> 42, all_1_ >> 41,
    //     all_1_ >> 40, all_1_ >> 39, all_1_ >> 38, all_1_ >> 37,
    //     all_1_ >> 36, all_1_ >> 35, all_1_ >> 34, all_1_ >> 33,
    //     all_1_ >> 32, all_1_ >> 31, all_1_ >> 30, all_1_ >> 29,
    //     all_1_ >> 28, all_1_ >> 27, all_1_ >> 26, all_1_ >> 25,
    //     all_1_ >> 24, all_1_ >> 23, all_1_ >> 22, all_1_ >> 21,
    //     all_1_ >> 20, all_1_ >> 19, all_1_ >> 18, all_1_ >> 17,
    //     all_1_ >> 16, all_1_ >> 15, all_1_ >> 14, all_1_ >> 13,
    //     all_1_ >> 12, all_1_ >> 11, all_1_ >> 10, all_1_ >> 9,
    //     all_1_ >> 8,  all_1_ >> 7,  all_1_ >> 6,  all_1_ >> 5,
    //     all_1_ >> 4,  all_1_ >> 3,  all_1_ >> 2,  all_1_ >> 1
    // }};
 
    // -------------------------------------------------------------------------
    //     Private Members
    // -------------------------------------------------------------------------
 
private:
 
    static uint8_t k_;
    WordType data_;
 
};
 
template<typename Tag>
uint8_t Kmer<Tag>::k_ = 0;
 
template<typename Tag>
const std::array<typename Kmer<Tag>::WordType, Kmer<Tag>::MAX_CHARS_PER_KMER> Kmer<Tag>::dna_char_mask_ =
{{
    1ul << 0  | 1ul << 1,  1ul << 2  | 1ul << 3,  1ul << 4  | 1ul << 5,  1ul << 6  | 1ul << 7,
    1ul << 8  | 1ul << 9,  1ul << 10 | 1ul << 11, 1ul << 12 | 1ul << 13, 1ul << 14 | 1ul << 15,
    1ul << 16 | 1ul << 17, 1ul << 18 | 1ul << 19, 1ul << 20 | 1ul << 21, 1ul << 22 | 1ul << 23,
    1ul << 24 | 1ul << 25, 1ul << 26 | 1ul << 27, 1ul << 28 | 1ul << 29, 1ul << 30 | 1ul << 31,
    1ul << 32 | 1ul << 33, 1ul << 34 | 1ul << 35, 1ul << 36 | 1ul << 37, 1ul << 38 | 1ul << 39,
    1ul << 40 | 1ul << 41, 1ul << 42 | 1ul << 43, 1ul << 44 | 1ul << 45, 1ul << 46 | 1ul << 47,
    1ul << 48 | 1ul << 49, 1ul << 50 | 1ul << 51, 1ul << 52 | 1ul << 53, 1ul << 54 | 1ul << 55,
    1ul << 56 | 1ul << 57, 1ul << 58 | 1ul << 59, 1ul << 60 | 1ul << 61, 1ul << 62 | 1ul << 63
}};
 
} // namespace sequence
} // namespace genesis
 
#endif // include guard