Coding for the wet lab

Large datasets are becoming the norm in modern biology, but I still see wet-lab students and postdocs who are stymied by anything that can’t be accomplished in Excel. Should every grad student learn how to write code? Probably not. But those who know even the basics of programming will find that their lives are oh so much easier. Not only will they be able to automate analyses to turn a daunting (or impossible) task into an easy afternoon’s work, they’ll also be able to “speak the language”, which will give them a leg up in describing problems and brainstorming solutions with dedicated computational biologists.

There are some nice pieces about this, but let me provide a few examples that should be more concrete to those involved in genome editing.

• Let’s say you’re doing a CRISPRi screen. You’ve done some sequencing, removed constant regions, and now you have a file containing a huge number of protospacers. These should all map to guides and associated genes, which are listed in another file with the format <gene> <guide sequence>. You can write a very simple program to map guide sequences to genes and count the number of times each one appears.
• Or maybe you’ve just edited a gene and you’re interested in the distribution of indels around the cut site. It’s not difficult to record the frequency of each insertion and deletion from the millions of reads in a next gen sequencing experiment.
• Even more simple: You want to target a whole bunch of genes in a pathway, so plan to make 50 guides targeting 25 distinct genomic regions. You can easily automate the design all of the primers needed for the T7E1 assays.

Of course, you must also be aware of the necessity for good statistics. The ability to automate tasks is no replacement for a good computational biologist or biostatistician, who can analyze whether measurements are actually significant. This is a good example of knowing enough to get you into trouble; be aware of when to seek out an expert! But a few hours’ work on a program can save days of trying to process data with inappropriate tools (e.g. Excel).

Personally, I think python is the way to go for most simple programming tasks. It’s reasonably fast, pretty easy to learn, deals well with text, has some good biology-oriented modules (e.g. Biopython), and has a growing user base (making it easy to find help). Python certainly isn’t the fastest, but the focus here is on writing quick one-offs to get a job done rather than making an ultra-fast production tool. General programming knowledge is probably sufficient for most tasks, and here are a few resources to get started:

• Learn Python the Hard Way (the best resource and a misnomer title – really “Learn Python by Doing It”)
• CodeAcademy Python (a very ingenious hands-on system)
• Coursera Programming for Everybody 
• Google’s Python Class
• QB3 Berkeley Python Course

To illustrate how simple a program can be, here’s a script to accomplish the CRISPRi (or any guide library) guide-gene mapping/counting tasks described above. It’s obsessively documented, to try to show how this can be quite straightforward. It is also written for ease of understanding by a beginner rather than speed or conciseness. Writing this took about 20 minutes.

#!/usr/bin/env python
import sys # load up a module that lets us read arguments on the command line
ref_fname = sys.argv[1] # first argument on the command line is the file containing reference guides in the library
expt_fname = sys.argv[2] # second argument on the command line is the file containing sequences of guides measured in the experiment
refguides = {} # set up a dictionary to hold reference guide-gene mapping. Key is guide sequence, value is gene
with open(ref_fname, mode="r") as f: # open the reference file
	for line in f: # read through each line
		gene, guide = line.rstrip().split(" ") # split out the gene name and the guide sequence from the line, assign them to variables
		refguides[guide] = gene # put guide sequence and gene name into the refguides dictionary
expt_genes = {} # set up dictionary to hold counts of genes measured during the experiment. Keys are genes, values are number of times observed.
unmapped_guides = {} # set up dictionary to hold counts of genes observed during experiment that don't map to genes. Keys are guide sequences, values are number of times observed.
with open( expt_fname, mode="r") as f: # open the file containing the guides sequenced during the experiment
	for line in f: # read through each line
		guide = line.rstrip() # assign the observed guide sequence to a variable
		if guide in refguides.keys(): # is this sequence in the reference dictionary?
			gene = refguides[guide] # get the associated gene from the reference dictionary
			if gene in expt_genes.keys(): # have we already seen this sequence?
				expt_genes[gene] += 1 # if so, increment
			else:
				expt_genes[gene] = 1 # if not, set
		else:
			# similar to the if/then increment/set above, but using a try/except block
			try: # we've seen this sequence before
				unmapped_guides[guide] += 1
			except KeyError: # we haven't yet seen this sequence
				unmapped_guides[guide] = 1
# output
print "COUNTED GENES"
for gene, count in expt_genes.items(): # list number of times we've seen guides that map to a gene
	print gene, count
print "UNMAPPED GUIDES"
for guide, count in unmapped_guides.items(): # list number of times we've seen guides that don't map to a gene
	print guide, count