fasta.bioch.virginia.edu/mol_evol

Molecular Evolution - Similarity Searching Exercises

These exercises use programs on the FASTA WWW Search page and the Molecular Evolution BLAST WWW Search page [pgm].

In the links below, [pgm] indicates a link with most of the information filled in; e.g. the program name, query, and library. [seq] links go to the NCBI, for more information about the sequence.

1. Use the FASTA search page [pgm] to compare Drosophila glutathione transferase GSTT1_DROME [seq] (gi|121694) to the PIR1 Annotated protein sequence database. Be sure to press , not .

1. Take a look at the output.
   1. How long is the query sequence?
   2. How many sequences are in the PIR1 database?
   3. What scoring matrix was used?
   4. What were the gap penalties?
   5. What are each of the numbers after the description of the library sequence? Which one is best for inferring homology?
   6. Looking at an alignment, where are the boundaries of the alignment (the best local region)?

   [show answer]

   Answer: (i) 209 aa; (ii) 13,351; (iii) BLOSUM50; (iv) gap open= -10; extend -2; one residue gap: -12; two residue: -14; (v) library sequence length; (band) optimized (opt) score (raw score), bit score, E()-value, percent identity, fraction identical, alignment length;

   Command line version.

2. What is the highest scoring non-homolog? (The non-homolog with the highest alignment score, or the lowest E()-value.) How would you confirm that your candidate non-homolog was truly unrelated?

3. Note that this drosophila glutathione transferase shares significant similarity with both sequences from bacteria (SSPA_ECO57, stringent starvation protein) and mammals. How might you test whether the stringent starvation protein is homologous to glutathione transferases? (Hint - compare your candidate non-homolog with SwissProt for a more reliable test.)

4. Compare the expectation (E()) value for the distant relationship between GSTT1_DROME and GSTP1_HUMAN (class-mu). How would you demonstrate that GSTT1_DROME is homologous to GSTP1_HUMAN?

   [show answer]

   Answer: Try searching PIR1 with an intermediate (but significant) mammalian sequence.

5. Examine how the expectation value changes with different scoring matrices (BlastP62, PAM250, MD40) and different gap penalties. (The default scoring matrix for the FASTA programs is BLOSUM50, with gap penalties of -10 to open a gap and -2 for each residue in the gap - e.g. -12 for a one residue gap).
   1. What happens to the E()-value for the 100% identical sequence with the different matrices and different gap penalties?
   2. What happens to the E()-value for distant homologs, like GSTA1_RAT with the different matrices and different gap penalties?
   3. What happens to the E()-value for the highest scoring unrelated sequence with the different matrices?

   [show answer]

   Best scores	BL50 -10/-2		BL50 -12/-2		BL62 -11/-1		PAM250 -10/-2
   GSTT1_DROME	1399	1.1e-98	1399	1.5e-109	1108	3.4e-117	1087	1.2e-75
   GSTF3_MAIZE	173	3.5e-06	163	1.8e-07	119	3e-07	152	1.1e-05
   GSTF1_MAIZE	151	7.1e-05	140	1.4e-05	96	0.0001	118	0.0035
   SSPA_ECO57	140	9.7e-05	134	4.2e-05	99	4.8e-05	140	7.8e-05
   GSTA1_RAT	139	0.00012	120	0.00064	78	0.011	118	0.037
   GSTA4_RAT	97	0.16	84	0.6	55	3.9	86	0.8
   GSTM2_RAT	93	0.32	85	0.48	61	0.83	75	5.9
   GSTP1_HUMAN	82	2.0	74	3.7	--	--	--	--
   Highest unrelated	79	3.5	76	2.6	58	3	93	1.7

6. Try the search with ssearch [pgm] (Smith-Waterman). Again, look at the E()-values for distant homologs and the highest scoring unrelated sequence.

2. Do the same search (121694) using the Course BLAST [pgm] WWW page.

1. Take a look at the output.
   1. How long is the query sequence?
   2. How many sequences are in the PIR1 database?
   3. What scoring matrix was used?
   4. What were the gap penalties?
   5. What are the numbers after the description of the library sequence? Which one is best for inferring homology?
   6. Looking at an alignment, where are the boundaries of the alignment (the best local region)?

2. What is the highest scoring non-homolog?

3. How do the blastp E()-values compare with the FASTA (blosum62) E()-values for the distantly related mammalian and plant sequences?

4. Try a limited search at the NCBI/BLAST web site. Compare Drosophila GSTT1_DROME (121694) to the Human Refseq proteins using BLASTP [pgm].
   1. What is the highest scoring unrelated sequence (hint, it has an E()-value < 0.001).
   2. How would you demonstrate that this "significant" unrelated sequence is in fact unrelated.

