Bad variable and function names

I teach a lot of beginner-targeted programming courses, and something that I’ve experimented with recently is trying to introduce the idea of self-documenting code early on in the learning process. I usually start off by talking about the difference between good and bad names for things (mostly functions and variables, though many of the same arguments apply to class names) , and I’ve noticed a few common patterns that tend to crop up in beginners code. I thought it might be useful to lay out these common errors in one place…. most examples are biological in nature, but hopefully they will still be comprehensible to non-biologists.

Single-letter names

OK, we’re writing a program and we need to create a new variable, but we can’t think of a good name….let’s just start with a and work our way through the alphabet. Later on we find that we have a bit of code like this:

a = 'acgatagc'
b = len(a) - 2
d = ""
for e in range(0,f,3):
    g = a[e:e+3]
    h = i.get(g.upper(), 'X')
    d = d + h

Which is well on its way to becoming completely incomprehensible. Sometimes single-letter names are the result of a desire to avoid typing or a worry about running out of space in your text editor. Neither of these are good enough reasons to write code as unreadable as the example above! Rest assured that using longer, more descriptive names will not make your program slow, or make any appreciable difference to the time it takes to type a line. Here’s a sane version:

dna = 'acgatagc'
last_codon_start = len(dna) - 2
protein = ""
for codon_start in range(0,last_codon_start,3):
    codon = dna[codon_start:codon_start+3]
    amino_acid = genetic_code.get(codon.upper(), 'X')
    protein = protein + amino_acid

which might now be interpretable as part of a program that translates DNA sequences.

Sometimes we might want to use non-meaningful variable names to illustrate a very generic bit of code: for example, if we want to demonstrate how to append a value to a list:

a = []
a.append(b)

but for these purposes it’s better to use the well-known set of metasyntatic variables:

foo = []
foo.append(bar)

You can find a handy list of such variable names in several languages at Wikipedia.

There are a couple of situations where single-letter variables do make sense; mostly where there are strong conventions for their use. For example, if we’re writing a program to deal with Cartesian co-ordinates then I won’t be too upset to see variables called x and y (though I might make a case for x_pos and y_pos). Similarly, i and j are the traditional names for variables used as counters in a loop:

for i in range(10):
    for j in range(20):

but remember that the most common use for these variables – to hold an index when iterating over a list – doesn’t often occur in Python because we generally iterate over the elements of lists directly, in which case there’s no excuse for not picking a meaningful variable name.

Naming thing after their type

This is a habit that people are most likely to fall into shortly after having learned of the existence of types, or shortly after having learned about a new type. The logic goes something like this: I’ve just been told that it’s important to remember whether a variable is a string or a number, so I’ll make that fact part of the name. This is not necessarily a terrible idea – in fact there is an entire system of variable naming based on it. Where it becomes a problem is when the type becomes the most important part of the name:

my_number = 20
my_string = "Homo sapiens"
the_list = [1,2,3]
a_file = open('foo.txt')
def my_function(dna):
    ...

This is obviously problematic: it’s generally much more important to know what values are stored in a variable:

minimum_orf_length = 20
species_name = "Homo sapiens"
reading_frames = [1,2,3]
input_file = open('foo.txt')
def translate_dna(dna):
    ...

There’s a more subtle problem with types-as-variable-names – the dynamic nature of Python means that it tends to work best when we worry about the various ways that a variable can be used, rather than its type (see Duck Typing). It’s this magic that allows us, for example, to iterate over lists, strings and files using a single syntax.

Extremely vague names

Often when we create a variable, or start writing a function, we’re not exactly sure what its job is going to be in our program. This is especially true when we first start writing code. Unfortunately, this can lead to some very unhelpful variable names – examples I have seen in the wild include data, input, output, do_stuff(), process_files(), params, object, and true_or_false. If you find yourself using a “placeholder” name like these during the process of coding, it’s a good idea to go back and change them once you’ve figure out what the function or variable is actually doing.

Sequential names

This is a perennial problem in the world of naming, whether we are talking about Python variables or word-processing files (how many people have a folder containing files with names like final_draft.doc, final_draft2.doc, final_draft2.1.doc, final_draft2_update.doc?) Thought of the perfect variable name, but then realized that you’ve already used it? No problem, just stick a “2″ on the end of that bad boy:

dna = 'agtcgnatgc'
dna2 = dna.upper()
dna3 = dna2.replace('N', '')

Hopefully it’s not necessary to point out why this can be confusing when you come back to read the code. We can rescue the above example in a couple of ways. One is to use more descriptive names:

dna = 'agtcgnatgc'
uppercase_dna = dna.upper()
clean_dna = uppercase_dna.replace('N', '')

Another way is to recognize that we’re probably not going to use dna or uppercase_dna in our program, and just do the whole thing in one step:

clean_dna = 'agtcgnatgc'.upper().replace('N', '')

When we find ourselves needing to create a new variable name by sticking a number onto the end of an existing one it’s often a good indication that the code in question should be turned into a function. One of the great thing about encapsulation using functions is that they provide a way for multiple variables with the same name to happily co-exist in a program without interfering with each other (have a look at the concept of a namespace for more on this idea.)

An even worse version of sequential names is….

Re-using names

Thought of the perfect variable name but you’ve already used it? Never mind, just overwrite it! Often this is a symptom of variables names that are too general to begin with:

# store a DNA sequence
sequence = 'ATGC'
...
# now we need to store a protein sequence
sequence = 'LSPV'

Re-using variables to store different things can make a program extremely difficult to understand. Solutions involve either using more descriptive names, which renders the above example fine:

dna_sequence = 'ATGC'
...
protein_sequence = 'LSPV'

or splitting up the code into functions where the same variable name can be re-used without fear of confusion.

Don’t confuse the idea of re-using variables for a different type of data, as in the dna/protein example above, with the idea of changing the value that’s stored in a variable, but keeping the type the same. For example, we often want to update a count, or append a character to a string:

count = count + 1
dna = dna + codon

This is fine, and doesn’t count as re-use, because the variables are still storing the same thing.

Names that only vary by case or punctuation

How many different ways can we write essentially the same variable name?

mydna
myDna
myDNA
Mydna
MyDna
MYDNA
my_dna
my_Dna
my_DNA
My_dna
My_Dna
MY_DNA

Quite apart from the Python style guidelines (all-caps names should be used only for global variables), using more than one of the above in a program will lead to madness…. don’t do it!

Names copied from example code

This is a trap that’s easy to fall into when working from examples in a course or a textbook. Imagine we are looking at a piece of example code that prints the number of each of the four bases in a DNA sequence:

dna = 'ATCGTGACTGTCAT'
for base in ['A', 'T', 'G', 'C']:
    print("count for " + base + " is " + str(dna.count(base)))

