This is the second in a series of three blog posts describing our solution to a bioinformatics problem from Rosalind.info, Problem BA1(i) (Find most frequent words with mismatches in a string). To solve this problem and generate variations of a DNA string as required, we implemented a recursive backtracking method in the Go programming language.

• Part 1: Counting Variations
• Part 2: Generating Variations (you are here)
• Part 3: Go Implementation of Recursive Backtracking

Table of Contents

• Problem Description
• Permutations vs Combinations vs Variations
• Recursion
  • Recursive Backtracking Pseudocode
• Appying to DNA Variations
  • Generating Visits with Binary Numbers
  • Assembling the Variation

Problem Description

The task at hand is to take a given input strand of DNA, and generate variations from it that have up to \(d\) differences (a Hamming distance of \(d\)) in the codons (base pairs).

In part 1 of this series, we walk through the construction of an analytical formula to count the number of variations of a given DNA string that can be generated, given the constraints of the problem.

In part 2 of this series, we cover several techniques to generate variations on a DNA string, and present pseudocode for the recursive backtracking method that we use here.

In part 3 of this series, we will cover our implementation of the recursive backtracking method in the Go programming language.

Permutations vs Combinations vs Variations

Before covering generation of variations of a DNA string, we should cover some terminology for clarification.

If we were to use the term permutations, as in, we are counting (or generating) permutations of the input DNA string, that would imply that we were doing some kind of rearrangement of the elements of the input DNA string (for example, swapping two codons). This is not the problem that we are solving, and requires different formulas. (See Permutations entry on Wolfram MathWorld.)

The variations that we are referring to are not exactly combinations, either, though. If we were to use the term combinations, it would imply that we were choosing a set of \(k\) integers from a set of \(d\) integers \({1, 2, \dots, d}\).

The variations that we are counting are similar to combinations, but with the additional act of swapping out each codon at the position (integer) selected with three other possible codons, so there are more variations than combinations (and many more permutations than variations).

Transforming the Problem Space

A surprisingly large variety of problems in combinatorics can be transformed into an equivalent problem involving binary numbers, which are usually easier to think about.

To generate variations, we can break up the process of producing a variation into two steps, or choices, and then convert these choices (and the process of making them) into an equivalent problem in terms of binary numbers.

We can decompose the cration of a DNA string variation into the first step of choosing which codons (indices) to edit, and the second step of cycling through every possible codon (ATGC) at the selected indices.

To translate this into an equivalent binary number problem, consider the input string of DNA "AAAAA" and let the Hamming distance that we are considering be \(d = 1\). Then we can code each index with a 0 (not chosen) or a 1 (chosen) and turn the problem into cycling throgh all binary numbers with 1 bit:

00001
00010
00100
01000
10000

The second step is to cycle through each alternate codon at the given position, so that 00001 would generate the variations:

AAAAC
AAAAG
AAAAT

and so on.

We saw this two-part technique already when counting the total number of variations that could be created in Part 1: Counting Variations. It resulted in a counting formula with two terms, a binomial term for step 1 and an exponential term for step 2.

We can think of the problem as forming a tree with several decision nodes that need to be explored; this type of problem structure is ideal for a recursive backtracking algorithm.

We will cover the use of recursive backtracking to actually explore the entire tree of possible outcomes (not just count it), starting with some review and background on recursive backtracking and how it works.

Recursion

Recursion is a common pattern to use for problems that require exploring a large problem space that requires us to make several selections.

A recursive backtracking algorithm is analogous to exploring a maze but laying out a rope as you go, so tht you can revisit each possible route. In this case, we are using backtracking to make the choice of which indices of the input DNA string to modify. We want to explore all possible choices to generate all possible variations of the input DNA string, and backtracking gives us the framework to do that.

For example, if we wanted to recursively generate codon choices for the case of an input DNA string like "AAAAA" and \(d = 2\), we would call a recursive method twice; the first time through, we would choose one of the five indices, and mark it as picked; then we would call the method again, and choose a second index (different from the first) and mark it as picked.

When unrolled, this is equivalent to a nested for loop,

for i in range( 0 .. len(dna_string) ):
    for j in range( 0 .. len(dna_string) ):
        if (i != j):
            Start with the binary number 00000
            Set the digit at index i to 1
            Set the digit at index j to 1