The default scoring matrices/gap-penalties for BLAST (BLOSUM62, -11/-1) and FASTA (BLOSUM50, 10/-2) are chosen to identify the most distant homologs possible. BLOSUM50 provides about 0.48 bits of significance per aligned amino acid, on average (ungapped). BLOSUM62 provides about 0.70 bits/aligned position (ungapped). In a search against SwissProt (~450,000 proteins, 1.7 x 10⁸ residues) with an "average" length query (400 aa), an alignment must score about 45 bits to produce E() < 0.001. For BLOSUM62, this means the alignment must be 45 bits/0.70 bits/aa = 64 residues long. While protein sequences are almost always longer than 64 residues, metagenomic DNA reads are often shorter than 200 nt (=66 aa max). For shorter DNA reads, shallower scoring matrices are required. BLAST (BLASTP, BLASTX) only offers one shallow scoring matrix, PAM30, which provides about 2.6 bits per position. FASTA/FASTX provides MD40 (2.2 bits/position, MD20 - 2.9, and MD10 - 3.4). A 70 nt DNA sequence can produce 60 bits with PAM30 or 68 bits with MD20.

1. Try the FASTA [pgm] search in question (1) using a fragment of gstt1_drome (51-100) (be sure to check the box "Use subset range" on the Query sequence" line.
   1. What is the most distantly related homolog with E() < 0.001? E()<0.01?
   2. How long is it's alignment? What is the percent identity?
   3. What would be the "right" scoring matrix for a 50 residue query to produce 40 bits of significance (40 bits, rather than 45 bits, because PIR1 only has 13,000, rather than 450,000 sequences - 1/32 as many, or 5 bits). (BLOSUM80 and PAM120 both produce about 0.91 bits of score per position.)

      Note that you can also use the "percent identity" to estimate the most effective matrix. 80% identity -> MD20, 60% identity -> MD40, 40% -> BLOSUM80.

   4. Do the same FASTA [pgm] search using the "right" scoring matrix.

2. Try the FASTA [pgm] search in question (1) using a fragment of gstt1_drome (51-80, or 90 nt). Again, be sure to check the box "Use subset range" on the Query sequence" line.
   1. Now what is the most sensitive scoring matrix?
   2. What is the percent identity of the most distant alignment?
   4. What sequences are significant if you use BLOSUM62 -11/-1 (BlastP62) as the scoring matrix? (BLAST default.)

In each of the exercises below, the BLASTP, BLASTX etc. links are pre-set to search Arabidopsis sequences.

1. BLASTP Compare the GSTT1_DROME [seq] (gi|121694) protein sequence to Arabidopsis proteins using BLASTP [pgm].

   What are the E()-values for Arabidopsis THETA1, phi, zeta-class 1and TAU

2. BLASTX

   Try the same search using the GSTT1_DROME cDNA DMGST [seq] (gi|8033) against Arabidopsis proteins using BLASTX [pgm]. If the search fails, try turning off low complexity filtering.

   What are the E()-values for Arabidopsis THETA-1, phi, zeta-class 1and TAU?

3. TBLASTN. Use GSTT1_DROME [seq] (gi|121694) against translated Arabidopsis DNA using TBLASTN [pgm].

   What are the E()-values for Arabidopsis THETA-1, phi, zeta-class 1and TAU?

4. Finally, try the DNA:DNA comparison. Use BLASTN [pgm] to compare dmgst (gi|8033) to the Reference mRNA sequences (refseq_rna) sequences in Arabidopsis.

   Are there detectable Arabidopsis homologues?

5. Use the PRSS shuffle [pgm] program to evaluate the statistical significance of a match.

1. Use PRSS [pgm] to compare GSTT1_DROME [seq] (gi|121694) to GSTA4_RAT [seq] (gi|121714)

   What is the E()-value? What equivalent database size is used to calculate the E()-value? Why isn't the database size 1? What should it be?

2. Use PRSS [pgm] to compare OPSD_HUMAN [seq] (gi|129207) to US27_HCMVA [seq] (gi|137159) Human Herpes GPCR homolog. Compare with uniform or window shuffling. How does the E()-value change? Why?

6. Exploring domains with local alignments --- Calmodulin

1. Use lalign [pgm] to examine local similarities between calmodulin CALM_HUMAN and itself.
2. Use plalign [pgm] to plot the same alignment. How many repeats are present in this sequence.

   Here are the two (lalign/plalign) alignments side by side.

3. What happens to the domain alignment plot when you use a shallower scoring matrix (try BlastP62, MD20).

7. Exploring domains with local alignments --- Cortactin (SRC8_HUMAN)

1. Use lalign [pgm] to examine local similarities between SRC8_HUMAN and itself. Note the average percent identity for the some of the most significant alignments.

2. Use plalign [pgm] to plot the same alignment. How many repeats are shown in this alignment? How long do they appear to be?

3. Based on the percent identities you saw in part (a), what would the appropriate scoring matrix be to accurately identify the cortactin domains?

4. What happens to the domain alignment plot when you use a shallower scoring matrix (try BlastP62, MD40). How do the BlastP62, PAM120, and MD40 alignments differ? Why are they different? Which matrix appears to best identify all the repeats found in the PFAM diagram?

5. You can look at the PFAM annotation of this protein at PFAM: SRC8_HUMAN [pfam] (Cortactin). How many repeats are shown in this diagram? How long are they?

For more complex domain alignments, try mwkw, or mouse RNA polymerase (rpb1_mouse residues 1500- ) against itself. Try the rpb1_mouse alignment using the MD20 scoring matrix as well as BLOSUM50.

Where to get the FASTA package: faculty.virginia.edu/wrpearson/fasta/

The "normal" FASTA WWW site:

Contact Bill Pearson: wrp@virginia.edu