Lab 5 - Gene Prediction!

Finding genes in a genome is not a trivial task. As we have seen in class, a signigficant level of sophistication is required to complete this task. In this laboratory, you will get the opportunity to work directly with genome files and implement a simple gene prediction algorithm. The performance of this algorithm will then be evaluated with the assumption that the output from GLIMMER3 is completely accurate.

As you will soon see, prokaryotic genomes can be treated as any other FASTA file. Since, you already have performed a conversion from FASTA to a Python string. It is therefore not difficult to regard a whole genome as a string.

Please use this MS-Word template to answer the questions for this assignment. Print your solution and submit a hard copy before the following lab.

Here is an overview of the lab:

This sounds like it should hurt, but it won't. Let's make a series of assumption on the structure of a gene:

1. Genes are defined as ORFs.
1. Begin with a ATG
2. End with either TAA, TAG, TGA
3. Both start and stop codon must be in the same frame.
4. START codons are more rare in a random sequence than either STOP codons.
5. Since there can be multiple ATG within a gene, the earliest in-frame ATG codon from a STOP codon is considered to be the real START codon.
2. Randomly occuring ORFs are very short. 

We'll work with the genome of Bacillus anthracis ames although any genome would work as well. For reference, the files below could be downloaded from the Bacillus anthracis ames folder on NCBI's FTP site.

1. GLIMMER3 output
2. Genome in FASTA format.

Save these text files in a new folder called "gene_prediction".

 

I will assume that you saved the genome as NC_003997.fna.

The following code has already been created and saved as the file geneXtract.py.

Download this script and save it in the same "gene_prediction" folder.

Please read through the code and explanation described below. You should understand what each line of code is doing.

Set two application-wide variables:

You will be able to change these experimental parameters later on, but they will do for now.

Read and clean the genome similarly to what was done in the previous lab.

Here is a plain English description of what we are going to do:


It would be great if we could just type this in a Python file and make it happen. However, this description contains ambiguities that would prevent a computer from executing as desired. The second step of designing an algorithm is to break down this description into small steps which are indivually simple to implement using a real programming language.


This translates in Python as:


The most arcane lines are these using the .find() string method. This method returns -1 if no substring are found, this will explain why we want to run the algorithm for as long as start is NOT equal to -1. The find method also begins searching at the index provided as the second parameter. If no second parameter is provided, Python assumes that this value is 0: which is the beginning of the genome.

Make sure you understand the role of each line in the bigger picture. If in doubt, ask questions.

Now, run the geneXtract.py script.

Compare the output of your program with the GLIMMER3 file. You may also find the .fna file useful to refer to. Notice how the algorithm failed to find the START codon of the first gene, although it correctly identified the STOP codon. Can you explain why?

It is now time to generate an outfile which we will use to evaluate how good is the prediction algorithm. The output from GLIMMER can be found in the genome folder as NC_003997.GLIMMER3. Let's have a look at a GLIMMER3 gene entry:

>gi|30260195|ref|NC_003997.3| Bacillus anthracis str. Ames, complete genome
orf00001 407 1747 +2 7.50
...

The first line contains a FASTA style header. Each subsequent line contains the following tokens:


We can reproduce this output with relative ease. By following these steps:

1. Before looping over the genome, open a file in writing mode called prediction.txt.
2. Also before the loop, add a bogus header to the file. 
3. After the loop, close the file.
4. Instead of the print statement found within the loop, use the code below to prepare a string called line and write lineto your output file.

Uncomment all lines of code in the genXtract.py script that activate these steps.

Run the script again and you should now have a newly created "prediction.txt" file that has the same format as the NC_003997.GLIMMER3 file.

For the sake of simplicity, we are disregarding all genes encoded on the complementary strand. We can identify these genes because their coding frame is a negative value. To compare GLIMMER3 and our method, we should remove all negative frame genes from the GLIMMER3 output. Instead of removing, we'll read the GLIMMER3 file and write to a new file all lines which have a positive frame (and thus contain the character +).

Open a new Python file, call it filter.py .

We'll be reading in this file using a different approach: on a line-per-line basis.


Execute this script and there should be a new file in your folder.

Today, you will be deriving the performance of the prediction manually. You may partner to complete the following two tables (because it is tedious):

Let's modify our gene model such that the codon TTG is also accepted as a START codon. Make a copy of your gene prediction program and call the copy geneXtract.alt.py .

The only modification to the code that is needed is to select as start value the smallest index returned for both ATG and TTG. Python has a handy function to compare two number and return the smallest:


In the example above, smallest will be assigned to either a or b, whichever is smaller. In our situation, a will be substituted by a the result of a find method call for the next ATG while b will be the index of the next TTG.



Note that this line will cause a bug near the end of the genome, but we will not worry about this in this laboratory.

Run the alternative script and obtain the gene prediction as before.




Next week, we'll see how we can analyse the results of the predictions to evaluate the performance over a larger span of the genome, as well as the evaluation of the prediction at the nucleotide level (and you will NOT count these by hand!).