Anyway, what's this post about? Well, while studying for my biology exam I came across the Polymerase Chain Reaction (PCR) and couldn't help my programmer senses from tingling. I'll first explain in a few words what the PCR method is and then how it can be implemented in code.
The Polymerase Chain Reaction, developed by Kary Mullis in the 60's, is a method of replicating a specific DNA sequence many, many (many, many, many) times. The point is, when you want to analyze a certain region of an individual's DNA, say to examine the gene responsible for eye color, you have to isolate that gene from the thousands of other genes found on your typical chromosome. So how do you go about doing that? Polymerase Chain Reaction. Let me throw in a few other necessary definitions and chunks of information about DNA:
A Nucleotide is the basic building block of the double helix that makes up our DNA. It consists of one of four organic bases (as well as a sugar and a phosphate molecule). Those four organic bases are Adenine, Cytosine, Guanine and Thymine (A, C, G, T for short). Because the DNA structure consists of two strands -- left and right -- there are always two complementary nucleotides next two each other. A always goes together with T and C always holds hands with G (aw). Therefore, when we speak of a strand of DNA, e.g. A-C-C-T, there is always a complementary strand, in this case T-G-G-A (always the complementary nucleotides). Moreover, it's important to note that DNA strands have a direction. One end of a DNA strand is called the 5' (five prime) end and the other the 3' (three prime) end (this naming stems from the chemical structure of the molecules). When new nucleotides are added to a DNA strand, this must always happen in the 5' to 3' direction. So if we wanted to add a new Guanine nucleotide to the DNA sequence from above, it would have to be like this: (5') A-C-C-T (3') + G = (5') A-C-C-T-G (3'). Note also that the left and right strand of a DNA structure run antiparallel (opposite directions), meaning the left strand has it's 5' prime end at the opposite end of the right strand's 5' end. Lastly, I need to mention primers and Polymerase. Polymerase is an enzyme that complements a strand, while primers are enzymes that tell the Polymerase at what position of a strand to start complementing.
Let's pull it all together. Have a look at this DNA strand:
(5') T-G-A-A-C-T-G-A-C-A-T-T (3')
Say we're really interested in the region between the 3rd (A) and the 7th (G) nucleotide (including those), because it tells us some information about a hereditary disease. The problem is that these nucleotides are about the size of Kim Kardashian's brain. So just a few molecules wide, around 2 nanometers actually. We can't do any meaningful research with just these four nucleotides. We need more of them. Millions and millions and millions more. This is what the Polymerase Chain Reaction does, by creating complementary strands with the Polymerase enzyme I just introduced. So let's have a further look at this problem. To copy the region between the 3rd and 7th locus (= position in a DNA segment), we first set a primer (P) to the 8th locus:
(5') T-G-A-A-C-T-G-A-C-A-T-T (3')
(P)
Then we throw some Polymerase into the equation. It'll look for the primer and synthesize a complementary strand from the primer onwards (that's why it's at locus 8, because 8 is then no longer included in the complementary strand). Note that this synthesis will go into the opposite direction of the old strand, because the two DNA strands must run antiparallel. Moreover, you'll notice that the synthesis goes all the way to the end of the old strand and doesn't just stop at locus 3, even though we'd want it to. Unfortunately, there is no "stop-primer":
(5') T-G-A-A-C-T-G-A-C-A-T-T (3')
(3') A-C-T-T-G-A-C (5')
Okey dokey, almost there. That was iteration one and we almost have our desired region between locus 3 and locus 7 isolated. However, we're still stuck with the old, super-long strand from the beginning and our new, still-too-long complementary strand. Let's ignore the old one for now and just focus on the one that's a bit shorter (the complementary strand). What we can do next is simple, we just set another primer to the new strand, this time on locus 2:
(3') A-C-T-T-G-A-C (5')
(P)
And let the Polymerase do its job:
(3') A-C-T-T-G-A-C (5')
(5') A-A-C-T-G (3')
Yes! We have our desired strand and can now eradicate cancer. Ok maybe not just yet. But anyway, notice that from now on, every single complementary strand we make of this new strand will be exactly perfect. You could almost call it ... a chain reaction.
Additionally, don't forget that we can always re-use the longer strands to make new "perfect" / "correct" ones, which will in turn again produce new correct strands and so on. Exponential!
So how can we implement this in code? Quite simple. We need to model the DNA as a class, consisting of two strands, which in turn consist of nucleotides, which in turn consist of an organic base as well as some information about its locus. This DNA model should give us access to the two strands -- left and right -- as well as make it possible to create a complementary strand for either of the two. Moreover, the DNA class' interface should make it easy to access the 3' and 5' end of either strands and provide a means of iterating over their nucleotides in the correct direction (5' to 3'). Lastly, there should also be an easy way of initializing the DNA class' left and right strands, either from another strand or directly from a sequence of bases. This description paves the way for the following C++ class:
Declaration:
#include <vector>
#include <deque>
class DNA
{
public:
enum class Base { A, C, G, T };
typedef std::size_t locus_t;
struct Nucleotide { Base base; locus_t locus; };
typedef std::vector<Base> sequence_t;
typedef std::deque<Nucleotide> strand_t;
typedef std::vector<strand_t> container_t;
typedef strand_t::iterator left_iterator;
typedef strand_t::const_iterator const_left_iterator;
typedef strand_t::reverse_iterator right_iterator;
typedef strand_t::const_reverse_iterator const_right_iterator;
DNA();
DNA(const sequence_t& left);
DNA(const sequence_t& left,
const sequence_t& right);
DNA(const strand_t& left);
DNA(const strand_t& left,
const strand_t& right);
void setLeftStrand(const strand_t& strand);
void setLeftStrand(const sequence_t& bases);
strand_t getLeftStrand();
const strand_t getLeftStrand() const;
left_iterator getLeftFivePrime();
left_iterator getLeftThreePrime();
const_left_iterator getLeftFivePrime() const;
const_left_iterator getLeftThreePrime() const;
strand_t getLeftComplementary(locus_t primer = 0) const;
void setRightStrand(const strand_t& strand);
void setRightStrand(const sequence_t& bases);
strand_t getRightStrand();
const strand_t getRightStrand() const;
right_iterator getRightFivePrime();
right_iterator getRightThreePrime();
const_right_iterator getRightFivePrime() const;
const_right_iterator getRightThreePrime() const;
strand_t getRightComplementary(locus_t primer = 0) const;
locus_t rightToLeftLocus(locus_t rightLocus) const;
locus_t leftToRightLocus(locus_t leftLocus) const;
private:
strand_t left_;
strand_t right_;
};
And definition:
#include "pcr.hpp"
#include <stdexcept>
DNA::DNA() = default;
DNA::DNA(const sequence_t& left)
{
setLeftStrand(left);
setRightStrand(getLeftComplementary(left.size() - 1));
}
DNA::DNA(const sequence_t& left,
const sequence_t& right)
{
setLeftStrand(left);
setRightStrand(right);
}
DNA::DNA(const strand_t& left)
: left_(left),
right_(getLeftComplementary(left_.size() - 1))
{ }
DNA::DNA(const strand_t& left,
const strand_t& right)
: left_(left),
right_(right)
{ }
void DNA::setLeftStrand(const strand_t& strand)
{
left_ = strand;
}
void DNA::setLeftStrand(const std::vector
The interface you see here matches the description I gave above. I'll nevertheless talk about two non-obvious things about this class: first, why the Nucleotide locus is hardcoded and, secondly, how the DNA::getLeftComplementary() and DNA::getRightComplementary() methods work.
Our DNA strands are std::deque<Nucleotide> (typedef strand_t, deque and not vector because we can push_front), where each Nucleotide consists of a member of the DNA::Nucleotide enum class, as well as a locus. The locus is hardcoded for each Nucleotide and not simply contained in the Nucleotide's index in the std::deque, because when DNA strands are complemented, they reduce in size and thus a primer locus wouldn't simultaneously be valid for the original and the complementary strand. As an example, if you set the primer locus to position 10 of a 15-nucleotide-long original strand, this would make the complementary strand 5 nucleotides long (indices 0-4) and thus locus 10 is invalid for the complementary strand. By hardcoding the locus for each nucleotide and not relying on the strand index, locus 10 would always work, because the complementary strand's nucleotides copy the original strand's hardcoded loci (original strand's hardcoded loci go from 0 to 14, complementary strand's loci then go from 10 to 14 → locus 10 valid for both).
Now to DNA::getLeftComplementary(). Its purpose is to take a locus and then return the complementary strand of the DNA object's left strand, starting from that locus (think of this locus as the primer position, just that it is actually included in the complementary strand, unlike the primer's position in real world PCR). At first, we search for the locus in the left strand using std::find_if(), to which we pass a lambda as a predicate which will make std::find_if look at the strand's nucleotide loci and return true if any locus matches the primer locus. We throw an exception if the locus isn't found. Else, we enter a for loop where we get the complementary base for each of the left strand's nucleotide bases (starting at the primer locus) and add this complementary base alongside the current locus (taken from left strand's nucleotide) as a nucleotide to the complementary strand (which is, don't forget, just a Nucleotide container). We iterate backwards from the primer locus because the complementary strand runs antiparallel, so it must move towards the 5' end of the original (left) strand.
Notice also that we call complementary.push_front and not push_back, to simulate this antiparallel behavior. You see, our right and left strand's are actually stored antiparallel in the DNA class, so we must append Nucleotides as well as iterate over them in different directions. Therefore, the reason why DNA::left_iterator is a typedef of strand_t::iterator (strand_t is the std::deque<Nucleotide>) and DNA::right_iterator is a typedef of strand_t::reverse_iterator, is that DNA::getLeftFivePrime() returns left_.begin(), so a normal forward iterator, while DNA::getRightFivePrime() returns right_.rbegin(), which is a reverse_iterator that actually points to the last valid Nucleotide in the strand (would be right_.end() - 1 as a normal iterator). We can then iterate over both strands from their respective 5' ends to their respective 3' ends in the same way.
So much to the DNA class. Let's talk about the algorithm for the actual PCR function. It'll be recursive. It'll take a DNA object and return a container of strands containing (mostly) the target sequence, which we create using the Polymerase Chain Reaction algorithm. It'll also take a "first" and a "last" locus, that tell us where the target DNA sequence starts and ends ("first" and "last" are included in the sequence). The PCR function's last parameter will be the number of iterations that should be performed for the algorithm.
Have a look a this tree diagram. It shows the development of the strands of each DNA strand pair for each iteration of the PCR function. A check mark stands for a strand that is "correct", which we defined before as any strand that is the target sequence. Conversely, a cross denotes an "incorrect" strand, which is any strand that is not the target sequence (such as the initial DNA pair).
The initial DNA strand pair is iteration 0 of the algorithm, the second row is iteration 1, the third iteration 2, the fourth is iteration 3 and the last is iteration 4. Regarding some notation to address single strands in this tree, we'll have one number stand for the iteration (0 is the top of the tree, 4 is the bottom here), a second number to to stand for the strand pair (left to right) and "l" or "r" is then either the left or the right strand of that pair, respectively. To give some examples:
00l) is in the 0th row of this tree (= 0th iteration), the 0th pair, the left strand (so the left strand of the initial DNA pair at the top of this tree)
37r) is in iteration 3, the 7th pair, the right strand (right most strand of this iteration)
22l) is in iteration 2, the 2nd pair, the left strand
... and so on.
We can make some generalizations for this algorithm. What we can say for sure, with little thought to it, is that after n iterations this algorithm produces a total of 2^(n+1) strands (both correct and incorrect). This is clear since the number of strands doubles for each iteration. The + 1 stems from the fact that we start with 2^1 strands and not 2^0 initially.
Next up, let's find out how many correct strands there are after each iteration. One solution is to find out how many incorrect strands there are after n iterations and to subtract that number from the total.
So how many incorrect strands do we get for n iterations? 2 * (n + 1) is the answer. Why? Have a look at the tree. Notice that during each iteration, two incorrect strands are added from the far left and far right (e.g. 10r creates 21r; 23l creates 37l etc.). The reason why is that the left and right most strands are always still as long as the original pair, so when you perform the PCR on them, you get a complementary strand that goes either from the first target locus to the end of the strand, or from the last target locus to the beginning of the strand. So they're still incorrect (though they yield a correct strand in the next iteration). Just to visualize, left and right most strands look something like this below, where | denotes either the last or first locus of the target sequence we want to copy:
00l, 10l, 20l, ...) 5' <---|-----|---> 3' (left most strand)
00r, 11r, 23r, ...) 3' <---|-----|---> 5' (right most strand)
When we perform the first PCR cycle on them, you'll get complementary strands that are almost correct. So for iteration 1, the two new incorrect strands are (don't forget complementary strands run antiparallel):
10r) 3' <---|-----| 5' (complementary to left most strand)
11l) 5' |-----|---> 3' (complementary to right most strand)
These are the two incorrect strands that are added during each iteration. In 2 * (n + 1), the 2 therefore stems from the fact that we progress in steps of 2 (since 2 incorrect strands are added for each n). The n + 1 is there because we already have 2 incorrect strands for iteration 0 (at the start).
So we know that the total number of strands, correct and incorrect, for each iteration is 2^(n+1) and that the number of incorrect strands for each iteration is 2(n + 1). Therefore, the first solution to finding the number of correct strand for each iteration is: 2^(n+1) - 2(n+1)
The second approach to finding the number of correct strands is to look directly at the correct strands produced per iteration and how they develop with time. Going back to the visualization from above, we notice that the two incorrect complementary strands will yield correct strands in the next iteration:
21r) 3' <---|-----| 5' (previously 10r)
21l) 5' |-----| 3' (complementary to 21r / 10r)
23l) 5' |-----|---> 3' (previously 11l)
23r) 3' |-----| 5' (complementary to the 23l / 11l)
We can say that starting with iteration 1 (n=1), every incorrect strand yields a correct strand, except for the the incorrect strands at the far left and far right, since those yield the two new incorrect strands we just spoke about (e.g. 10r or 11l), which keep on wandering down the tree producing those correct strands (10r -> 21r -> 32r ...). But just to clarify, those newly created incorrect strands like 10r 11l on the far left and far right mean that the number of correct strands produced by incorrect strands also increments by 2 per iteration (new incorrect strand in n means new correct strand in n+1). Therefore, we can say that for n iterations where n is greater zero, 2(n-1) correct strands are added solely by complementing incorrect ones.
So those are the correct strands produced from incorrect ones. What about correct strands produced from correct ones? Those will be correct again, of course:
21l) 5' |-----| 3'
=
32l) 3' |-----| 5' (complementary of 21l)
32r) 3' |-----| 5' (previously 21l)
This time, the number of correct strands doubles. For the next iteration, they will double again, thereafter they double again and so on. This is clearly a case for a power of two. So we know two that the number of newly added correct strands per iteration is 2(n-1). We also know that these 2(n-1) strands will subsequently double for each iteration. Therefore, two iterations on, there will be 2 * 2 * 2(n-1) = 2^3(n-1). We can generalize and say that after m iterations, where m equals 1 initially, there will be 2^m * (n-1) strands created from those initial 2(n-1). Notice that this also works at the very beginning when m is 1: 2^1 * (n-1) = 2(n-1).
What we can conclude from this is that to calculate the number of correct strands directly (i.e. without just subtracting the number of incorrect from the total) for a single iteration n, we have to count from n back to the first iteration and perform this operation of doubling the single, newly created correct strands as we just did above, but for each of the previous iterations. Let me put it like this: in iteration n we have 2(n-1) newly added correct strands. Those haven't doubled yet. In the previous iteration, so n-1, we had 2(n-2) newly added correct strands. Those have doubled, so in iteration n those 2(n-2) correct strands have turned into 2*2(n-2) = 2^2(n-2) strands. In iteration n-2 we had 2(n-3) correct strands, those doubled in n-1 and then doubled again for n. So in n, those 2(n-3) have turned into 2*2*2(n-3) = 2^3*(n-3) strands. If we do this for every iteration starting at n, going back to 1, we get the total number of correct strands present in iteration n:
Phew. Hope you survived that. To answer the question you're asking yourself: no, that knowledge won't really be of any use for the C++ PCR algorithm I'm about to show you. But thanks for the fun discussion!
Here's the actual PCR non-member non-friend function:
DNA::container_t pcr(const DNA& dna,
DNA::locus_t first,
DNA::locus_t last,
std::size_t n)
{
// If we reached the maximum number of iterations,
// return the DNA strand pair
if (! n) return {dna.getLeftStrand(), dna.getRightStrand()};
else
{
// Create left branch of this position in the tree
DNA left {dna.getLeftStrand(), dna.getLeftComplementary(last)};
// rvalue reference to capture the return value
DNA::container_t&& leftReturned = pcr(left, first, last, n - 1);
// Same thing for the right hand side
DNA right {dna.getRightComplementary(first), dna.getRightStrand()};
DNA::container_t&& rightReturned = pcr(right, first, last, n - 1);
// Tie together right and left returned values
leftReturned.insert(leftReturned.end(), rightReturned.begin(), rightReturned.end());
// Then return the two
return leftReturned;
}
}
Ok, maybe I'm actually wrong. Those last few paragraphs are of use. They should make this function really, really easy to understand. What we do is: move left down the tree whenever we can, else back and then right. Essentially, we create complementary DNA strand pairs moving down the tree as long as we can, then return all the DNA pairs from the very bottom (when the parameter n is zero).
Test run:
#include "pcr.hpp"
#include <iostream>
int main(int argc, char * argv [])
{
// Our whole DNA. Note that this is
// the left strand, the constructor
// initializes the right as the left's
// complementary
DNA dna({DNA::Base::A,
DNA::Base::C,
DNA::Base::G,
DNA::Base::T,
DNA::Base::C});
// Do the PCR using our dna, with the first
// locus of our target sequence at 1 (= C)
// and the last locus at 3 (= T). Perform
// 4 iterations.
DNA::container_t ret = pcr(dna, 1, 3, 4);
std::size_t correct = 0, incorrect = 0;
for (auto i : ret)
{
// Suffices to do the check
if (i.size() == 3) ++correct;
else ++incorrect;
for (auto j : i)
{
switch (j.base)
{
case DNA::Base::A: std::cout << "A\t"; break;
case DNA::Base::C: std::cout << "C\t"; break;
case DNA::Base::G: std::cout << "G\t"; break;
case DNA::Base::T: std::cout << "T\t"; break;
}
}
std::cout << std::endl;
}
std::cout << std::endl << "Correct: " << correct << std::endl;
std::cout << "Incorrect: " << incorrect << std::endl;
}
Output:
A C G T C T G C A C G T T G C A C G T G C A C G T T G C A C G T G C A C G T G C A C G T G C A C G T T G C A C G T C G C A C G T G C A C G T G C A C G T G C A C G T C G C A C G T G C A C G T C G C A C G T C T G C A G Correct: 22 Incorrect: 10
Note that it doesn't matter if a strand is the complementary of the target sequence, it still counts as correct because biologically speaking complementary strand and original strand are the same thing. The bases you see here are the bases of the DNA strands at the bottom of the tree from before, just turned by 270°. As discussed above, the far left and far right (here very top and very bottom) strands are the original strands (ACGTC), the other strands of those two pairs as well as a few others in the tree are the "almost correct ones" (TGCA / CGTC) and all those strands with three nucleotides are our target sequence: CGT / GCA (same thing).
What about the numbers? As we said (n is 4, don't forget):
Number of correct and incorrect strands after n iterations:
2^(n+1) = 2^(4+1) = 32 ... 22 + 10 = 32 ✓
Number of incorrect strands after n iterations:
2(n+1) = 2(4+1) = 2 * 5 = 10 ✓
Number of correct strands after n iterations:
2^(n+1) - 2(n+1) = 32 - 10 = 22 ✓
or
2^1*(n-1) + 2^2*(n-2) + 2^3*(n-3) = 2 * 3 + 4 * 2 + 8 * 1 = 22 ✓
We can thus conclude:
Total number of blog posts about modeling the Polymerase Chain Reaction in C++: 1
Thanks for reading, hope you learnt something :) kthxbye
No comments :
Post a Comment