Later on, we want to implement the same idea for counting the number of times each amino acid occurs in a protein sequence, so we copy and paste the example code and modify it. We replace the dna sequence with a protein sequence, and replace the bases with amino acid codes:

dna = 'MSRSLLLRFLLFLLLLPPLP'
for base in ['M', 'L', 'F']:
    print("count for " + base + " is " + str(dna.count(base)))

The program works, but the variable names are now very misleading. Using a piece of example code as a starting point for your own programs is an excellent way to learn – but be sure to go back once you’ve finished modifying it and check that the variable names still make sense.

Sign up for occasional updates and useful Python tips.

0

Advanced Python for Biologists

book

Advanced Python for Biologists is now finished and available to buy as either an instant download, or as a physical book. It’s filled with all the cool stuff that didn’t fit into Python for Biologists, including things like:

  • object-oriented Python
  • functional programming
  • comprehensions
  • exceptions
  • testing
  • recursive functions

and lots more useful stuff. Take a look at the the new book by clicking here.

Sign up for occasional updates and useful Python tips.

0

Why readable, documented code is especially important for scientists (and a three-step plan for getting there)

During my most recent teaching engagement I spent some time talking specifically about code readability and documentation. As often happens, presenting these ideas to a roomful of novice programmers helped to crystallize my thoughts on the topic, and made me realize that I’d never written about it – plus I thought that it would be an ideal topic for first post of the new year, since documentation is something that many programmers constantly resolve to do better at!

There’s no shortage of articles and book chapters explaining the general importance of documenting your code – if you learned to program from a book or an online tutorial (such as Python for Biologists) then it will almost certainly have been mentioned. The arguments in favour of documentation are well-rehearsed: it makes it easier for you to work on your own code over a long period of time, it makes it easier for others to contribute fixes and features, it forces you to think about the purpose of each section, etc. In this post, I want to talk about why documentation is particularly important for you – somebody who is using programming for carrying out scientific work. The basis of my argument is that for us as scientists, code serves two important features over and above simply being executed: it acts as both the ultimate supplementary methods, and it’s the way in which you express your original ideas.

Code as supplementary methods

It’s probably fair to say that most users of most programs don’t care about how they actually work, as long as they do what they’re supposed to. When you fire up an image editing program to brighten up a dull digital photo or rotate one where the horizon isn’t straight, you probably aren’t interested in exactly what transformation is being applied to the RGB values, or what trigonometry is being used to straighten the image – you’re only interested in the end result.

Scientific software is different: the end-users are often extremely interested in how the program works internally, since understanding that is a part of understanding the results. And the ultimate way to resolve questions or disagreements about what a program is doing is to examine the source code. This is a great advantage we have when working in bioinformatics. For wet-lab work, there is only so much information you can give in the pages of a journal about how an experiment was carried out. Using supplementary information can help, but even then you’re limited to what the authors thought important enough to write down. For a bioinformatics experiment, however, one can always see exactly where the data came from and what happened to it, providing one has access to the source code. You can read about a piece of bioinformatics software in a journal, listen to a talk on it, discuss it with the authors, but at the end of the day if you still have questions about how it works, you can always go back to the source code.

The vast majority of programmers don’t have to worry about their users wanting to read the source, but we do – so we should make readability and documentation a priority to make sure that it’s as useful as possible.

Code as a way of expressing original ideas

The vast majority of software projects don’t implement any ideas that are particularly original. This isn’t a problem, it’s just a reflection of the fact that many pieces of software do very similar things to other pieces of software, and do them in similar ways. There are fairly standard ways of writing a blog engine, a stock management program, an image manipulation program etc. We could make an argument, therefore, that for those categories of software it’s not super-important that the code is really well-documented, since it’s unlikely to be doing anything surprising, and a reader can probably work out what’s going on in each section by referring to other software that carries out the same task.

Scientific software is different. Yes, we tend to write scripts to carry out tedious everyday tasks like tweaking file formats and correcting sequence headers, but we also use it for implementing entirely new ideas about how to assemble genomes, or how to correct frameshift mutations, or how to pick species for phylogenetic analysis. We’re far more likely than other programmers to write code that does something entirely new. As such our programs (at least the ones that do something interesting) are going to be harder to understand than yet another text editor or chat program.

As programmers, we’re very lucky in that the language we use to implement our original ideas – code – is also an excellent way to communicate them to other researchers. But the usefulness of that language depends on whether we write it in a readable way and document it well.

Three steps to readable, documented code

Documentation is a skill that is learned over the course of a career, but here’s an exercise that I often have my students do. Using a framework like this can make documenting your code less daunting if you’ve no idea where to start.

Step one: make sure your variable and function names are meaningful

Programmers are fond of talking about self-documenting code – i.e. code that doesn’t require external documentation to be understood. A large part of this is using meaningful variable names. Examples of bad variable and function/method names include:

  • Single-letter names e.g. a, b, f (with the exception of variable names that follow common conventions such as x and y for co-ordinates or i for an index)
  • Names that describe the type of data rather than the contents e.g. my_list, dict
  • Names that are extremely generic e.g. process_file(), do_stuff(), my_data
  • Names that come in multiples e.g. file1, file2
  • Names that are excessively shortened e.g. gen_ref_seq_uc
  • Multiple names that are only distinguished by case or punctuation e.g. input_file and inputfile, DNA_seq and dna_seq
  • Names that are misspelled – the computer does not care about spelling but your readers might

Go through your code and look for any instances of the above, and replace them with good names. Good variable names tell us the job of the variable or function. This is also a good opportunity to replace so-called magic numbers – constants that appear in the code with no explanation – with meaningful variable names e.g. 64 might be replaced by number_of_codons.

Example: we want to define two variables which hold the DNA sequence for a contig and a frame, then pass them to a method which will carry out protein translation and store the result. Here’s how not to do it, even though the code is perfectly valid Python:

a = 2
b = 'ATGCGATTGGA'
c = do_stuff(a, b)

This is much better:

frame = 2
contig_dna_seq = 'ATGCGATTGGA'
contig_protein_seq = translate(frame, contig_dna_seq)

Step two: write brief comments explaining the reasoning behind particularly important or complex statements

For most programs, it’s probably true to say that the complexity lies in a very small proportion of the code. There tends to be a lot of straightforward code concerned with parsing command-line options, opening files, getting user input, etc. The same applies to functions and methods: there are likely many statements that do things like unpacking tuples, iterating over lists, and concatenating strings. These lines of code, if you’ve followed step one above, are self-documenting – they don’t require any additional commentary to understand, so there’s no need to write comments for them.

This allows you to concentrate your documentation efforts on the few lines of code that are harder to understand – those whose purpose is not clear, or which are inherently difficult to understand. Here’s one such example – this is the start of a function for processing a DNA sequence codon-by-codon (e.g. for producing a protein translation):