Recursive Backtracking Pseudocode

Basic pseudocode for a backtracking method:

explore method:
    base case:
        visit this solution
    recursive case:
        for each available choice:
            make a choice
            explore outcomes
            unmake the choice
            move on to the next choice

Applying to DNA Variations

There are actually two places where we need to apply backtracking to our problem.

Generating Visits with Binary Numbers

The first application of recursive backtracking is to carry out step 1, choosing indices in the original DNA string to modify or cycle through altnerate codons. We showed above how generating variations on a kmer of length \(k\) at a distance \(d\) from the original kmer was equivalent to generating binary numbers with \(d\) bits set to 1.

We can use recursive backtracking to generate these numbers. By creating a method that recursively selects an index to switch to 1, and passing that (and all prior choices) on to further recursive calls, the function can recurse to a given depth \(d\) and visit all possible binary numbers with \(d\) bits set to 1.

The base case of this recursive method would be reached when all \(d\) choices had been made and \(d\) bits were set to 1. Then the choice of indices to swap out with alternate codons would be passed on to a recursive method that would carry out Step 2 (see below).

For example, to generate variations of the 5-mer AAAAA, we would start by selecting a Hamming distance \(d\), then generate a binary number with \(d\) bits set to 1 to select indices to modify. Suppose \(d = 2\); then the first few binary numbers are:

AAAAA
11000
10100
10010
10001
...

To expand on the pseudocode a bit more, to generate a binary number with \(d\) bits flipped to 1 we will want to call a recursive method with a depth of \(d\), making a choice at each recursive call of which index to set to 1 next.

The \(n^{th}\) recursive call picks the \(n^{th}\) index for 1. Each index can only be chosen once in the stack of recursive calls, and the indices that have been chosen by prior recursive function calls are passed along.

Thus we need a minimum of two parameters: an integer indicating the depth level of this recursive function call, and an integer array of index choices.

function generate_binary_numbers( depth, choices[], ... ):

    if depth is 0,
        base case
        no more choices left to make
        choices[] is full
        pass along choices[] to assemble the variations

    else,
        recursive case
        for each possible index,
            if this index is not already in choices,
                add this index to choices
                generate_binary_numbers( depth+1, choices[] )
                remove this index from choices

Assembling the Variation

Each binary number is then turned into variations by substituting every combination of 3 codons in every position with a 1 possible, so the first binary number for \(d=2\) would generate the variations:

AAAAA
11000
-----
CCAAA
GCAAA
TCAAA
CGAAA
GGAAA
TGAAA
CTAAA
GTAAA
TTAAA

This would be repeated for all Hamming distances up to the maximum specified Hamming distance.

Like the generation of binary numbers, the substitution of all possible combinations of codons at these positions is a task conducive to a recursive backtracking algorithm.

Like the prior task's recursive method, this task's recursive method will have one parameter for depth (number of choices left to make) and a range of choices to try (codons).

function assemble_variations( depth, choices[], ... ):

    if depth is 0,
        base case
        no more choices left to make
        choices[] is full
        pass along choices[] to assemble the variations

    else,
        recursive case
        for each possible index,
            if this index is not already in choices,
                add this index to choices
                generate_binary_numbers( depth+1, choices[] )
                remove this index from choices

In the final part, Part 3, of this blog post, we will cover the actual Go implementation of these functions.

This is the first in a series of three blog posts describing our solution to a bioinformatics problem from Rosalind.info, Problem BA1(i) (Find most frequent words with mismatches in a string). To solve this problem and generate variations of a DNA string as required, we implemented a recursive backtracking method in the Go programming language.

• Part 1: Counting Variations (you are here)
• Part 2: Generating Variations
• Part 3: Go Implementation of Recursive Backtracking

Table of Contents

• Problem Description
• Useful Functions
  • Binomial Function
  • Factorial Function
• Counting Permutations
  • Deriving the Formula
  • Term 1: Picking DNA Indices
  • Term 1: Side Note on Ordering
  • Term 2: Modifying DNA Codons
• Final Counting Formula
• Implementing in Go
  • Binomial and Factorial Functions in Go
  • Variations Counting Function in Go

Problem Description

