Блог

Предположим, я ученый-геном, пытающийся хранить чрезвычайно длинные строки символов, каждый из которых представляет два бита информации (т. Е. Каждый элемент является либо G, A, T, либо C). Поскольку строки невероятно длинные, мне нужно иметь возможность хранить строку длиной N ровно в 2N битах (или, скорее, N / 4 байта).

Имея в виду эту мотивацию, я ищу обобщение std::bitset (или boost::dynamic_bitset<> ), которое работает с двухбитовыми значениями вместо однобитовых значений. Я хочу сохранить N такие двухбитовые значения, каждое из которых может быть 0, 1, 2 или 3. Мне нужны данные, упакованные как можно плотнее в памяти, поэтому vector<char> они не будут работать (поскольку они тратят впустую в 4 раза больше памяти).

Каков наилучший способ достижения моей цели? Один из вариантов — обернуть существующие шаблоны наборов битов настраиваемыми operator[] , итераторами и т. Д., Но я бы предпочел использовать существующую библиотеку, если это вообще возможно.

std::bitset<> это фиксированная длина, и вы, вероятно, этого не хотите.

Я думаю, вам следует продолжить и завершить std::vector<bool> .

Обратите внимание, что std::vector<bool> это оптимизировано для пространства, но имеет то преимущество, что оно динамическое по размеру. Предположительно, вам нужно откуда-то прочитать геном произвольной длины.

Подумайте, нужно ли вам много API для доступа к нему; вам может понадобиться всего пара методов.

Ответ @Jefffrey уже охватывает соответствующий код, если для bitset<> .

[Я не знаком boost::dynamic_bitset<> и что это может дать vector .]

Еще одна мысль заключается в том, может ли вам быть удобно работать с квадратиками букв, квадратиком, красиво заполняющим символ в пространстве.

 class Genome
{
public:
    enum class Letter {A,C,G,T};
    Genome(const std::stringamp; source)
    {
        code_.resize(source.size() * 2);
        for (unsigned index = 0; index != source.size();   index)
        {
            char text = source[index];
            Letter letter = textToLetter(text);
            set(index, letter);
        }
    }  
    static Letter textToLetter(char text)
    {
        // Or search through the array `letterText`.
        // Or come up with a neat but unintelligible one liner ...
        Letter letter = Letter::A;
        switch (text)
        {
        case 'A':
            letter = Letter::A;
            break;
        case 'C':
            letter = Letter::C;
            break;
        case 'G':
            letter = Letter::G;
            break;
        case 'T':
            letter = Letter::T;
            break;
        default:
            // Invalid - handle error.
            break;
        }
        return letter;
    }
    static char letterToText(Letter l) 
    {
        return letterText[(unsigned)l];
    }
    // Add bounds checking
    Letter get(unsigned index) const
    {
        unsigned distance = index * 2;
        char numeric = code_[distance]   code_[distance   1] * 2;
        return Letter(numeric);
    }
    // Add bounds checking
    void set(unsigned index, Letter value)
    {
        unsigned distance = index * 2;
        bool low = (unsigned)value amp; 1;
        bool high = (bool)((unsigned)value amp; 2);
        code_[distance] = low;
        code_[distance   1]  = high;
    }
    unsigned size()
    {
        return code_.size() / 2;
    }
    // Extend by numLetters, initially set to 'A'
    void extend(unsigned numLetters)
    {
        code_.resize(code_.size()   numLetters * 2);
    }
private:

    static char letterText[4];
    std::vector<bool> code_;
};

char Genome::letterText [4] = { 'A', 'C', 'G', 'T' };

int main()
{
    Genome g("GATT");
    g.extend(3);
    g.set(5, Genome::Letter::C);
    for (unsigned i = 0; i != g.size();   i)
        std::cout << Genome::letterToText(g.get(i));
    std::cout << std::endl;
    return 0;
}
  

У вас есть два варианта.

Учитывая:

 enum class nucleobase { a, c, g, t };
  

У вас есть два варианта. Вы можете:

• используйте один std::bitset и играйте с индексацией
• использовать std::bitset в сочетании с другим контейнером

Для первого вы можете просто определить пару функций, которые нацелены на правильное количество битов для каждого set / get:

 template<std::size_t N>