for codon_start in range(0, len(dna)-2, 3):
codon_stop = codon_start+3
codon = dna[codon_start:codon_stop]
    ...

The first line is not trivial to understand, so we want to write a comment explaining it. Here’s an example of how not to do it:

# iterate over numbers from zero to the length of
# the dna sequence minus two in steps of three
for codon_start in range(0, len(dna)-2, 3):
...

The reason that this is a bad comment is that it simply restates what the code does – it doesn’t tell us why. Reading the comment leaves us no better off in knowing why the last start position is the length of the DNA sequence minus two. This is much better:

# get the start position for each codon
# the final codon starts two bases before the end of the sequence
# so we don't get an incomplete codon if the length isn't a multiple of three
for codon_start in range(0, len(dna)-2, 3):
...

Now we can see from reading the comment that the reason for the -2 is to ensure that we don’t end up processing a codon which is only one or two bases long in the event that there are incomplete codons at the end of the DNA sequence.

Go through your code and look for lines whose function isn’t obvious just from reading them, and add explanations

Step three: add docstrings to your functions/methods/classes/modules

Functions and methods are the way that we break up our code into discrete, logical units, so it makes sense that we should also document them as discrete, logical units. Everything in this section also applies to methods, classes and modules, but it keep things readable I’ll just refer to functions below.

Python has a very straightforward convention for documenting functions: we add a triple-quoted string at the start of the function which holds the documentation e.g.

def get_at_content(dna):
  """return the AT content of a DNA string.
     The string must be in upper case.
     The AT content is returned as a float"""
  length = len(dna)
  a_count = dna.count('A')
  t_count = dna.count('T')
  at_content = float(a_count + t_count) / length
  return at_content

This triple-quoted line is called a docstring. The advantage of including function documentation in this way as opposed to in a comment is that, because it uses a standard format, the docstring can be extracted automatically. This allows us to do useful things like automatically generate API documentation from docstrings, or provide interactive help when running the Python interpreter in a shell (take a look at the chapter on testing and documentation in Advanced Python for Biologists for an in-depth look at how this works).

There are various different conventions for writing docstrings. As a rule, useful docstrings need to describe the order and types of the function arguments and the description and type of the return value. It’s also helpful to mention any restrictions on the argument (for instance, as above, that the DNA string must be in upper case). The example above is written in a very unstructured way, but because triple-quoted strings can span multiple lines, we could also adopt a more structured approach:

def get_at_content(dna):
  """return the AT content of a DNA string.

     Arguments: a string containing a DNA sequence.
                The string must be in upper case.

     Returns: the AT content as a float"""
  ...

If you think it’s helpful, you can also give examples of how to use the function in the docstring. Notice that we’re not saying anything in the docstring about how the function works. The whole point of encapsulating code into functions is that we can change the implementation without worrying about how it will affect the calling code!

Summary

These three steps represent the minimum amount of work that you should do on any code that you plan on keeping around for more than a few weeks, or that you plan on showing to anybody else. As always, if you have questions or suggestions, leave a comment.

Sign up for occasional updates and useful Python tips.

1

What you have in common with the Wright brothers

Warning: vast historical oversimplification below in pursuit of a point :-)

Famously, the Wright brothers built and flew the first aircraft capable of sustained, powered flight in 1903. Looking at the famous photos with eyes used to seeing modern aircraft, it looks pretty airworthy:

220px-1904WrightFlyer

 

There were plenty of other people working on heavier-than-air flying machines around that time, many with much more money and far more resources. So what was the key to the Wright brothers’ success? Did they invent a new type of engine? A new type of wing? Not really – their greatest invention was this:

wright_tunnel

This unprepossessing-looking box is a wind tunnel, which the Wright brothers – realizing that it was far too time-consuming to test wing designs by building them full scale – used to test their aeronautical designs using models. The innovation that prompted their break-through was not an improvement to aircraft, but an improvement in the process for designing aircraft. By using a wind tunnel, they were simply able to make their mistakes faster than anyone else, and to learn from them. Others had to learn by building, and crashing, full-size aircraft.

This is far from an original observation, but I think it has some connection with programming. The story of the Wright brothers illustrates the power of rapid iterative improvement – their approach would probably be called “agile” if it were being used today. The difference between the Wright brothers and their contemporary rivals mirrors one that I often see between the different approaches to writing code I see being used by my students.

On the one hand, you have people who favour small, incremental improvements when writing a program or a function, testing each bit of code as soon as possible and uncovering bugs and mistakes early. Students who program in this way end up with programs and functions that resemble the Wright Flyer pictured above: crude and primitive, perhaps, but certainly fit-for-purpose and relatively unlikely to result in broken bones.

On the other hand, you have people who try to write an entire program or function all in one go, never testing any bit of it until the whole thing is written. Students who program in this way end up with programs and function that resemble other products of early aviation:

images

 

As the picture above attests, this is a recipe for pain.

 

Sign up for occasional updates and useful Python tips.

0

Improving database web interfaces with userscripts

Those of us who spend most of our day working on the command line have generally got into the habit of writing small, simple scripts to solve everyday problems and annoyances. The web-browser equivalent of these everyday problem-solving scripts are userscripts – small snippets of Javascript which run on specific web pages. In this post I’m going to show a quick example of how we can use a userscript to add a missing feature to the NCBI taxonomy browser.

Two quick notes before we get started…..

What is a userscript? A userscript is just a piece of Javascript plus a few bits of metadata stored as comments.

How do I use a userscript? You have to install an extension for your web browser – follow the instructions here.

The problem

When I’m in a phylogenetic frame of mind, I often find myself browsing the NCBI taxonomy. This is the database that stores taxonomic information and relationships for all organisms with sequence data in GenBank and we can use it to view a page for any taxonomic group – here’s a screenshot of a bit of the page for Endopterygota (click to embiggen):

Selection_006

This page shows us the various taxonomic groups that belong to Endopterygota laid out in hierarchically. Clicking on the name of one of these groups will take us to the appropriate page, so we can navigate around the tree of life quite easily. We can also display some extra information on this page – checking the “Nucleotide” box above the tree and then clicking “Display” will cause the number of nucleotide records belonging to each group to be displayed after the name:

Selection_007

This is pretty useful – we can go straight to the list of nucleotide records for a given group by clicking the number, but we can also just use the numbers to survey the distribution of sequence data for the groups we’re looking at. For example, in the above view, there are lots of sequence for butterflies and moths, and beetles. The trouble is that the view presented above isn’t a very intuitive way to look at the relative numbers – reading the counts and comparing requires a fair amount of mental overhead. What would be great is if there were an option on the NCBI website to display the number of nucleotide records for each group visually – say, as a bar whose width corresponds to the number of records:

Selection_008