The task at hand is to take a given input strand of DNA, and generate variations from it that have up to \(d\) differences (a Hamming distance of \(d\)) in the codons (base pairs).

In part 1 of this series, we walk through the construction of an analytical formula to count the number of variations of a given DNA string that can be generated, given the constraints of the problem.

In part 2 of this series, we cover several techniques to generate variations on a DNA string, and present pseudocode for the recursive backtracking method that we use here.

In part 3 of this series, we will cover our implementation of the recursive backtracking method in the Go programming language.

Useful Functions

It's always useful to review some basic mathematics useful for combinatorics applications. We'll review the factorial and binomial functions, which will show up in our final formula for the total nubmer of variations we will be generating.

Factorial Function

The factorial function for an integer \(n\) is written \(n!\) and is defined for \(n \geq 1\) as:

for example, \(5!\) would be:

Binomial Function

The binomial function has many applications in combinatorics. It is the number of ways of independently selecting \(k\) items from a set of \(n\) items, and is written:

Counting Permutations

What we want is a formula to count the number of permutations

To derive a formula, it helps to think through the problem starting with smaller special cases, and generalize from there in terms of the problem parameters.

Deriving the Formula

The problem we're trying to solve is generating all perms with hamming distance less than or equal to d, but let's start with a simpler problem: generating all perms with hamming distance of exactly d.

Then we can just sum up over each d.

Start with a simple situation: string of dna with 3 codons. Case of hamming distance 0 too trivial, so start with case of hamming distance of 1.

There are two terms in our combinatorics formula that we need to think about:

• Term 1: We have a certain number of codons to modify (this is fixed by the Hamming distance d that we pick). Term 1 counts up the number of ways of selecting which indices of the original DNA string to modify.

• Term 2: Once we've picked out the indices we are going to modify, we have several variations for each index (4 total codons, so 3 variations). Term 2 is a count of the number of variations that are possible, given the choice of indices in the original DNA string to modify.

The approach here is to think about these two terms independently and separately. Each term has a formula to count the number of possibilities indexed by each. Then, because these are independent choices, the total number of combined choices is the product of these two terms.

Term 1: Picking DNA Indices

The first term in our formula for number of variations will be the term representing the number of ways of choosing which indices in the original DNA input string to edit.

Given a Hamming distance of \(d\), and the fact that we have one and only one edit (Hamming distance unit) per base pair, Term 1 counts the number of ways of picking \(d\) items from a set of \(n\) items. Order does not matter.

This problem is equivalent to having a row of \(n\) on/off switches, all in the off position, and counting the number of ways of throwing exactly \(d\) of them into the on position.

Likewise, it is equivalent to having \(d\) identical colored balls, and counting the number of ways of placing them into \(n\) slots, one ball per slot.

We can see how the problem has a kind of triangular structure. Returning to the scenario of \(n\) on/off switches:

• If we have \(d = n\) switches to throw, or if we have \(d = 0\) switches to throw, in either case we have only 1 possible outcome.

• If we have \(d = n-1\) switches to throw, or if we have \(d = 1\) switch to throw, either way we have \(n\) possible outcomes

• If we have \(d = n-2\) or \(d = 2\) switches to throw, there are \(n (n-1)\) possible outcomes; etc.

In fact, this problem - choosing \(d\) things from a set of \(n\) things - is common enough that there is a special function just to describe it, and that's the binomial function (covered above).

The binomial function is defined as:

In the scenarios posed above, the order of our choices did not matter - the balls were not numbered, the order in which we threw each switch did not affect the outcome.

If the order did matter, if the order in which the on/off switches were thrown mattered or if the balls that were placed into slots had sequential numbers on them, then we would need a different function - the expression above to count the number of outcomes would not have a \(k!\) in the denominator.

Term 1: Side Note on Ordering

If the order of the index choices does not matter, the \(k!\) term in the denominator must be included to cancel out double-counting in the situations where (for example) \(i\) is chosen first and \(j\) is chosen second, and then the situation where \(j\) is chosen first and \(i\) is chosen second.

If the \(k!\) term is present in the denominator, it says that the order in which items are selected does not matter, in which case we are generating combinations.

To count combinations, use the "n choose k" function. See the Combination article on Wolfram MathWorld.

