Constructing a De Bruijn Graph

In this problem we are asked to find the  adjacency list corresponding to the De Bruijn graph constructed from a set of short reads and their reverse complements. To do this was actually a lot easier than I initially thought. After a bit of research on De Bruijn graphs and a bit of doodling, I figured out the following three steps to complete the task:

1. Make the reverse complements and combine with the initial reads.
2. Add all unique reads (whether they are reverse complements or the original reads) to a set.
3. For each of the reads in the set, print the prefix followed by the suffix.

I set about to write a program that performs these three steps and the following is what I managed to come up with:

from Bio.Seq import Seq

from Bio.Alphabet import IUPAC

def rev_comp(read):

    seq = Seq(read, IUPAC.unambiguous_dna)

    rev_seq = seq.reverse_complement()

    return(str(rev_seq))

def mk_set(seqs):

    S = set()

    for i in seqs:

        if i not in S:

            S.add(i)

        C = rev_comp(i)

        if C not in S:

            S.add(C)

    return(S)

data = []

with open('rosalind_dbru.txt','r') as f:

    for line in f:

        data.append(line.strip())

s = mk_set(data)

for j in s:

    print("(" + j[:-1] + ", " + j[1:] + ")")

The code can also be found on my GitHub along with the sample data and problem description.

Error Correction in Reads

In this problem we are looking at the errors that occur during genome sequencing. We are given a set of reads of equal length, some of which contain an error (one nt is exchanged). The reads are either of the following:

1. The read was correctly sequenced and appears in the dataset at least twice (possibly as a reverse complement)
2. The read is incorrect, it appears in the dataset exactly once, and its Hamming distance is 1 with respect to exactly one correct read in the dataset (or its reverse complement)

From this dataset we are expected to find the incorrect reads and correct them.

Sample Dataset

>Rosalind_52

TCATC

>Rosalind_44

TTCAT

>Rosalind_68

TCATC

>Rosalind_28

TGAAA

>Rosalind_95

GAGGA

>Rosalind_66

TTTCA

>Rosalind_33

ATCAA

>Rosalind_21

TTGAT

>Rosalind_18

TTTCC

Expected Output

TTCAT->TTGAT

GAGGA->GATGA

TTTCC->TTTCA

I really enjoyed working on this problem and I managed to solve it pretty quickly. A problem I had though was that I would get the correct result for the sample dataset, but when I tried with a 'real' dataset I got way too many reads in my answer, something like 4 times as many as the input reads. I looked at my output and could conclude that a lot of the reads occurred multiple times. I figured something was wrong in the part of my program where I build the list of correct reads, and sure enough, I managed to find the error and fix it (highlighted below). The following code is what I finally ended up with:

from Bio import SeqIO                                                 

from Bio.Seq import Seq                                               

from Bio.Alphabet import generic_dna                                  

reads = []                                                            

handle = open('sampledata.fasta', 'r')                                

for record in SeqIO.parse(handle, 'fasta'):                           

    reads.append(str(record.seq))                                     

handle.close()                                                        

right = []                                                            

wrong = []                                                            

for i, j in enumerate(reads):                                         

    read = Seq(j, generic_dna)                                        

    rev_read = read.reverse_complement()                              

    for k in range(i + 1, len(reads)):                                

        if read == reads[k] or rev_read == reads[k]:                  

            if read not in right and rev_read not in right:           

                right.append(str(read))                               

                right.append(str(rev_read))                           

for l in reads:                                                       

    if l not in right:                                                

        wrong.append(l)                                               

for incorrect in wrong:                                               

    for correct in right:                                             

        hamming = 0                                                   

        for nt1, nt2 in zip(incorrect, correct):                      

            if nt1 != nt2:                                            

                hamming += 1                                          

                if hamming > 2:                                       

                    break                                             

        if hamming == 1:                                              

            with open('answer.txt', 'a') as textfile:                 

                print(incorrect, '->', correct, sep='', file=textfile)

Genome Assembly as Shortest Superstring

Before reading any further in this post, I feel that I must warn you: my solution to this problem is very inefficient, not to say very ugly and probably hard to follow for anyone who isn't me. I know, great programming practice right here! I am still a bit proud of the solution though, because even though it might be the most hideous thing you have seen this month (year?) I managed to come up with it almost entirely without looking things up! And, hey, at least it managed to finish within the five minute limit! So here it goes... Please try not to laugh too hard at me...

Problem Description:
For a set of sequences in FASTA format, we are asked to find the shortest superstring that they make up. The criteria on the sequences is that there exists a unique way to reconstruct the superstring from them by putting together pairs of reads that overlap by more than half their length.

Sample Dataset:
>Rosalind_56
ATTAGACCTG
>Rosalind_57
CCTGCCGGAA
>Rosalind_58
AGACCTGCCG
>Rosalind_59

GCCGGAATAC

Expected output:
ATTAGACCTGCCGGAATAC

Final code:
from Bio import SeqIO
reads = []
handle = open('sampledata.fasta', 'r')
for record in SeqIO.parse(handle, 'fasta'):
    reads.append(str(record.seq))
handle.close()

overlaps = []
overlapping = []
for i in range(len(reads)):
    curr_read = reads[i]
    for j in range(len(curr_read) // 2, len(curr_read)):
        curr_suffix = curr_read[-(j + 1):]
        for k in range(len(reads)):
            curr_comp = reads[k]
            for l in range(len(curr_comp) // 2, len(curr_comp)):
                curr_prefix = curr_comp[:l]
                if curr_suffix == curr_prefix:
                    overlaps.append(k)
                    overlapping.append([len(curr_suffix), i, k])

s = set(overlaps)
first_read = ''
count = len(overlapping)
for m in range(len(overlapping)):
    suf_index = overlapping[m][1]
    if suf_index not in s:           #find first read and initialise new str
        first_read = suf_index
        new_str = reads[overlapping[m][1]] + reads[overlapping[m][2]][
            overlapping[m][0]:]
        count -= 1
        pref_index = overlapping[m][2]
        while count > 0:                       #when the first read is found, add 
            for n in range(len(overlapping)):  #the remaining in the correct order
                if pref_index == overlapping[n][1]:
                    new_str += reads[overlapping[n][2]][overlapping[n][0]:]
                    pref_index = overlapping[n][2]
                    count -= 1

print(new_str)



Yeah, so... Sorry about that.