The NCBI website doesn’t have such a feature but, as you can probably guess from the above screenshot, we are going to add it ourselves using a userscript.

The Solution

Before we start coding, we can break the problem down into a few steps. First, we need to get a list of all the counts that appear on the original web page. Then, we need to figure out what the largest count is so that we can scale the bars to a sensible size. Finally, we need to replace each count with a bar of the appropriate width.

Getting the list of counts is pretty easy – if we look at the source HTML for the page we can see that each of the nucleotide counts is an a element with the title attribute set to “Nucleotide”. We can use the JQuery library to grab the list of count elements and store it as nuc_counts:

// get list of nucleotide counts
var nuc_counts = $('[title="Nucleotide"]');

Calculating the maximum count is a little bit trickier. We need to take each count, strip out all the commas, turn the count into an integer, then grab the maximum value from the list of integers. I won’t spend time here going into the craziness required to get the maximum value from an array in Javascript: suffice it to say that we’ll use JQuery’s map function to turn our list of string counts into a list of integers, then find the maximum and store it in a variable called max_nuc_count:

// calcaulate maximum nucleotide count
max_nuc_count = Math.max.apply(Math, $.map(nuc_counts, function(x,i){return parseInt(x.text.replace(/,/g,""))}))

Now for the main body of the script. We’ll iterate over our array of count elements, and for each one use JQuery to construct a new element to replace it. The new element will be a div, and we’ll need to set its width to a value that reflects the original count. To do this, we’ll take the count, multiply it by five hundred, then divide the result by the maximum count that we calculated earlier – in other words, we’ll scale all the bars so that the widest one is five hundred pixels wide. The only other tricky bit is making sure that the bar gets displayed inline with the taxon name, rather than on a line of its own – to do this, we set the “display” attribute of the bar to “inline” or “inline-block”:

// for each count element...
for (var i=0; i<nuc_counts.length; i++){
    var count_element = nuc_counts[i];

    // remove the commas from the number and turn it into an integer
    var count = parseInt(count_element.text.replace(/,/g,""));

    // use jquery to create a new div element which will be the bar representing the nucleotide record count
    bar = $('<div>&nbsp;</div>')	// the div needs to contain a non-breaking space; if it is completely empty then it will not be displayed
    	.css('margin-bottom', 2)	// add a tiny space at the bottom so that there's a little gap between bars
    	.css('display', 'inline-block')	// force the div to display as an inline element so that it can share a line with the taxon name
    	.css('background-color', 'RoyalBlue') // pick a nice colour for the bar
    	.css('width', (count * 500) / max_nuc_count);	// calculate the width for the bar, scaled to the max

    // replace the original count element with the new bar
    $(count_element).replaceWith(bar);
}

So far so good: this gives us a nicely scaled set of bars and makes sure that the widest bar (i.e. the one at the top, which is the sum of all the others) fits on the screen. We could easily make the bar scale bigger or smaller by changing the 500 in the above code to something else – we could even take into account the width of the browser window if we wanted.

Finally, let’s add a couple of finishing touches. There are two things missing from the above solution: firstly, there’s no way to see the actual numbers, and there’s no way to click through to the list of records themselves. We can solve the second problem by creating an anchor element to wrap around the bar, with the target url copied from the original count. And we can solve the first problem by giving the anchor a “title” attribute which contains the original count, so that when we hover the mouse cursor over a given bar, it will display the exact number of nucleotide records. JQuery does most of the hard work here:

// get list of nucleotide counts
var nuc_counts = $('[title="Nucleotide"]');

// calcaulate maximum nucleotide count
max_nuc_count = Math.max.apply(Math, $.map(nuc_counts, function(x,i){return parseInt(x.text.replace(/,/g,""))}))

// for each count element...
for (var i=0; i<nuc_counts.length; i++){
    var count_element = nuc_counts[i];
    
    // remove the commas from the number and turn it into an integer
    var count = parseInt(count_element.text.replace(/,/g,""));
    
    // use jquery to create a new anchor element which will link to the nucleotide records
    anchor = $('<a></a>')
    	.attr('href', count_element.href)	// use the original count as a tooltip
    	.attr('title', count_element.text); // grap the nucleotide search url from the original element
    
    // use jquery to create a new div element which will be the bar representing the nucleotide record count
    bar = $('<div>&nbsp;</div>')	// the div needs to contain a non-breaking space; if it is completely empty then it will not be displayed
    	.css('margin-bottom', 2)	// add a tiny space at the bottom so that there's a little gap between bars
    	.css('display', 'inline-block')	// force the div to display as an inline element so that it can share a line with the taxon name
    	.css('background-color', 'RoyalBlue') // pick a nice colour for the bar
    	.css('width', (count * 500) / max_nuc_count);	// calculate the width for the bar, scaled to the max
    
    // put the bar inside the anchor so that you can click on
    anchor.append(bar);	
    
    // replace the original count element with the new anchor/bar
    $(count_element).replaceWith(anchor);
}

And there we have it. To turn this into a userscript, all we have to do is add a set of specially-formatted comments at the top which can be parsed by whichever browser extension we want to use. In particular, we need to specify which web pages the script should run on using a regular expression (the @match line below). Here’s the script in full:

// ==UserScript==
// @name       NCBI Taxonomy nucleotide record count barchart
// @namespace  http://pythonforbiologists.com/
// @version    0.1
// @description replace nucleotide record counts in NCBI taxonomy with bars, see http://pythonforbiologists.com/index.php/adding-features-ncbi-taxonomy/
// @match      http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi*
// @copyright  2012+, You
// ==/UserScript==

// @require  http://ajax.googleapis.com/ajax/libs/jquery/1.8.3/jquery.min.js

// get list of nucleotide counts
var nuc_counts = $('[title="Nucleotide"]');

// calcaulate maximum nucleotide count
max_nuc_count = Math.max.apply(Math, $.map(nuc_counts, function(x,i){return parseInt(x.text.replace(/,/g,""))}))

// for each count element...
for (var i=0; i<nuc_counts.length; i++){
    var count_element = nuc_counts[i];
    
    // remove the commas from the number and turn it into an integer
    var count = parseInt(count_element.text.replace(/,/g,""));
    
    // use jquery to create a new anchor element which will link to the nucleotide records
    anchor = $('<a></a>')
    	.attr('href', count_element.href)	// use the original count as a tooltip
    	.attr('title', count_element.text); // grap the nucleotide search url from the original element
    
    // use jquery to create a new div element which will be the bar representing the nucleotide record count
    bar = $('<div>&nbsp;</div>')	// the div needs to contain a non-breaking space; if it is completely empty then it will not be displayed
    	.css('margin-bottom', 2)	// add a tiny space at the bottom so that there's a little gap between bars
    	.css('display', 'inline-block')	// force the div to display as an inline element so that it can share a line with the taxon name
    	.css('background-color', 'RoyalBlue') // pick a nice colour for the bar
    	.css('width', (count * 500) / max_nuc_count);	// calculate the width for the bar, scaled to the max
    
    // put the bar inside the anchor so that you can click on
    anchor.append(bar);	
    
    // replace the original count element with the new anchor/bar
    $(count_element).replaceWith(anchor);
}

