10
\$\begingroup\$

Introduction

I have this small C++ program that packs a DNA string over alphabet ACGT into a bit vector, two bits per character.

Code

dnapack.hpp

#ifndef NET_CODERODDE_DNAPACK_HPP
#define NET_CODERODDE_DNAPACK_HPP

#include <string>
#include <vector>

namespace net {
namespace coderodde {
namespace dna {

    class packed_dna_sequence {
    public:
        packed_dna_sequence(std::string& dna_sequence);
        std::string&& unpack() const;
        std::vector<bool> get_packed_sequence() const;

    private:
        std::vector<bool> m_packed_bits;
    };

} // End of namespace net::coderodde::dna.
} // End of namesapce net::coderodde.
} // End of namespace net.

#endif // NET_CODERODDE_DNAPACK_HPP

dnapack.cpp

#include "dnapack.hpp"
#include <algorithm>
#include <cctype>
#include <iostream>
#include <stdexcept>
#include <string>
#include <utility>

static std::string&& dna_to_upper_case(std::string& unpacked_dna_sequence)
{
    std::string ret;
    std::for_each(unpacked_dna_sequence.cbegin(),
                  unpacked_dna_sequence.cend(),
                  [&ret](char c) { ret.push_back(std::toupper(c)); });

    return std::move(ret);
}

namespace net {
namespace coderodde {
namespace dna {

    using std::string;

    packed_dna_sequence::packed_dna_sequence(string& unpacked_dna_sequence)
    {
        std::string upper_case_unpacked_dna_sequence =
        dna_to_upper_case(unpacked_dna_sequence);
        m_packed_bits.reserve(2 * upper_case_unpacked_dna_sequence.length());

        for (char c : upper_case_unpacked_dna_sequence) {
            switch (c) {
                case 'A':
                    m_packed_bits.push_back(false);
                    m_packed_bits.push_back(false);
                    break;

                case 'C':
                    m_packed_bits.push_back(false);
                    m_packed_bits.push_back(true);
                    break;

                case 'G':
                    m_packed_bits.push_back(true);
                    m_packed_bits.push_back(false);
                    break;

                case 'T':
                    m_packed_bits.push_back(true);
                    m_packed_bits.push_back(true);
                    break;

                default:
                    throw std::invalid_argument("Unknown character.");
            }
        }
    }

    std::vector<bool> packed_dna_sequence::get_packed_sequence() const
    {
        // Copy. We don't want the client to tamper with the bits.
        return m_packed_bits;
    }

    std::string&& packed_dna_sequence::unpack() const
    {
        std::string ret;
        ret.reserve(m_packed_bits.size() / 2);

        for (size_t index = 0; index < m_packed_bits.size(); index += 2)
        {
            bool bit1 = m_packed_bits[index + 1];
            bool bit2 = m_packed_bits[index];
            size_t num = (bit2 << 1) | bit1;

            switch (num)
            {
                case 0:
                    ret += 'A';
                    break;

                case 1:
                    ret += 'C';
                    break;

                case 2:
                    ret += 'G';
                    break;

                case 3:
                    ret += 'T';
                    break;
            }
        }

        return std::move(ret);
    }

} // End of namespace net::coderodde::dna.
} // End of namespace net::coderodde.
} // End of namespace net.

main.cpp

#include "dnapack.hpp"
#include <algorithm>
#include <iostream>
#include <iterator>
#include <string>

bool same_dna(std::string& dna1, std::string& dna2)
{
    if (dna1.length() != dna2.length())
    {
        return false;
    }

    for (size_t i = 0; i < dna1.length(); ++i)
    {
        if (std::toupper(dna1[i]) != std::toupper(dna2[i]))
        {
            return false;
        }
    }

    return true;
}

