function.hpp

Go to the documentation of this file.

#ifndef GENESIS_SEQUENCE_KMER_FUNCTION_H_
#define GENESIS_SEQUENCE_KMER_FUNCTION_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/sequence/kmer/kmer.hpp"
 
#include <cassert>
#include <cstdint>
#include <iostream>
#include <stdexcept>
#include <string>
#include <type_traits>
 
namespace genesis {
namespace sequence {
 
// =================================================================================================
//     Helper Functions
// =================================================================================================
 
inline size_t number_of_kmers( size_t k, size_t alphabet_size = 4 )
{
    // We do an explicit loop instead of using std::pow,
    // in case that double precision is not enough here.
    size_t n = 1;
    for( size_t i = 0; i < k; ++i ) {
        n *= alphabet_size;
    }
    return n;
}
 
// =================================================================================================
//     Kmer Operators
// =================================================================================================
 
template<typename Tag>
bool operator==( Kmer<Tag> const& lhs, Kmer<Tag> const& rhs )
{
    return lhs.value() == rhs.value();
}
 
template<typename Tag>
bool operator!=( Kmer<Tag> const& lhs, Kmer<Tag> const& rhs )
{
    return lhs.value() != rhs.value();
}
 
template<typename Tag>
bool operator<=( Kmer<Tag> const& lhs, Kmer<Tag> const& rhs )
{
    return lhs.value() <= rhs.value();
}
 
template<typename Tag>
bool operator>=( Kmer<Tag> const& lhs, Kmer<Tag> const& rhs )
{
    return lhs.value() >= rhs.value();
}
 
template<typename Tag>
bool operator<( Kmer<Tag> const& lhs, Kmer<Tag> const& rhs )
{
    return lhs.value() < rhs.value();
}
 
template<typename Tag>
bool operator>( Kmer<Tag> const& lhs, Kmer<Tag> const& rhs )
{
    return lhs.value() > rhs.value();
}
 
template<typename Tag>
std::string to_string( Kmer<Tag> const& kmer )
{
    auto result = std::string( kmer.k(), 'X' );
    for( size_t i = 0; i < kmer.k(); ++i ) {
        switch( kmer[i] ) {
            case 0x00: {
                result[i] = 'A';
                break;
            }
            case 0x01: {
                result[i] = 'C';
                break;
            }
            case 0x02: {
                result[i] = 'G';
                break;
            }
            case 0x03: {
                result[i] = 'T';
                break;
            }
        }
    }
    return result;
}
 
template<typename Tag>
std::ostream& operator<< ( std::ostream& os, Kmer<Tag> const& kmer )
{
    os << to_string( kmer );
    return os;
}
 
// =================================================================================================
//     Nucleotide Kmer Functions
// =================================================================================================
 
template<typename Tag>
Kmer<Tag> reverse_complement( Kmer<Tag> const& kmer )
{
    // Function is written for a specific bit width.
    static_assert(
        std::is_same<typename Kmer<Tag>::WordType, std::uint64_t>::value,
        "Kmer::WordType != uint64_t"
    );
 
    // Adapted from Kraken2 at https://github.com/DerrickWood/kraken2/blob/master/src/mmscanner.cc
    // which itself adapted this for 64-bit DNA use from public domain code at
    // https://graphics.stanford.edu/~seander/bithacks.html#ReverseParallel
 
    // Reverse bits (leaving bit pairs intact, as those represent nucleotides):
    // Swap consecutive pairs, then nibbles, then bytes, then byte pairs, then halves of 64-bit word
    auto value = kmer.value();
    value = (( value & 0xCCCCCCCCCCCCCCCCUL ) >>  2) | (( value & 0x3333333333333333UL ) <<  2 );
    value = (( value & 0xF0F0F0F0F0F0F0F0UL ) >>  4) | (( value & 0x0F0F0F0F0F0F0F0FUL ) <<  4 );
    value = (( value & 0xFF00FF00FF00FF00UL ) >>  8) | (( value & 0x00FF00FF00FF00FFUL ) <<  8 );
    value = (( value & 0xFFFF0000FFFF0000UL ) >> 16) | (( value & 0x0000FFFF0000FFFFUL ) << 16 );
    value = ( value >> 32 ) | ( value << 32);
 
    // Finally, complement, and shift to correct position, removing the invalid lower bits.
    return Kmer<Tag>( (~value) >> ( Kmer<Tag>::BIT_SIZE - kmer.k() * Kmer<Tag>::BITS_PER_CHAR ));
}
 
template<typename Tag>
Kmer<Tag> canonical_representation( Kmer<Tag> const& kmer )
{
    auto rev_comp = reverse_complement( kmer );
    return kmer < rev_comp ? kmer : rev_comp;
}
 
} // namespace sequence
} // namespace genesis
 
#endif // include guard