If you want to install this extension and try it out, I’ve added it to the userscripts.org repository – you should be able to install it by going here, once you have installed the browser extension. If you come up with any improvements to the code, or have any suggestions for other database web interface fixes or features, shout out in the comments!

Sign up for occasional updates and useful Python tips.

0

The world’s worst genome assembler in six lines of Python

So, after I posted my new business cards the other day I got a comment to the effect that I should have made one with an aligner. That got me thinking about the biggest thing that I could conceivably fit on a business card, if I didn’t care about readability. So I decided that I could probably fit an incredibly bad sequence assembly program on one. Just for fun, I wrote the whole thing in a purely functional style – a total of six lambda expressions.

lp (4)

Warning: the rest of this post contains discussions of horrible code, rampant abuse of Python features, and the complete opposite of all good programming practise. Read on with caution!

The API is quite straightforward: give the ah() expression a list of DNA sequences, and it will return a consensus string. Here’s an ultra-short example:

reads = ['TCCCAGTGAACCCA', 'TTCCGTGCGGTTAAG', 'GTCCCAGTGAACCCACAA', 'TGAACCCACAAAACG', 'ACCCACAAAACGTGA', 'GAACCCACAAAACGTGA', 'TCCGTGCGGTTAAGC', 'TGAACCCACAAA', 'CCGTGCGGTTAAGCGTGA', 'TGACAGTCCCAGTGAA', 'AACCCACAAAACGTGA', 'AGTGAACCCACAAAACGT', 'GTTAAGCGTGA', 'CCGTGCGGTTAAGCGTGA', 'AGCGTGACAGT', 'TGCGGTTAAGCG', 'ACAAAACGTGATG', 'ACAGTCCCAGTGAACC', 'TAAGCGTGACAGTCCCA', 'TCGAATTCCGT', 'TTCTCGAATTCCGTGCG', 'ACAAAACGTG', 'CCACAAAACGTG', 'TGCGGTTAAG', 'GAACCCACAAAACGTGA', 'TCTCGAATTCC', 'ATTCCGTGCGGTTAA', 'ACCCACAAAAC', 'CGTGCGGTTAAGCGTGA', 'CCAGTGAACCCACAA', 'TGCGGTTAAGCGTG', 'CCCACAAAACG', 'TCTCGAATTC', 'AATTCCGTGCGGTT', 'ACAGTCCCAGTGA', 'GTCCCAGTGAACCCA', 'TGAACCCACAAA', 'CCCACAAAACGTG', 'TCCCAGTGAACCCACA', 'CTCGAATTCCGTGCG']
print(ah(reads))
# TGCGGACAAAACGTGTGAACGTGAGGTTAAGCGTGACAGTCCCAGTGAACCCACAAAACGT

And here’s a local alignment to the sequence from which the fake reads were generated, just to prove that it’s actually doing something slightly better than picking the longest read:

 6 GAATTCCGTGCGGTTAAGCGTGACAGTCCCAGTGAACCCACAAAACGT     53
   |||   ||||.|||||||||||||||||||||||||||||||||||||
17 GAA---CGTGAGGTTAAGCGTGACAGTCCCAGTGAACCCACAAAACGT     61

The longest read is 20 bases, and the local match region is 48 bases.

How it works

Let’s start by putting each expression on its own line:

cm = lambda d,ds: max([m(d,e) for e in ds if e != d])
m = lambda d,e: max([(s(d,e,o),o,e,d) for o in range(1-len(e),len(d))])
s = lambda d,e,o: sum([1 for p in range(max(0-o,0), min([len(e)-o, len(e), len(d)-o])) if e[p] == d[p+o]])
con = lambda x,o,s,c : c[0:max(0,o)] + s +  c[len(s)+o:]
a = lambda s, o : con(*cm(s, o)) if len(o) == 1 else a(con(*cm(s, o)), [ y for y in o if y != cm(s, o)[2]])
ah = lambda d : a(d[0],d[1:])

Next we’ll convert each expression to an equivalent function:

def cm(d,ds):
    return max([m(d,e) for e in ds if e != d])

def m(d,e):
    return max([(s(d,e,o),o,e,d) for o in range(1-len(e),len(d))])

def s(d,e,o):
    return sum([1 for p in range(max(0-o,0), min([len(e)-o, len(e), len(d)-o])) if e[p] == d[p+o]])

def con(x,o,s,c):
    return c[0:max(0,o)] + s +  c[len(s)+o:]

def a(s, o):
    return con(*cm(s, o)) if len(o) == 1 else a(con(*cm(s, o)), [ y for y in o if y != cm(s, o)[2]])

def ah(d):
    return a(d[0],d[1:])

And change the names to be slightly more descriptive (this makes some lines extremely long, so you’ll have to scroll right to read them), and rearrange them slightly:

# given two sequences and an offset, count the number of matching bases
def score(sequence1,sequence2,offset):
    return sum([1 for position in range(max(0-offset,0), min([len(sequence2)-offset, len(sequence2), len(sequence1)-offset])) if sequence2[position] == sequence1[position+offset]])

# given two sequences, find the offset which gives the best score
def find_best_offset(sequence1,sequence2):
    return max([(score(sequence1,sequence2,offset),offset,sequence2,sequence1) for offset in range(1-len(sequence2),len(sequence1))])

# given a single sequence and a collection of others, find the other sequence with the best match score
def find_best_match(sequence,others):
    return max([find_best_offset(sequence,sequence2) for sequence2 in others if sequence2 != sequence])

# given two sequences and an offset, calculate the consensus
def consensus(score,offset,sequence1,sequence2):
    return sequence2[0:max(0,offset)] + sequence1 +  sequence2[len(sequence1)+offset:]

# given a sequence and collection of others, return the complete consensus using recursion
def assemble(sequence, others):
    return consensus(*find_best_match(sequence, others)) if len(others) == 1 else assemble(consensus(*find_best_match(sequence, others)), [ y for y in others if y != find_best_match(sequence, others)[2]])

# given a collection of sequences, call assemble() to start the recursion
def assemble_helper(dnas):
    return assemble(dnas[0],dnas[1:])

Now we can look at each function in more detail.

score() is the basic scoring function. It takes two sequences and an offset, and returns the number of characters that match. It’s the simplest possible function for scoring an ungapped alignment between two sequences. It works by calculating the first and last positions of the overlapping region relative to sequence2, then counts up the number of positions for which the base in sequence2 is the same as the base in sequence1 at that position plus the offset. To get everything in one expression, it uses a list comprehension to build a list of 1′s for each matching position, then sums the list. Here it is written out a bit more conventionally. The only complicated thing going on here is the calculation of the start and stop positions.