int main(int argc, const char * argv[]) {
    using net::coderodde::dna::packed_dna_sequence;

    std::string dna;
    std::cout << "> ";
    std::cin >> dna;

    packed_dna_sequence pds(dna);

    std::vector<bool> rep = pds.get_packed_sequence();

    std::cout << "Binary representation: ";
    std::ostream_iterator<int> out_it(std::cout);
    std::copy(rep.cbegin(), rep.cend(), out_it);
    std::cout << '\n';

    std::string unpacked_dna = pds.unpack();
    std::cout << "Unpacked DNA: "
              << unpacked_dna
               << '\n';

    std::cout << "Agrees: "
              << std::boolalpha
              << same_dna(dna, unpacked_dna)
              << '\n';
}

Critique request

I would like to hear comments regarding:

  • const correctness,
  • move semantics,
  • resource management,
  • coding style,
  • API design.
\$\endgroup\$
  • \$\begingroup\$ Just for the sake of saying it, have you heard of the Burrows-Wheeler transform and its use in the enhanced suffix array? \$\endgroup\$ – miscco Aug 13 '17 at 21:30
  • \$\begingroup\$ Using pairs of vector<bool> elements may compile to much less efficient code than if you implemented 2-bit elements yourself, on top of uintmax_t or something. (Or on top of uint8_t, but uintmax_t or uint32_t is more likely to compile to more efficient asm when copying multiple elements). It's likely that some compilers (even with max optimization enabled) will still compile m_packed_bits[index + 1] and m_packed_bits[index]; to mask the bits separately and then merge them, rather than doing the asm equivalent of (bits[index/8] >> bit_offset) & 0x3. \$\endgroup\$ – Peter Cordes Aug 14 '17 at 4:05
  • \$\begingroup\$ You're worrying about move semantics when you're already shooting yourself in the foot by processing string data. \$\endgroup\$ – user997112 Jun 21 at 0:29
10
\$\begingroup\$
  1. static std::string&& dna_to_upper_case(std::string& unpacked_dna_sequence). That's wrong. It returns a dangling reference. You can't return an rvalue reference to a temporary object. Return it by value. It doesn't change the unpacked_dna_sequence, so it would be better to pass it by a const reference.

  2. The same is true for std::string&& packed_dna_sequence::unpack() const. It's also broken because it returns a reference to a temporary object.

  3. bool same_dna(std::string& dna1, std::string& dna2). Once again, if you don't change the input objects, pass them by a const reference. You can also simplify this function using std::equal function with a custom predicate.

  4. This piece of code is a little bit confusing:

    std::string ret;
    std::for_each(unpacked_dna_sequence.cbegin(),
                  unpacked_dna_sequence.cend(),
                  [&ret](char c) { ret.push_back(std::toupper(c)); });
    

    std::for_each is usually used for something with side effects. If you want to transform one sequence into another, std::transform is an idiomatically better choice.

\$\endgroup\$
8
\$\begingroup\$

First things first:

Undefined behavior

