Programming Project 1.4: Developing an Algorithm

Most of the discussion of Programming Project 1.4 is applicable to any language. Our focus is on Algorithm 1.3, in which we take a different perspective on the purpose of the algorithm and on the code produced.

Algorithm 1.3, like most of the algorithms in this text, takes a perspective that the algorithm (a) accept[s] an input string from the user, (b) does some computation, and then (c) outputs some values. Because Python provides an interactive computing environment, our model will more often be that an algorithm (a) takes input as a parameter, (b) does some computation, and then (c) returns some values.

One benefit of this alternative approach is that we can more easily chain together algorithms. For example, we can use one function to compute the second DNA strand and then another function to process it.

So, let's look at our new definition of the goals of the algorithm. We're actually going to start with a simpler algorithm, one that, given a nucleotide, returns the paired nucleotide.

Algorithm:	Base Pair
Input:	nucleotide, a string or character representing one of the four nucleotides (A, C, G, or T).
Output:	pair, a string representing the paired nucleotide (T for A, G for C, C for G, or A for T).

This algorithm is fairly straightforward:

• If the nucleotide is A, return T
• If the nucleotide is C, return G
• If the nucleotide is G, return C
• If the nucleotide is T, return A

There is, however, one subtlety. What if the nucleotide is neither A nor C nor G nor T? There are at least two responses. We can ignore that case, since it violates the preconditions of the algorithm (that is, we've only said that the algorithm works if the input is A, C, G, or T). We could also return a special character to indicate something else. We'll do the latter, and return a '*' to indicate a problem.

Okay, what does the Python look like? Remarkably similar to the English we might write, except with a few more symbols.

def base_pair(nucleotide):
    """Convert a nucleotide (represented as a one-character string) to 
       the paired nucleotide."""
    if nucleotide == 'A':
        return 'T'
    elif nucleotide == 'C':
        return 'G'
    elif nucleotide == 'G':
        return 'C'
    elif nucleotide == 'T':
        return 'A'
    else:
        return '*'

Okay, what's going on here?

The def base_pair(nucleotide): is a function definition. It means we're defining a procedure, base_pair, that takes one input, called nucleotide.

The text surrounded by triple quotes is a documentation comment. It allows the programmer to find out about a procedure using Python's built-in help system.

>>> help(base_pair)
base_pair(nucleotide)
    Convert a nucleotide to the paired nucleotide.

The if nucleotide == 'A': is a test. We use the double equal signs for comparison (since a single equal sign is used to assign a value to a variable).

The return 'T' says to end the function, and return a value (either to another function, or to the top-level interpreter).

The following elif (short for else if) commands represent subsequent tests.

Having found a way to deal with individual nucleotides, we're now ready to turn toward the problem of dealing with a whole strand of DNA, rather than just a single nucleotide.

Algorithm:	Opposite Strand
Input:	strand, a list corresponding to a strand of DNA. The elements of strand are only 'A', 'C', 'G', and 'T'.
Output:	result, a list corresponding to a strand of DNA. The elements of result represent the appropriate pairs to the elements of the first. For example, an 'A' in strand will become a 'T' in result.

Now, let's consider how we might make this translation.

• Initialize the result to be the empty list.
• For each base in strand, add the corresponding base to the result using the rules of the Base Pair algorithm.
• Return the result.

Here's one way to implement that algorithm in Python, similar to the one from p. 12 of Exploring Bioinformatics.

def opposite_strand(strand):
    """Given a single DNA strand, reprented as a sequence, compute the 
       opposite strand."""
    result = []
    for nucleotide in strand:
        result.append(base_pair(nucleotide))
    return result

Let's consider the parts of that algorithm.

The def opposite_strand(strand): is a procedure definition. It says that we're defining a procedure (a chunk of code) that takes as input a strand.

The text surrounded by three double quotation marks is a documentation comment. As you've already seen, documentation comments permit us to specify text that appears when someone asks for help about our procedure.

The result = [] initializes a variable called result to the empy sequence.

The for nucleotide in strand: is an instance of a looping structure. It grabs each value from a sequence, in turn, calls it nucleotide, and then does any subsequent indented lines.

So, we use each nucleotide in the strand, find its pair, and then append it to our result sequence.

Finally, the return result returns the sequence we've computed.

If you think this all makes sense (and even if you don't), you're ready to go on to the corresponding lab.