void set(std::bitset<N>amp; bits, std::size_t i, nucleobase x) {
    switch (x) {
        case nucleobase::a: bits.set(i * 2, 0); bits.set(i * 2   1, 0); break;
        case nucleobase::c: bits.set(i * 2, 0); bits.set(i * 2   1, 1); break;
        case nucleobase::g: bits.set(i * 2, 1); bits.set(i * 2   1, 0); break;
        case nucleobase::t: bits.set(i * 2, 1); bits.set(i * 2   1, 1); break;
    }
}

template<std::size_t N>
nucleobase get(const std::bitset<N>amp; bits, std::size_t i) {
    if (!bits[i * 2])
        if (!bits[i * 2   1]) return nucleobase::a;
        else                  return nucleobase::c;
    else
        if (!bits[i * 2   1]) return nucleobase::g;
        else                  return nucleobase::t;
}
  

Live demo

Приведенный выше пример — всего лишь ужасный пример (здесь почти 4 часа утра, и мне действительно нужно поспать).

Для второго вам просто нужно сопоставить аллели и биты:

 bit_pair bits_for(nucleobase x) {
    switch (x) {
        case nucleobase::a: return bit_pair("00"); break;
        case nucleobase::c: return bit_pair("10"); break;
        case nucleobase::g: return bit_pair("01"); break;
        case nucleobase::t: return bit_pair("11"); break;
    }
}

nucleobase nucleobase_for(bit_pair x) {
    switch (x.to_ulong()) {
        case 0: return nucleobase::a; break;
        case 1: return nucleobase::c; break;
        case 2: return nucleobase::g; break;
        case 3: return nucleobase::t; break;
        default: return nucleobase::a; break; // just for the warning
    }
}
  

Live demo

Конечно, если вам нужна длина времени выполнения, вы можете просто использовать boost::dynamic_bitset and std::vector .

Вот что я использую для k-мер фиксированной длины.

 #include <cstdint>
#include <cstdlib>
#include <ostream>

enum class nucleotide { A, C, G, T };

inline std::ostreamamp;
operator<<(std::ostreamamp; pOut, nucleotide pNt)
{
    switch (pNt) {
        case nucleotide::A: pOut << 'A'; break;
        case nucleotide::C: pOut << 'C'; break;
        case nucleotide::G: pOut << 'G'; break;
        case nucleotide::T: pOut << 'T'; break;
    }
    return pOut;
}

class kmer_base;

class nucleotide_proxy {
public:
    operator nucleotide() const {
        return nucleotide((*mWord >> (mPosition * 2)) amp; 3);
    };

    nucleotide_proxyamp; operator=(nucleotide pNt) {
        uint64_t word = *mWord;
        word amp;= ~(uint64_t(3) << (mPosition*2));
        word |= uint64_t(pNt) << (mPosition*2);
        *mWord = word;

        return *this;
    };

private:
    friend class kmer_base;

    nucleotide_proxy(uint64_t* pWord, uint8_t pPosition)
        : mWord(pWord), mPosition(pPosition)
    {
    }

    uint64_t* mWord;
    uint8_t mPosition;
};


class kmer_base {
protected:
    nucleotide_proxy access(uint64_t* pWord, size_t pPosition)
    {
        return nucleotide_proxy(pWord   (pPosition / 32), (pPosition amp; 31));
    }

    const nucleotide_proxy access(uint64_t* pWord, size_t pPosition) const
    {
        return nucleotide_proxy(pWord   (pPosition / 32), (pPosition amp; 31));
    }
};


template<int K>
class kmer : public kmer_base
{
    enum { Words = (K   31) / 32 };
public:
    nucleotide_proxy operator[](size_t pOutdex) {
        return access(mWords, pOutdex);
    }

    const nucleotide_proxy operator[](size_t pOutdex) const {
        return access(mWords, pOutdex);
    }

private:
    uint64_t mWords[Words];
};
  

Расширение этого до k-mere динамической длины оставлено в качестве упражнения; это довольно легко, если у вас есть nucleotide_proxy в вашем распоряжении. Эффективная реализация оператора обратного дополнения также остается в качестве упражнения.