If the \(k!\) term is not present in the denominator, it says that the order in which items are selected does matter, in which case we are generating permutations.

To count permutations, use the "n pick k" function. See the Permutation article on Wolfram MathWorld.

Term 2: Modifying DNA Codons

Once we've selected the \(d\) indices in the original DNA string that we are going to modify, we have to count the number of ways those base pairs can be modified.

We have \(d\) base pairs to modify, and \(c = 4\) total codons (ATGC). Each base pair that we are modifying has \(c-1\) possible codons that it we can swap it out with, and each choice is independent, so the number of possibile outcomes (Term 2) is:

Final Counting Formula

To write the final formula for counting the number of variations \(V\) of a given DNA string of length \(n\) that are a Hamming distance of less than or equal to \(d\), with \(c\) possible codons (A, T, G, C), we will need to sum over Hamming distances from 0 to \(d\):

Implementing in Go

Now, let's look at how we would implement this counting formula in Go. This will be useful, since programs run much faster when they are able to allocate all the sapce they need in memory ahead of time. Counting the number of variations of our DNA input string will allow us to do just that.

Binomial and Factorial Functions in Go

We'll start with binomial and factorial functions in Go: continuing with our theme of recursion, we implement a recursive factorial function.

// Compute the factorial of an integer.
func Factorial(n int) int {
    if n < 2 {
        // base case
        return 1
    } else {
        // recursive case
        return n * Factorial(n-1)
    }
}

The factorial function will behave correctly for the case of \(n=1\) and \(n=0\), and will return 1 if \(n\) is negative (which is reasonable behavior for our purposes.)

The binomial function utilizes the factorial function:

// Returns value of the binomial coefficient Binom(n, k).
func Binomial(n, k int) int {

    result := 1

    // Since C(n, k) = C(n, n-k)
    if k > (n - k) {
        k = n - k
    }

    // Calculate value of:
    //
    // ( n * (n-1) * ... * (n-k+1) )
    // -----------------------------
    //   ( k * (k-1) * ... * 1 )
    // 
    for i := 0; i < k; i++ {
        result *= n - i
        result /= i + 1
    }

    return result
}

(Note that we might want to add some additional error checks to the Binomial() function.)

Variations Counting Function in Go

Now we can put everything together into a function to count the number of "Hamming neighbors" - variations on a given DNA string that are a Hamming distance of up to \(d\) away from the original DNA string.

To count the number of Hamming neighbors, we implement the formula above. We leave out the error checks on the parameter values here, for brevity.

The inputs are:

• \(n\) - length of DNA input string
• \(d\) - maximum Hamming distance
• \(c\) - number of codons (4, ATGC)

// Given an input string of DNA of length n,
// a maximum Hamming distance of d,
// and a number of codons c, determine
// the number of Hamming neighbors of
// distance less than or equal to d
// using a combinatorics formula.
func CountHammingNeighbors(n, d, c int) (int, error) {

    // We require the following:
    // n > 0
    // d >= 0
    // c > 0

    // Use combinatorics to calculate number of variations
    nv := 0
    for dd := 0; dd <= d; dd++ {

        // Binomial(n,d) => number of ways we can
        //                  pick codons to edit
        next_term := Binomial(n, dd)

        // (c-1)^d => number of ways that the codons
        //            we picked to edit can be edited
        for j := 0; j < dd; j++ {
            next_term *= (c - 1)
        }
        nv += next_term
    }
    return nv, nil
}

We can run this with a few values of k and d to verify it returns the expected values:

For kmer AAA k = 3:
d = 0, count = 1
d = 1, count = 10
d = 2, count = 37
d = 3, count = 64

for a kmer of length 3, we can compute the first 3 values (1, 10, 37) by hand. The last value, when \(d = k\), is a special case where every base pair in the DNA strand can be changed to any codon. Since there are 4 possible codons, this leads to \(4^k = 2^{2k}\) possibilities.

For \(d = k = 3\), we have \(2^6 = 64\) possible DNA strings.

Now, moving on to \(k=5\):

For kmer AAAAA k = 5:
d = 0, count = 1
d = 1, count = 16
d = 2, count = 106
d = 3, count = 376
d = 4, count = 781
d = 5, count = 1024