def score(sequence1,sequence2,offset):
    start_of_overlap = max(0-offset,0)
    end_of_overlap = min([len(sequence2)-offset, len(sequence2), len(sequence1)-offset])
    total_score = 0
    for position in range(start_of_overlap, end_of_overlap):
        if sequence2[position] == sequence1[position+offset]:
            total_score = total_score + 1
    return total_score

find_best_offset() is the function that tries to maximize the score for a pair of sequences by trying every possible offset. It works by first calculating the range of possible offsets, then using a list comprehension to build a list of tuples, one tuple for each possible offset. Each tuple contains the score, the offset, and the two sequences – this slightly weird way of storing the results is necessary so that the information can be passed to the other functions, as we’ll see in a minute. To find the single best offset, we take advantage of the fact that in Python, sorting a list of tuples sorts them by their first element. Since the first element of each of our tuples is the score, if we simply ask for the max() of the list we get the tuple with the highest first element i.e. the one representing the best score. Here’s the sensible version:

def find_best_offset(sequence1,sequence2):
    lowest_offset = 1-len(sequence2)
    highest_offset = len(sequence1)
    all_offsets = []
    for offset in range(lowest_offset,highest_offset):
        # add the 4-tuple for this offset
        all_offsets.append((score(sequence1,sequence2,offset),offset,sequence2,sequence1))
    return max(all_offsets)

find_best_match() is probably the most straightforward function of the bunch. Given a single sequence and a list of other sequences, it finds the other sequence that has the best match by calling find_best_offset() for each of them in turn. It uses the same tuple-sorting trick as before to figure out which match is the best:

def find_best_match(sequence,others):
    all_matches = []
    for sequence2 in others:
        if sequence2 != sequence:
            all_matches.append(find_best_offset(sequence,sequence2))
    return max(all_matches)

The consensus() function gave me quite a bit of trouble. Its job is to take two sequences plus a given offset, and return the consensus sequence of the two. Of course, it doesn’t do anything like what we normally mean by consensus – it simply concatenates the relevant bits of the two sequences to make a longer one. The logic behind how it works is a little bit hard to follow. We construct the consensus sequence by taking the full length of sequence one, and sticking any left-hand overhang from sequence two on the left end and any right-hand overhang from sequence two on the right end. In other words, you should read the return line as “return any bits of sequence2 that stick out to the left, followed by the whole of sequence1, followed by any bits of sequence two that stick out to the right”. For most overlapping pairs of sequences, either the first or last bit of the returned string will be zero-length, which is why the thing works as a single expression in the compact version.

def consensus(score,offset,sequence1,sequence2):
    sequence2_left_overhang = sequence2[0:max(0,offset)]
    sequence2_right_overhang = sequence2[len(sequence1)+offset:]
    return sequence2_left_overhang + sequence1 + sequence2_right_overhang

The assemble() function is probably the most complicated (and certainly the most inefficient). I cheated a little bit to get it onto a single line by using the ternary operator “x if y else z”. It’s a recursive function that takes a single sequence and a collection of other sequences. It finds the best match for the sequence among the others and calculates the consensus of the sequence and the best-matching other. If that’s the only member of others (i.e. the others list has just one element) it simply returns the consensus. If the others list has more than one element, it removes the best-matching one and calls itself recursively with the newly-built consensus as the single sequence. Here it is expanded:

def assemble(sequence, others):
    # remember, best_matching_other is a 4-tuple
    best_matching_other = find_best_match(sequence, others)
    # the * expands the elements of the tuple so we can use them as arguments to consensus()
    consensus_sequence = consensus(*best_matching_other)
    if len(others) == 1:
        return consensus_sequence
    else:
        # get the second element of the best_matching_other tuple, which is the sequence
        best_matching_sequence = best_matching_other[2]
        others.remove(best_matching_sequence)
        return assemble(consensus_sequence, others)

assemble_helper() is, as the name suggests, just a helper function which kicks off the recursion by calling assemble() with the first element as the single sequence and the remainder of the elements as the list of other sequences.

Let’s sum up the algorithm then (described iteratively, even though it’s written recursively). To assemble a list of N DNA sequences, we take the first sequence, and find the member of the remaining N-1 sequences which has the best match. We remove this best-matching member from the list (leaving N-2 sequences) and calculate the consensus of these two sequences. We then append the newly-built consensus onto the end of the list (bringing the sequence count back up to N-1), then go back to the start and begin again. Hopefully it’s clear that, since the number of sequences in the list decreases by one in each iteration, we will eventually end up with a list of just a single sequence, which is our result.

A discussion of the performance of this algorithm is both beyond the scope of this post, and entirely unnecessary. Suffice it to say that it has horrible performance in terms of both computation time and results! At around 430 characters I think that the compact version is pushing the limits of what can fit on a business card in 12 point text. If anyone can think of a way to squeeze a few characters, please let me know in the comments!

Sign up for occasional updates and useful Python tips.

0

New business cards!

I’ve been meaning for a while to get round to making some business cards to hand out to folks who ask me about learning to program. Normally I just tell people to google “python for biologists” and they’ll end up in the right place, but it would be nice to have a physical reminder to give out. At first I though about having some USB memory stick business cards made  - there are some really cool ones that are the shape of a normal business card but can fold in half to reveal a set of USB contacts. Unfortunately they’re way expensive, and the minimum order is far more than I need.

Next I thought about making a “cheat sheet” style business card – the type with contact information on the front and some useful quick-reference information (e.g. a list of regular expression characters) on the back. I guess the idea would be that the recipient is more likely to hang onto the card if it has useful information on it. But I couldn’t think of anything that would fit in well with my website – after all, the emphasis of pythonforbiologists is on learning to program, not simply the practice of programming itself.

Finally, I had an idea; I would put a tiny biology-themed programming exercise on the back of each of my business cards, along with a link to a web page giving the solution.

cards1

This would hopefully mean that when somebody gets hold of one of my cards they can see straight away what kind of material and training I provide, and can head over to the website for more information. I wrote five different nano-exercises on five different biological topics:

  • parsing FASTQ file format
  • counting the number of occurrences of short motifs in DNA sequences
  • calculating AT content using a sliding window
  • generating the reverse complement of a DNA sequence
  • calculating restriction fragment lengths

Fitting the sample code onto the business cards was quite difficult. I wanted to make sure that the code would be readable and not too hard to understand – I even found room for a few comments – but it also had to be very concise. I only had about ten lines to work with, so I had to use very short variable names.

cards2

You can see images of all the reverse sides at this link.

