Genome 540 Homework Assignment 1

Due Sunday Jan. 19

Late homework policy is described on course web page.

  1. Download and begin reading Initial sequencing and analysis of the human genome. The Genome International Sequencing Consortium. Nature 409, 860-921 (15 February 2001) . (To print this out, I would recommend the pdf format which corresponds exactly to the printed version, rather than the html format.) For next week, read
  2. Write a program that implements the method described in Lecture 1 to find the longest exactly matching subsequence between two sequences. Specifically, your program should read two input files in "FASTA" format (i.e. having a header line which starts with the character ">" and includes the sequence name, with the sequence itself following on subsequent lines), each of which contains one of the sequences to be compared. As a data check, your program should count the number of bases of each type in each sequence. It should then use the described algorithm to find the longest subsequence present in both sequences. If there are several different perfectly repeated subsequences of the same length, find all of them. When you are looking for these, consider both strands of each sequence (i.e. both the sequence given in the FASTA file, which is sometimes called the 'forward' strand, and its reverse complement, which is sometimes called the 'reverse' strand) simultaneously, so that if the repeated subsequence happens to occur on different strands in the two sequences you will still find it. The easiest way to do this is to create and store in memory a single sequence of length 2N + 2M + 4 (where the two input sequences have lengths N and M bases respectively) constructed by concatenating together the two sequences and their reverse complements (each terminated by a null character to terminate the strings, so that the lexicographic sort does not read across the sequence boundaries), and then create a single list of pointers to each position in this merged sequence, which you then sort. If a longest subsequence is present multiple times in either genome, please report all locations of the subsequence.
  3. Using this program, you should then find the longest exactly matching sequences in the Candidatus Carsonella ruddii genome and the Sorangium cellulosum genome - the smallest (165Kb) and largest (15Mb) known bacterial genomes, respectively. Once you find the longest match(es) using your program, you should then try to figure out what biological feature they correspond to. You can do this by exploring the features tables for both organisms here and here.
  4. To test whether your program is working correctly, run it first on the test example indicated below (with two different bacterial genomes) to see whether you get the right answer. The FASTA files for the test examples can be found here and here. The Genbank files are here and here. In general, FASTA files for these and other bacterial genomes can be found by going to the NCBI web site and following appropriate links. (On the NCBI site, the FASTA files containing the full genome sequences have the suffix .fna). (To find the biological features, look in the 'Genbank format' files for the organisms (which have the suffix .gbk on the NCBI site) and find the annotated 'feature' in each genome that overlaps the matching segment you found. You don't need to write a program to do this -- just can just read the .gbk file on the web site.)
  5. You must turn in your results and your computer program, using the template (file format) described below. Please put everything into ONE file - do not send an archive of files or a tar file. After creating a plain text file (NOT a word processing document file) in this format, compress it (using either Unix compress, or gzip -- if you don't have access to either of these programs let us know), and send it as an attachment to both Phil at phg@u.washington.edu, and Benjamin at bvernot@gmail.com.
Here is an example of a homework file with the fields filled in with the correct answers for a test case, between the "=============" (several items are explained in detail below the template):

=============
Assignment: GS 540 HW1
Name: Benjamin Vernot
Email: bvernot@gmail.com

Fasta 1: NC_000912.fna
>gi|13507739|ref|NC_000912.1| Mycoplasma pneumoniae M129 chromosome, complete genome
A=249211
C=162920
G=163703
T=240560
N=0
*=816394

Fasta 2: NC_004829.fna
>gi|294660180|ref|NC_004829.2| Mycoplasma gallisepticum str. R(low) chromosome, complete genome
A=349623
C=159237
G=159490
T=344450
N=0
*=1012800

Match length: 122
Number of match strings: 2

Match string: GGGGTCTTTCCGTCCCGTTGCGCCTAACGGGTGTCTTCACCCGTACCTGGATTTCACCGAGTCTATAGCCGAGACAGTCAAGAGATGGTTACACCATTCAAGCGGGACGGAATTTACCCGAC
Description: This sequence comes from the 23S ribosomal RNA gene.

Fasta: NC_000912.fna
Position: 122006
Strand: reverse

Fasta: NC_004829.fna
Position: 338335
Strand: reverse

Fasta: NC_004829.fna
Position: 82469
Strand: reverse

Match string: GTCGGGTAAATTCCGTCCCGCTTGAATGGTGTAACCATCTCTTGACTGTCTCGGCTATAGACTCGGTGAAATCCAGGTACGGGTGAAGACACCCGTTAGGCGCAACGGGACGGAAAGACCCC
Description: This sequence comes from the 23S ribosomal RNA gene.

Fasta: NC_000912.fna
Position: 122006
Strand: forward

Fasta: NC_004829.fna
Position: 82469
Strand: forward

Fasta: NC_004829.fna
Position: 338335
Strand: forward

Program:

int main() {
	do_analysis();
	return 0;
}

=============

Details:

Fasta: put the name of the fasta file, along with the first line.

Nucleotide histogram: count the total number of bases ("*") and the number of times each specific base occurs. i.e.:

A=349623
C=159237
G=159490
T=344450
N=0
*=1012800

Description: Put a short identification of the shared DNA subsequence, for instance: "This DNA sequence is the first 20 bases of an RNA polymerase gene"

Position: Put the location of the matching subsequence within the input sequence here. The location should be the index of the DNA base in the sequence that is closest to the beginning of the forward strand. Use a coordinate system starting at 1 rather than 0. For example, if the two chromosomal strands are:

			 5'-ACTGA-3'
			 3'-TGACT-5'
and you found the subsequence TCA on the reverse strand to be the longest match to the other sequence, then the location should be reported as 3. If instead you found CTG on the forward strand, then the location should be reported as 2.