We can calculate 1 and 16 by hand, verifying those two numbers. As before, the case of \(k = d = 5\) gives a total of \(4^5 = 2^{10} = 1024\) possible DNA strings.

Next: Recursive Backtracking in Go for Bioinformatics Applications: 2. Generating Variations

Basic Data Structures in Go: Maps

Continuing with our series of blog posts on what we've been learning about Go in the process of solving problems on Rosalind.info, this post will cover how some basic data structures work in Go, and how we used each to solve problems from the Chapter 1 Rosalind problems.

Maps

The simplest way to describe a map is to say it is a structure useful for storing key-value pairs.

Before we walk through what maps look like in Go, let's talk about what a map is (in the data structure sense). And to do that, it's useful to talk a bit about mathematical functions and maps in the mathematical sense.

What is a map

The term "map" is taken from mathematics. A map is just a relationship - a general term, but still useful. Most mathematics courses deal with functions, which are formally defined as maps from one set of numbers onto another, such that one input corresponds to one output.

A (data structure) map, similarly, is a relationship between two sets - a key set, and a value set. Each key corresponds to only one value.

Maps are typically stored under the hood as either a hash map (which does not sort keys and has very fast O(1), or constant time, operations) or a tree map (which sorts keys using a binary tree and has slower O(log N) operations).

Map notation in Go

In Go, map types are denoted map[keyType]valueType.

For example, to map strings to integers, we would use a map of type map[string]int.

We create a map variable by declaring its type:

var m map[string]int

However, this does not allocate any space for the map, and trying to add keys to the map at this point would result in an error.

We need to actually allocate space for the map. In Go, you can allocate space for a map two ways: first, using Go's built-in make() function; and second, by creating and populating the map in one line.

Using make() with maps

To allocate space for a map, you can use the make() function and pass it the map type. This will actually create space in memory for the map, and allow you to add items or look up keys in the map.

var m map[string]int
m = make(map[string]int)

Creating and populating

If you want to create and populate the map in one line, you can specify the type, then have trailing brackets containing the items you want to add:

// Create an empty map
m := map[string]int{}

// Create and populate a map
m := map[string]int{"A":10, "T":15, "G":20, "C":25}

Zero values

One feature of maps that makes them really easy to work with is, if you try and look up a key, and the key does not exist in the map, the map will not raise an exception (which Python does), it will return the zero value of the value type.

For example, the zero value of the int type is 0, so if we create a map like m := map[string]int{"A":10} and we then look up a key that isn't in the map, like m["G"], Go will return 0.

Similarly, the zero value for booleans is false, so you can utilize the zero value behavior to create a set data structure using maps. By creating a map[keyType]bool map, you can use the boolean value to indicate membership of a key in the given set. Then, if you look up keys that do not exist in the map, Go will return the zero value of booleans by default, which will be false.

Easy iterating over maps

It is easy to iterate over maps using the range keyword in Go, which will return the keys (optionally, both keys and values) of the map in a loop:

var m map[string]int{
        "ABC":10, 
        "DEF":20, 
        "GHI":30
}

for k,v := range m {
    fmt.Println("Key:",k," --> Value:",v)
}

Example: Assembling kmer histogram

Here is an example using maps: this function assembles a kmer histogram from a strand of DNA base pairs. To do this, it loops over every codon in the DNA strand and increments a counter in a map. Finally, this histogram map is returned. (This function is useful for determining the most frequent kmer.)

Here's the whole function that helps solve Rosalind problem BA1B. We'll look at it piece by piece below:

// Return the histogram of kmers of length k 
// found in the given input
func KmerHistogram(input string, k int) (map[string]int,error) {

    result := map[string]int{}

    if len(input)<1 {
        err := fmt.Sprintf("Error: input string was not DNA. Only characters ATCG are allowed, you had %s",input)
        return result, errors.New(err)
    }

    // Number of substring overlaps
    overlap := len(input) - k + 1

    // If overlap < 1, we are looking
    // for kmers longer than our input
    if overlap<1 {
        return result,nil
    }

    // Iterate over each position,
    // extract the string,
    // increment the count.
    for i:=0; i<overlap; i++ {
        // Get the kmer of interest
        substr := input[i:i+k]

        // If it doesn't exist, the value is 0
        result[substr] += 1
    }

    return result,nil
}