Kraskevich pointed at the most obvious one. Yet, there is one more (although it will work on X86, from what I've been told). Basically, do not pass chars directly to any function in cctype, convert them to unsigned char first. More info.

vector bool

Very interesting decision. I'd rather pack it myself using fallthrough switch statement into std::uint8_t. Duff's device never gets old.

Supported operations

I believe scientists would want to do more operations on it. I'm not sure which ones, but definitely building some algorithms on top of this would be useful. This leads to the next point:

What does it represent?

If one wants to pack it, I'd guess they want to store it in some persistent storage, but there is no serialization/deserialization functions. Passing around into functions? May be, but it is not clear yet if it is beneficial. So, what I want to say is, the class doesn't have clear purpose.

Smaller things:

  • Constructor can take the sequence by const reference, since it is not modifying it.

  • Return by value, let people/compiler decide if they want to copy. I believe in the worst case they'll need to add std::move(), though at the return site it will prevent NRVO.

  • No equivalence operator. I believe it will be useful.

Conclusion

The code looks like it was written in a rush. It needs rethinking.

\$\endgroup\$
  • \$\begingroup\$ What comes to persistence, I was sort of anticipating it, and, therefore, added the get_packed_sequence so the client programmer could, say, send it over network, etc. \$\endgroup\$ – coderodde Aug 13 '17 at 13:41
  • \$\begingroup\$ @coderodde, I think adding a more straightforward way would be great. People can still mess with the vector before serializing it. Some people could accidentally do something erroneous, so preventing it would be great. Also, it will allow you to write more efficient serialization/deserialization. \$\endgroup\$ – Incomputable Aug 13 '17 at 13:43
  • \$\begingroup\$ Equivalence operators = operator==, operator!=, ...? \$\endgroup\$ – coderodde Aug 13 '17 at 14:42
  • \$\begingroup\$ Also, can you tell me some other methods I should implement on behalf of your 3rd point? \$\endgroup\$ – coderodde Aug 13 '17 at 14:43
  • \$\begingroup\$ @coderodde, yes, those are equivalence operators. I'd want vector behavior, but I'm not bioinformatics student. \$\endgroup\$ – Incomputable Aug 13 '17 at 16:18
7
\$\begingroup\$

The Interface

class packed_dna_sequence {
public:
    packed_dna_sequence(std::string& dna_sequence);
    std::string&& unpack() const;
    std::vector<bool> get_packed_sequence() const;

private:
    std::vector<bool> m_packed_bits;
};

Does not look very useful. You have a full dna_sequence as a string already (so you have the space). When you unpack the sequence you have to build the full string again.

So does it really save space? Not really; you must have the full size string to create it and you can only get values from it after de-compressing the whole thing back to a string.

So the only thing this structure is good for is pre-materalizing a packed structure before serialization to some storage.

This is not a great design -- memory is much more scarce than offline storage. So your design seems to fail on all these fronts. If you wanted to optimize for storage then you should write a compressing stream.

If you want to optimize for memory then you should have an object that allows random (or serialized) access to the members while it is in the compressed form.

Code Review

Const Reference Parameters

Sure pass by reference (it avoids the copy).

        packed_dna_sequence(std::string& dna_sequence);

But if you are not going to modify the input then pass by const reference.

R-Value return

Looks like a good move.

        std::string&& unpack() const;

But in practice does nothing. Return by value and RVO (or NRVO) will automatically kick in and build the object in place at the destination (so technically a tick faster as you don't even need to swap pointers).

Return by Value.

        std::vector<bool> get_packed_sequence() const;

Should that not return a const reference to the internal object.

Leaky Abstraction

        std::vector<bool> get_packed_sequence() const;

Also, why do you want to expose the raw packed sequence? This only exposes your implementation in binds you to maintaining it. This is known as a leaky abstraction. Never leak your internal implementation details.

Copying the string?

Your uppercasing function makes a copy of the complete string. That seems rather wasteful. Why not do it as you use the letters. That would seem a lot less costly and would not require you to loop over the string multiple times.

static std::string&& dna_to_upper_case(std::string& unpacked_dna_sequence)
{
    std::string ret;
    std::for_each(unpacked_dna_sequence.cbegin(),
                  unpacked_dna_sequence.cend(),
                  [&ret](char c) { ret.push_back(std::toupper(c)); });

    return std::move(ret); // Not needed NRVO will kick in here.
}
\$\endgroup\$
  • \$\begingroup\$ Compact data representations are useful for cache footprint / memory bandwidth. Available RAM limitations are not the only reason for compact data structures. Or maybe you want to free that string before loading another from disk, so you only have one unpacked string in memory at any one time, but hundreds of packed DNA sequences (which you can then compare against each other for equality with 1/4 the memory / cache bandwidth of one base per byte). \$\endgroup\$ – Peter Cordes Aug 14 '17 at 4:11
  • \$\begingroup\$ The OP hasn't yet implemented any operators like that yet, but even being able to decompress any chosen sequences on the fly into buffers that stay hot in cache could be better than the memory bandwidth of reading a couple std::string or char[] sequences. Anyway, the idea has potential, but this implementation and API is not good. Great points, especially about random access or iterators over the compressed data. \$\endgroup\$ – Peter Cordes Aug 14 '17 at 4:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.