You will implement software to assemble a genome from synthetic reads. We supply Python code snippets that you might use or adapt in your solutions, but you don't have to. Read all the hints, highlighted in bold.
10% of the points for this question
Download the reads:
http://www.cs.jhu.edu/~langmea/resources/f2014_hw4_reads.fa
All the reads come from the same synthetic genome and each is 100 nt long. For simplicity, these reads don't have any quality values.
The following Python code will download the data to a file called reads.fa in the current directory. (Caveat: I don't think the code below works in Python 3. Sorry about that. Go here for details on how to fix: http://python-future.org/compatible_idioms.html#urllib-module.)
In [1]:
import urllib
# Download the file containing the reads to "reads.fa" in current directory
urllib.urlretrieve("http://www.cs.jhu.edu/~langmea/resources/f2014_hw4_reads.fa", "reads.fa")
# Following line is so we can see the first few lines of the reads file
# from within IPython -- don't paste this into your Python code
! head reads.fa
Or you can download them manually from your browser. Or you can use wget or curl from the command line.
This code will help you parse the FASTA:
In [2]:
def parse_fasta(fh):
''' Parse FASTA into a dictionary '''
fa = {}
name = None
# Part 1: compile list of lines for each sequence
for ln in fh:
if ln[0] == '>': # new sequence
name = ln[1:].split()[0]
fa[name] = []
else:
# append nucleotides to current sequence
fa[name].append(ln.rstrip())
# Part 2: join lists into strings
for name, nuc_list in fa.iteritems():
fa[name] = ''.join(nuc_list) # join into one long string
return fa
40% of the points for this question
Goal: Write a file containing each read's best buddy to the right. Let's define that.
For each read $A$, find the other read $B$ that has the longest suffix/prefix match with $A$, i.e. a suffix of $A$ matches a prefix of $B$. $B$ is $A$'s best buddy to the right. However, if there is a tie, or if the longest suffix/prefix match is less than 40 nucleotides long, then $A$ has no best buddy to the right. For each read, your program should output either (a) nothing, if there is no best buddy to the right, or (b) a single, space-separated line with the IDs of $A$ and $B$ and the length of the overlap, like this:
0255/2 2065/1 88
This indicates a 88 bp suffix of the read with ID 0255/2 is a prefix of the read with ID 2065/1. Because of how we defined best buddy, it also means no other read besides 2065/1 has a prefix of 88+ bp that is also a suffix of read 0255/2. A corrolary of this is that a particular read ID should appear in the first column of your program's output at most once. Also, since we require the overlap to be at least 40 bases long, no number less than 40 should every appear in the last column.
Notes:
Hint 1: the following function groups reads such that you can avoid comparing every read to every other read when looking for suffix/prefix matches.
In [3]:
def make_kmer_table(seqs, k):
''' Given dictionary (e.g. output of parse_fasta) and integer k,
return a dictionary that maps each k-mer to the set of names
of reads containing the k-mer. '''
table = {} # maps k-mer to set of names of reads containing k-mer
for name, seq in seqs.iteritems():
for i in range(0, len(seq) - k + 1):
kmer = seq[i:i+k]
if kmer not in table:
table[kmer] = set()
table[kmer].add(name)
return table
Hint 2: here's a function for finding suffix/prefix matches; we saw this in class:
In [4]:
def suffixPrefixMatch(str1, str2, min_overlap):
''' Returns length of longest suffix of str1 that is prefix of
str2, as long as that suffix is at least as long as min_overlap. '''
if len(str2) < min_overlap: return 0
str2_prefix = str2[:min_overlap]
str1_pos = -1
while True:
str1_pos = str1.find(str2_prefix, str1_pos + 1)
if str1_pos == -1: return 0
str1_suffix = str1[str1_pos:]
if str2.startswith(str1_suffix): return len(str1_suffix)
50% of the points for this question
Goal: Write a program that takes the output of the overlap program from part 1 and creates uniquely assemblable contigs (unitigs), using the best buddy algorithm described below.
We currently know each read's best buddy to the right. I'll abbreviate this as bbr. We did not attempt to compute each read's best buddy to the left (bbl), but we can infer this. Consider the following output:
A B 60
E A 40
C B 70
D C 40
$A$'s bbr is $B$. But $B$'s bbl is $C$, not $A$! Your program should form unitigs by joining together two reads $X$ and $Y$ if they are mutual best buddies. $X$ and $Y$ are mutual best buddies if $X$'s bbr is $Y$ and $Y$'s bbl is $X$, or vice versa. In this example, we would join $D$, $C$, and $B$ into a single unitig (and in that order), and would join reads $E$ and $A$ into a single unitig (also in that order).
Your program's output should consist of several entries like the following, with one entry per unitig:
START UNITIG 1 D
C 40
B 70
END UNITIG 1
START UNITIG 2 E
A 40
END UNITIG 2
The first entry represents a unitig with ID 1 consisting of 3 reads. The first (leftmost) read is D. The second read, C, has a 40 nt prefix that is a suffix of the previous read (D). The third (rightmost) read in the contig (B) has a 70 bp prefix that is a suffix of the previous read (C).
Each read should be contained in exactly one unitig. The order of unitigs in the file is not important, but the unitig IDs should be integers and assigned in ascending order.
Hint: the correct solution consists of exactly 4 unitigs.
This part is optional. You can submit your solution if you like. No extra credit will be awarded.
Goal: Assemble the genome! Report the sequence of the original genome as a FASTA file.
This requires that you compare the unitigs to each other, think about what order they must go in, and then put them together accordingly. Submit your solution as a single FASTA file containing a single sequence named "solution". The FASTA file should be "wrapped" so that no line has more than 60 characters. You can use the following Python code to write out your answer.
In [5]:
import sys
def write_solution(genome, per_line=60, out=sys.stdout):
offset = 0
out.write('>solution\n')
while offset < len(genome):
nchars = min(len(genome) - offset, per_line)
line = genome[offset:offset+nchars]
offset += nchars
out.write(line + '\n')
Here is an example of what your output should look like. Note how the sequence is spread over many lines.
In [6]:
import random
random.seed(5234)
write_solution(''.join([random.choice('ACGT') for _ in xrange(500)]))
Hint 1: the correct genome is 7959 nucleotides long
Hint 2: to learn how the unitigs should go together, try overlapping them with each other