After I’d designed the code samples and exercises I wrote web pages for each of the solutions. I decided to put the exercise description and the link to the solution pages on the front of the card, as I’d used up all the room on the back with the code samples.

cards3

I tried to make each solution page interesting to read. As well as giving an answer to the exercise, I included extra material about useful bits of the Python language that some people don’t know about. For example, in the solution page to the FASTQ parser exercise I talked about generator functions, and in the sliding window exercise solution I talked about higher-order functions.

You can browse all of the exercises along with links to their solution pages here. Comments appreciated!

Sign up for occasional updates and useful Python tips.

3

The role of instructors in programming training

I’ve been spending a bit of time recently arranging to run some instructor-led training courses early next year (see my training page if this is something you’re interested in), which has got me thinking about the role of the instructor in teaching/learning programming. This is a pretty important question – essentially, why is an instructor-led course better than self-directed learning? – so I decided to write down my thoughts on the matter.

Roadblocks to learning

I want to address the question in a kind of roundabout way, by first talking about a phenomenon I have tentatively called roadblocks. These are the moments that occur whenever you’re learning any new skill in which there’s something that you don’t quite understand, which prevents you from making further progress. When you’re learning a new skill and you hit a roadblock, it’s not necessarily a bad thing – indeed, the point when you overcome a roadblock by correcting your understanding is probably the instant we would point to if asked exactly when “learning” occurs. The nature of programming, though, means that these roadblocks tend to come thick and fast, especially for beginners, and are particularly hard to get around.

Programming is unintuitive

It’s very easy for those of us who have been programming for a long time to forget just how unintuitive programming is. When you’re first learning to write code, there’s no obvious reason why variables or loops or dictionaries behave the way they do – these things seem profoundly arbitrary. Ironically, it’s only once you dig deeper into programming, and understand something about how these things are implemented, that their behaviour seems to flow predictably from the way they work. Once you know about stack frames, it’s easy to see why scoping works the way it does – but that knowledge usually arrives too late to save you from having to struggle with scoping issues earlier in your education. This unintuitive behaviour means that roadblocks arise more frequently than in other fields of learning.

Programming is relentlessly progressive

The practise of programming revolves around building up small, simple pieces of functionality (like statements) into bigger, more complicated ones (like functions, objects and entire programmes). The process of learning to program follows the same pattern – you start by learning the most basic, atomic blocks (variable assignment, simple statements) then use this knowledge to bootstrap your ability to use more complicated things (loops, functions, etc.). The practical upshot of this is that you have to understand each of the simpler building blocks in order to make sense of the more complicated ones. In other words, when a roadblock comes along, you can’t continue to make progress by simply going around it and moving onto a different topic.

By way of analogy, imagine you’ve never cooked a meal in your life, and you decide to learn cooking from scratch using an online tutorial. You happily work your way through the various chapters with titles like Sauces and Pickles and Vegetable Soups until you come to the section on Bread Making. At this point, you run into a roadblock – none of your breads will rise, and you have no idea what you are doing wrong1 . No matter – you can just skip the bread section for now, and move onto the chapter on Rice and Pasta.

Contrast this with the situation where you’re learning programming from scratch. You make it through Variable Assignment and Printing Strings, and are making good progress until you encounter a mental roadblock in the chapter on Loops. Somehow you just can’t get your head around the way that the loop variables acquire different values in each loop iteration. Just like in the cooking example, you decide to skip the troublesome section and work on something else for a while, so you move onto the next chapter, Processing Files. Unfortunately, the very first example involves using a loop to parse each line of an input file, and you can’t understand the example because you don’t understand loops. You try again, and skip forward to the chapter on Dictionaries, but again, half of the examples in this chapter use loops to either construct or iterate over dictionaries. Your forward progress is stalled until you can go back and really get to grips with the concept of loops.

Because programming is progressive in a way that cooking is not (and this analogy is in no way meant to belittle cooking as a skill!), the roadblocks that beginners encounter are harder to get around.

It’s hard to ask the right questions

This is another aspect of programming that experience tends to render invisible: when you encounter a roadblock in programming and need to ask for help, very often it’s difficult to know how to phrase the question. The lack of understanding that causes the student to need help in the first place also ensures that they’re unlikely to know what question to ask, or the right way to phrase it. I’m certainly not suggesting that this problem is unique to programming – it’s found in pretty much every technical field – but the highly abstract nature of the things that we talk about (“method calls”, “return values”, “function pointers”) make it particularly tricky. This difficulty in communication explains why, even though the internet is bursting with forums, message boards and mailing lists populated by helpful people who are happy to assist novice programmers, it can take a long time to pin down the root cause of a student’s misunderstanding.

It’s hard to stay motivated

There are two main reasons why students find themselves on my courses – either they want to use programming to solve a problem, or they’ve been told to attend the course by some higher authority (an employer, a PhD supervisor). Overwhelmingly, it’s the ones who have a concrete problem to solve that tend to make better progress, not because of any intrinsic ability, but because the problem provides the motivation for them to persist with learning a difficult skill. Especially in the early stages, learning to program can seem a relatively thankless task, where the only payoff from successfully understanding a tricky new concept is the prospect of moving onto the next, even trickier one.

Of course, later on in the learning process it usually becomes clear how mindblowingly useful programming is (and I purposefully structure my courses to get students to this point sooner rather than later). Nevertheless, one of the biggest problems many students have when learning programming is simply running out of steam and becoming demoralized – a process that is usually triggered by encountering yet another roadblock.

Getting over roadblocks

The point that I am trying to make under the headings above is not simply that programming is hard, but rather that it’s hard in a particular way that is amenable to being solved by the presence of an instructor. The chief role of the instructor, as I see it, is to get students over these roadblocks as quickly and painlessly as possible. To illustrate the value of this, let’s consider a typical-case scenario facing the self-taught programmer…..

Imagine that you have set aside a week from your busy schedule to get started with learning to program, something that you’ve been meaning to do for ages. You sit down at your desk on Monday morning with your chosen learn-to-program book in hand, and start working your way through the exercises. Give that you’re fairly computer-savvy – you know how to use a text editor and a command line – you should have no problem getting a good grounding in the basics by the end of the week.

Just before lunchtime, you run into a roadblock – some concept or example that you can’t seem to figure out. You’ve obviously misunderstood something, because even when you look at the example solution to the exercise, it doesn’t make sense. You decide to take a break for lunch. When you get back from lunch you still don’t understand the example, and resolve to read the chapter again from the start to see if you’ve missed something important. When this doesn’t help, you decide it’s time to ask for help. You post a question on the mailing list for your language of choice, explaining your problem. Then you alt-tab over to Facebook and kill some time while you wait for a reply……