Function Walkthrough

The first thing you'll notice is the comment style: there is a comment right before each function, which is common practice in Go, because comments that come before a function are picked up by godoc (the Go documentation tool) and turned into documentation.

// Return the histogram of kmers of length k 
// found in the given input
func KmerHistogram(input string, k int) (map[string]int,error) {

Go has internalized the idea that comments are a part of the documentation, so comments don't need to be formatted in any special way (like /// this or /** this */ business) to end up being picked up by godoc.

Next, we create an empty map (kmer strings to integer frequency counters) and stride over the entire input string with a window the size of the kmers we are interested in, adding or incrementing each corresponding counter in the map as we go.

The overlap variable is the number of possible kmers of length k in the entire input string, which requires a bit of algebra to gt the indices right:

    // Number of substring overlaps
    overlap := len(input) - k + 1

The for loop utilizes Go's slice notation to take a slice of the string (which does not require creating or duplicating any string data), and uses the extracted kmer as a key to add or increment the counter:

    // Iterate over each position,
    // extract the string,
    // increment the count.
    for i:=0; i<overlap; i++ {
        // Get the kmer of interest
        substr := input[i:i+k]

        // If it doesn't exist, the value is 0
        result[substr] += 1
    }

This is where the behavior of maps for non-existent keys comes in handy - in Go, if you ask for a key that does not exist in the map, the map will return the zero value of the specified type.

This statement:

        result[substr] += 1

can also be written as:

        result[substr] = result[substr] + 1

So, the first time a kmer (stored in substr) is encountered, the kmer will not exist as a key in the map, and the key lookup on the right hand side will return 0, incrementing the counter to 1 the first time the kmer is encountered.

Subsequent times the kmer is encountered, the value will be found and substituted on the right side so it will be incremented by 1.

Finally, when we return from the function, we can use a Python-like syntax of returning multiple values separated by commas, yet another great feature of Go:

    return result,nil
}

By convention, the return types of functions will include an error type at the very end, so if you had a function named ExampleFunction that returned three integers, the function definition would look like this:

func ExampleFunction() (int, int, int, error) {
    ...    
}

Additionally, we also see from the function's return statement above that we can use the reserved keyword nil to set a variable to a null value, and that the convention is to return nil in place of the error, so that in the calling function we can set up a structure like this:

a, b, c, err := ExampleFunction()

if err != nil {
    err := "Please read this important message"
    return -1,-1,-1,errors.New(err)
}

Summary

Maps are my favorite data structure, so I'm glad that they're easy to use in Go. Some important points about using maps in Go:

• Declaring a variable as a map type does not allocate any space for the map; saying var m map[keyType]valueType and then trying to access keys will cause an exception.

• To allcoate space for a map, use make(map[keyType]valueType) or instantiate with the {} notation, like m := map[keyType]valueType{"asdf":1}

• To ask for a key, use the square bracket notation. To set a value for a key, use the square bracket notation on the left and assign a value on the right: m[my_key] = my_value

• Asking for missing keys will return the zero value of whatever type the map values are.

• Iterating over key-value pairs in a map using a for loop is easy using the built-in range keyword.

Addendum: Check if a Key is in a Map

Because of the default behavior of maps, where they return a zero value for keys that do not exist in the map, it is not immediately obvious how to differentiate between the case where a key is in a map already and has a zero value, versus the case where the key does not yet exist in the map and the zero value is only being returned because the key can't be found.

To resolve this, we can use two Go features: error-checking, and the underscore - a variable that serves as a one-way sink for information, and serves a similar purpose to the underscore in Python.

First, we mentioned above that various operations that can return errors will return an error type as the last return type, along with any other return values. This includes the operation of accessing a key. To assign the value and the error to variables at once:

v, err := m[my_key]

Now, to check if the key exists in the map, we are only concerned with the variable err and we don't really need the variable v. Instead of assigning v to a variable that we never use, and then having the Go compiler complain about it, we can use the underscore as a placeholder:

_, err := m[my_key]

Now, we just add a check for whether err is nil, and voila, we have our check of whether a key is in the map or not:

_, err := m[my_key]

if err != nil {
    // This key is missing from the map

} else {
    // This key already exists in the map

}