Handle exception

As I mentioned in last blog, there is a problem unsolved in the program. That is , each genome must have unique primers. If not, the program is endless.
There exceptions have been considered:
1. the file is empty
2.the size of genome is too small to suggest a primer of length 20
3.one genome A is subset of another genome B, there is no way to find a unique primer for genome A.

The program was tested to handle successfully these three exceptions.
test the program on the  three strains of Bacillus anthracis:

Bacillus anthracis str.A0248.fa : ATAAGCGAATAAGAATGTAG

Bacillus anthracis str.CDC684.fa : TGAAGTTGCTGAGAAGAGAT

Bacillus anthracis str.Sterne.fa : CATCAGTTACAAAAAGAATG

real 168.97

user 167.91

sys 0.16

So far,the problems mentioned before have been solved. The next work is to write a scientific report for this project.

The problem of checking the primer is unique has been solved

I found another way to check if the primer is unique.


for k in range(0,len(genome)):
    while True:
        mer=substring(genome[k])

        if len(check(mer,genome))>1:
            continue
        elif check(mer,genome)[0]!=k:
            continue
        else:
            break

    primer.append(mer)

This new method to check uniqueness of the primer is quite different from the previous one.
I was always thinking in a way to search the substring (primer sequence) in other genome sequences and delete it if it is also present in other genomes. The flow were like this:
1. suggest a primer from genome[i]
2. store it in the primer list as primer[i]
3. check if primer[i] is present in other genomes except genome[i]
4. delete it if yes and go back to 1

while in the new method, the flow control were like this:
1.suggest a primer from genome[i]
2.assign the primer to a string(dose not store it in the primer list immediately)
3.search the string in the genome list including genome[i] and return the found index in a new list
4.if the length of index list is larger than 2, go back to 1
5.if the length of index list is 1, check if the only element is equal to i. if no, go back to 1. if yes, assign the primer to the primer list.

While there is problem in the new method, all the genomes must have unique perimers. If any of them does not, the program is endless. The solution could be that the program will check all genomes before suggesting the unique primer.

Mistakes were found in the second round of test

My first test ran without error is because there was an error in my script. Since k+1 was written as i+1, loop2 which was supposed to check the uniqueness of the primer did not run at all.

for k in range(0,len(genome)): #this loop1 is to suggest a primer for genome[k] and store in primer[].
    primer.append(substring(genome[k]))
    for j in range(k+1,len(genome)): # this loop2 is to check if the primer is unique to genome[j]
        if genome[j].find(primer[-1])!=-1:
            del primer[-1] #if primer[-1] is also present in other genomes, remove it!

But after I corrected it, another error came up.

Traceback (most recent call last):

  File "./primer.py", line 34, in <module>

    if genome[j].find(primer[-1])!=-1:

IndexError: list index out of range

Besides that, I found loop2 was wrong to check if the primer is unique. It always delete the last element of the list primer[] each time the primer was found in other genome. I will work on the problem.

first draft of the program and test

Writing the script for the program is finished. Now start to test it.

1. Control group. First, 3 testing files are made specially for the program. File 1 is empty, file 2 contains a sequence (length<20), file 3 contains a short sequence(20<length<100).
these 3 files are given as input to the program

./primer.py empty.fa test1.fa test2.fa
empty.fa : your genome file is empty or the size of genome is smaller than 20
test1.fa : your genome file is empty or the size of genome is smaller than 20
test2.fa : the unique primer is AGAGAGAGCTTTTTTTTTTT

2. Experimental group. Three strains of Bacillus anthracis, the bacteria causing anthrax.

./primer.py Bacillus\ anthracis\ str.Sterne.fa Bacillus\ anthracis\ str.CDC684.fa Bacillus\ anthracis\ str.A0248.fa 
Bacillus anthracis str.Sterne.fa : the unique primer is AGCTTCACTAGTTTCGATCA
Bacillus anthracis str.CDC684.fa : the unique primer is ACGACAGCGTGGTTCATACG
Bacillus anthracis str.A0248.fa : the unique primer is ATAATCCAGTTTTGTACTAA

3. display information about the source used running the program

Control        E-group
real 0.02       0.16
user 0.01       0.12
sys 0.00         0.03

4. testing repeatedly the program using the same input files will give different results
The program only report one unique primer to each genome, so if there are more than one primers  unique to one genome, the primer which will be reported is randomly selected.

start project: finding unique primers

One way to identify organisms and their specific strains quickly is to construct primers that are specific to a strain and an organism.

The goal of the project is to write a program that suggest regions of length 20 which are unique to the input genomes. It means a primer for genome A will not be present in genome B or C.

construct project directory
-DD2397 project
      -data
      -scripts
      -articles
      -results(project notebook, a link to project blog)