Lab 7, May 16

In Lab 7, you will continue to work on your toy BLAST program. If you didn't get everything working in the Lab 6 part of the that program, continue working on that as well as the new parts. If you absolutely cannot get the Lab 6 part finished and working by the time it is due, and hence can't continue of with Lab 7, let us know, and we will give you some parts of the code for Lab 6. But it is best if you do it all yourself.

In Lab 7, you will also learn about the Prosite data base of regular expressions.

If you have not already, please also read the third notes on Perl, although no excercises will be assigned from them. Please let me know if you find any errors or ambiguities in those notes.

What you need to turn in are your answers to the questions and exercises below. For Perl programs use script, or some other means, to print out the program and show how it runs on data.

Continuation of the toy BLAST program

First, when a database string matches a substring in the query string by w characters, say, where k < w, then the same w-length substring will be found multiple times by the previous program, and if w > t, the threshold for reporting, it will get reported many times; that is not desireable. Real BLAST has ways to avoid *finding* the same substring multiple times, but we will only make sure that the substgring is not *reported* multiple times. We can do it using a hash called stringhash, like this: Whenever BLAST finds a reportable substring in a database, starting at position $i, say, in the database string, it looks to see if $stringhash{$i} is defined. If so, it does not report the string again. Otherwise it assigns the string to $stringhash{$i} and does report the string. Implement this change in the program.

We would like to process strings that are more than a single line long. So in the file each string will be held in consecutive lines, with strings seperated by blank lines. That is like saying that each string is a paragraph instead of just a single line. To read in a whole paragraph, put the line

$/ = "";

somewhere in the program before the reading begins. Read about this on page 102 of Johnson or some other Perl reference.

Read kmer4.pl and get it running, and try to understand how it works. That will be in preparation for our final change to the toy BLAST program, in lab 8.

Prosite:

http://www.expasy.ch/prosite/

Go to the prosite site and familiarize yourself with the various services available there, and scan the list of other pattern-based utilities.

Be sure to study completely the section called "The optimal way to develop patterns". I am not sure where it is, but it may be in the documents page somewhere. That material will be on the final exam!

We discussed in class that the regular expression that characterises the granin family of proteins is a false positive for BRCA1. Prosite has now been changed so that the search with BRCA1 does not return a granin protein.

Exercise:

• Go to the Prosite homepage and follow the link for the ScanProsite tool. ScanProsite allows you to perform a comparison of your sequence vs. all of the regular expressions in the Prosite database. A match will tell you that your protein is a likely candidate for a particular feature or family. Read the ScanProsite Manual.
• We will be using the SwissProt sequence for BRCA1 as our query sequence. Its accession number is P38398. To search against all the Prosite patterns, leave the window called "Motifs to Scan for" empty. There is an option on the query form for excluding motifs that occur with high probability. First do not use that option. How many hits are found? Explain what the score means - read the documentation sufficiently to understand it. Find the Prosite pattern, regular expression, for the first match reported. Copy that pattern and put it aside.
  Now go back and check the option to exclude motifs that occur with high probability. Do the search again. How many hits are found this time? Explore each hit to see what other information is given about the hits.
• A different way to use ScanProsite is essentially in reverse, i.e., to specify a regular expression pattern and then find all the sequences in the database that have a substring that matches that pattern. Use the pattern you copied earlier to search for all the proteins that have a substring that matches that pattern. Report what you found.
  There is an option to randomize the database in several ways. Find out from the documentation exactly what ``shuffle" does. Run the search on the shuffled database (this may take some time). If it is taking too much time, try the reverse option. Report your results compared to using the original database. Explain.
• Another tool available at Prosite is PRATT, which can create patterns found in a set of sequences. Download the file packed globins and get this tool working for these sequences, if you can. There may be too many. In that case, get it working on a subset of the sequences. You will have to enter the sequences in FASTA format. You should read the Prosite documentation or manual, or look on the web, to learn FASTA format requirements. I tried PRATT for a few minutes but couldn't get it working - but I am sure you will be able to. Report your results.
  

  Extra credit project: Some people have suggested that rather than doing alignment, a better way to measure the biological similarity of two sequences S and S' is to count the number of kmers that they have in common, divided by the total length of S and S'. In more detail, if a kmer is found in S t times, and in S' t' times, then that kmer contributes Min[t,t'] to the total count of common kmers.

  Write a program that asks the user for the value of k, and then finds all the kmers in S and the number of times each one occurs in S, and finds all the kmers in S' and the number of times each occcurs in S'. The program then computes the ``similarity" of S and S' as above.

  Download the file packed globins

  Extend your program to take in a set of sequences and compute the above similarity for each pair of strings in the input. Try out different values of k. Also extend and run the NW program to compute the Needleman-Wunsch similarity for each pair of sequences. Do the two different methods give the same relative order of sequence pairs, i.e, which pairs are declared to be the most similar, etc? When do you think the kmer approach might be better at finding biologically related sequences, than are alignment based approaches? Which k should one use?