Sometime near the end of the day you get a response. Excellent! it’s from an experienced programmer who’s prepared to help you work through the problem. They’ve emailed you back with a couple of clarifying questions that will help to figure out exactly which bit of the code you don’t understand. Unfortunately, they live on the other side of the world, so you will only be able to communicate during the brief period when you’re both awake. By the time you go home, you’re feeling a bit demotivated; you had planned on getting through at least two chapters per day but you’re still stuck on the first one, and you have a feeling that there’ll be many more roadblock moments to come.

Now, admittedly this is a bit of a gloomy picture – things will not always be this bad! There may be experienced colleagues you can talk to, or you may be able to get over your roadblock with the help of a second tutorial that explains the concept in a slightly different way, etc. But the overall pattern will, I think, be familiar to anyone (myself included) who describes themselves as a self-taught programmer.

Contrast this with what might happen in an instructor-led course. Just as before, you encounter a roadblock in one of the exercises and you can’t understand why your code isn’t working. After puzzling over it for a couple of minutes, you stick your hand up and the instructor comes over. Because the instructor has taught this material many times in the past, they’ve probably seen this exact problem before, and rather than just fixing the code for you, they can use this experience to quickly figure out the cause of your confusion by asking a couple of questions. They can then write a couple of lines of code illustrating the problem that you’re having and explaining how to solve it, while simultaneously clearing up the original source of your confusion, while you watch and ask questions in real-time.

In this way, five minutes after you encountered the roadblock you’ve already overcome it and can move on to the next section and keep making progress. Rather than feeling demotivated that you’ve wasted a bunch of time, instead you feel like you understand the material better for having struggled with it, and are increasingly confident that this programming business might actually turn out to be quite tractable.

I have exaggerated the two scenarios above to make the point, but the central idea remains: the main job of the instructor is to ensure that when a student encounters a roadblock, they overcome it rapidly and don’t simply give up2 .

Having read this far, if you think that instructor-led training could be useful to your organization, get in touch.


  1. Let’s imagine, for the purposes of this analogy, that you are using an expired batch of yeast. 

  2. Of course, there’s lot of other stuff that instructors do: they choose which content to teach and in what order, create learning material, tailor examples to the audience, etc. 

Sign up for occasional updates and useful Python tips.

0

How to count non-DNA bases in a sequence using Python

I noticed recently that two particular questions are popping up quite regularly in my search logs: “how to count non-DNA bases in a sequence” and “how to tell if a sequence contains DNA” (presumably as opposed to protein). It struck me that the second question is really a special case of the first – once we have a way to count the number of DNA bases in a sequence, we can simply apply a rule that if more than 80% (or any other number we choose) of bases in a sequence are A,T,G or C, then it is probably DNA.

Let’s start with the simplest thing that we think will work – we’ll simply count the number of A, T, G and C characters in a sequence, then divide by the length and multiply by 100 to get a percentage. For this example I’m using a DNA sequence that has three non-ATGC characters: one each of N, Y and R. I’ve included the division fix at the start of the code in case you want to run this on Python 2:

from __future__ import division
seq = "ACTNGTGCTYGATRGTAGC"
dna_count = seq.count("A") + seq.count("T") + seq.count("G") + seq.count("C")
dna_fraction = dna_count / len(seq)
print(dna_fraction * 100)

The output from this bit of code shows that it’s working as expected:

84.2105263158

However, in some circumstances, we might want to allow characters other than A,T,G and C in our DNA sequences. Take a look at this table showing the set of standard IUPAC ambiguity codes:

Nucleotide Code:  Base:
----------------  -----
A.................Adenine
C.................Cytosine
G.................Guanine
T (or U)..........Thymine (or Uracil)
R.................A or G
Y.................C or T
S.................G or C
W.................A or T
K.................G or T
M.................A or C
B.................C or G or T
D.................A or G or T
H.................A or C or T
V.................A or C or G
N.................any base
. or -............gap

Depending on which subset of these we want to allow, we might want to count as many as sixteen different characters. Rather than cram sixteen different calls to count() into one line, it’s probably better to loop through the allowed characters and build up the count one at a time. Here’s a bit of code to do that, using a list to define the set of allowed characters. For this example I’m allowing the four standard bases plus purines (R) and pyrimidines (Y):

from __future__ import division
seq = "ACTNGTGCTYGATRGTAGC"
allowed_bases = ["A", "T", "G", "C", "Y", "R"]
total_dna_bases = 0
for base in allowed_bases:
    total_dna_bases = total_dna_bases + seq.count(base)
dna_fraction = total_dna_bases / len(seq)
print(dna_fraction * 100)

As expected, the answer is higher than in our first example because we are now counting the R and Y as DNA bases:

94.7368421053

This seems like a perfect bit of code to turn into a function. We’ll make the DNA sequence and the list of allowed bases into function arguments, and use a sensible default of counting just ATGC characters.

def count_dna(seq, allowed_bases=['A','T','G','C']):
    seq = seq.upper()
    total_dna_bases = 0
    for base in allowed_bases:
        total_dna_bases = total_dna_bases + seq.count(base.upper())
    dna_fraction = total_dna_bases / len(seq)
    return(dna_fraction * 100)

Notice how we’ve changed both the input sequence and the allowed bases to upper case, to make sure that the function will work regardless of the case of the inputs. Here are a few quick tests:

print(count_dna("ACTRGATCYGATCGANTCGATG"))
print(count_dna("ACTRGATCYGATCGANTCGATG", ['A','T','C','G','N']))
print(count_dna("actgratcygtganctttgacg"))

along with their output:

86.3636363636
90.9090909091
86.3636363636

Having written this function, it’s pretty straightforward to define a function to test if a sequence is DNA. To make the function as flexible as possible, we’ll assign sensible defaults to both the allowed bases and the minimum percentage of bases that must match. We’ll pass the input sequence and the list of allowed bases through to the count_dna function, and then compare the result of that call to the minimum. Here’s the function along with a couple of lines to test it:

def is_dna(seq, allowed_bases=['A','T','G','C'], minimum=80):
    return count_dna(seq, allowed_bases) > minimum

print(is_dna("ACTRGATCYGATCGANTCGATG"))
print(is_dna("ACTRGATCYGATCGANTCGATG", minimum=90))
print(is_dna("ACTRGATCYGATCGANTCGATG", minimum=90, allowed_bases=['A','T','G','C','R','Y']))

As you can see, the function is very concise – we simply ask whether the percentage of DNA bases returned by our earlier function is greater than the minimum, and return the result. As the output shows, we can make the test more stringent by increasing the minimum, or more lenient by allowing some ambiguous bases:

True
False
True

Another, much more concise way to write the counting function would be to use a list comprehension to select just the characters that are in some group:

total_count = len(1)

see the chapter on list comprehensions in Advanced Python for Biologists for an in-depth look at similar examples.

Sign up for occasional updates and useful Python tips.

3

Powered by WordPress. Designed by Woo Themes