Sequence Database Searching (Basic Tools and Advanced Methods)

1. I519 Introduction to Bioinformatics, Fall 2012 Sequence Database Searching (Basic Tools and Advanced Methods)

2. What’s database searching • Goal: find similar (homologous) sequences of a query sequence in a sequence of database • Input: query sequence & database • Output: hits (pairwise alignments)

3. Basics of database searching • Core: pairwise alignment algorithm • Speed (fast sequence comparison) • Relevance of the search results (statistical tests) • Recovering all information of interest • The results depend of the search parameters like gap penalty, scoring matrix. • Sometimes searches with more than one matrix should be preformed • Specificity and sensitivity

4. The different strategies for database searching Faster Higher sensitivity

5. What program to use for searching? • BLAST is fastest and easily accessed on the Web • A suite; BLASTP, BLASTN, BLASTX • lag behind the development of sequencing techniques • FASTA • more sensitive for DNA-DNA comparisons • FASTX and TFASTX can find similarities in sequences with frameshifts • Smith-Waterman is slower, but more sensitive • known as a “rigorous” or “exhaustive” search • SSEARCH in GCG and standalone FASTA

6. FASTA • Derived from logic of the dot plot • compute best diagonals from all frames of alignment • Word method looks for exact matches between words in query and test sequence • hash tables (fast computer technique) • DNA words are usually 6 bases • protein words are 1 or 2 amino acids • only searches for diagonals in region of word matches = faster searching

7. FASTA algorithm

8. Makes longest diagonal • after all diagonals found, tries to join diagonals by adding gaps • computes alignments in regions of best diagonals

9. FASTA alignments

10. !!SEQUENCE_LIST 1.0 (Nucleotide) FASTA of: b2.seq from: 1 to: 693 December 9, 2002 14:02 TO: /u/browns02/Victor/Search-set/*.seq Sequences: 2,050 Symbols: 913,285 Word Size: 6 Searching with both strands of the query. Scoring matrix: GenRunData:fastadna.cmp Constant pamfactor used Gap creation penalty: 16 Gap extension penalty: 4 Histogram Key: Each histogram symbol represents 4 search set sequences Each inset symbol represents 1 search set sequences z-scores computed from opt scores z-score obs exp (=) (*) < 20 0 0: 22 0 0: 24 3 0:= 26 2 0:= 28 5 0:== 30 11 3:*== 32 19 11:==*== 34 38 30:=======*== 36 58 61:===============* 38 79 100:==================== * 40 134 140:==================================* 42 167 171:==========================================* 44 205 189:===============================================*==== 46 209 192:===============================================*===== 48 177 184:=============================================* FASTA results: histogram

11. The best scores are: init1 initn opt z-sc E(1018780).. SW:PPI1_HUMAN Begin: 1 End: 269 ! Q00169 homo sapiens (human). phosph... 1854 1854 1854 2249.3 1.8e-117 SW:PPI1_RABIT Begin: 1 End: 269 ! P48738 oryctolagus cuniculus (rabbi... 1840 1840 1840 2232.4 1.6e-116 SW:PPI1_RAT Begin: 1 End: 270 ! P16446 rattus norvegicus (rat). pho... 1543 1543 1837 2228.7 2.5e-116 SW:PPI1_MOUSE Begin: 1 End: 270 ! P53810 mus musculus (mouse). phosph... 1542 1542 1836 2227.5 2.9e-116 SW:PPI2_HUMAN Begin: 1 End: 270 ! P48739 homo sapiens (human). phosph... 1533 1533 1533 1861.0 7.7e-96 SPTREMBL_NEW:BAC25830 Begin: 1 End: 270 ! Bac25830 mus musculus (mouse). 10, ... 1488 1488 1522 1847.6 4.2e-95 SP_TREMBL:Q8N5W1 Begin: 1 End: 268 ! Q8n5w1 homo sapiens (human). simila... 1477 1477 1522 1847.6 4.3e-95 SW:PPI2_RAT Begin: 1 End: 269 ! P53812 rattus norvegicus (rat). pho... 1482 1482 1516 1840.4 1.1e-94 FASTA results: list

12. SCORES Init1: 1515 Initn: 1565 Opt: 1687 z-score: 1158.1 E(): 2.3e-58 >>GB_IN3:DMU09374 (2038 nt) initn: 1565 init1: 1515 opt: 1687 Z-score: 1158.1 expect(): 2.3e-58 66.2% identity in 875 nt overlap (83-957:151-1022) 60 70 80 90 100 110 u39412.gb_pr CCCTTTGTGGCCGCCATGGACAATTCCGGGAAGGAAGCGGAGGCGATGGCGCTGTTGGCC || ||| | ||||| | ||| ||||| DMU09374 AGGCGGACATAAATCCTCGACATGGGTGACAACGAACAGAAGGCGCTCCAACTGATGGCC 130 140 150 160 170 180 120 130 140 150 160 170 u39412.gb_pr GAGGCGGAGCGCAAAGTGAAGAACTCGCAGTCCTTCTTCTCTGGCCTCTTTGGAGGCTCA ||||||||| || ||| | | || ||| | || || ||||| || DMU09374 GAGGCGGAGAAGAAGTTGACCCAGCAGAAGGGCTTTCTGGGATCGCTGTTCGGAGGGTCC 190 200 210 220 230 240 180 190 200 210 220 230 u39412.gb_pr TCCAAAATAGAGGAAGCATGCGAAATCTACGCCAGAGCAGCAAACATGTTCAAAATGGCC ||| | ||||| || ||| |||| | || | |||||||| || ||| || DMU09374 AACAAGGTGGAGGACGCCATCGAGTGCTACCAGCGGGCGGGCAACATGTTTAAGATGTCC 250 260 270 280 290 300 240 250 260 270 280 290 u39412.gb_pr AAAAACTGGAGTGCTGCTGGAAACGCGTTCTGCCAGGCTGCACAGCTGCACCTGCAGCTC |||||||||| ||||| | |||||| |||| ||| || ||| || | DMU09374 AAAAACTGGACAAAGGCTGGGGAGTGCTTCTGCGAGGCGGCAACTCTACACGCGCGGGCT 310 320 330 340 350 360 FASTA results: alignment

13. FASTA on the web • Many websites offer FASTA searches • Various databases and various other services • Be sure to use FASTA 3 • Each server has its limits • Be aware that you are depending on the kindness of strangers.

14. BLAST • [BLAST= Basic Local Alignment Search Tool] • The central idea of the BLAST algorithm is that a statistically significant alignment is likely to contain a high-scoring pair of aligned words. • Uses word matching like FASTA • Similarity matching of words (3 aa’s, 11 bases) • Does not require identical words. • If no words are similar, then no alignment • won’t find matches for very short sequences • Original BLAST does not handle gaps well; the “gapped” blast is better

15. Word match Extend hits

16. MEAAVKEEISVEDEAVDKNI MEA EAA AAV AVK VKE KEE EEI EIS ISV ... BLAST: word matching Break query into words: Break database sequences into words: Compare word lists by Hashing (allow near matches)

17. BLAST: extend hits Use two word matches as anchors to build an alignment between the query and a database sequence.

18. The NCBI’s BLAST website now both use “gapped BLAST” This algorithm is more complex than the original BLAST It requires two word matches close to each other on a pair of sequences (i.e. with a gap) before it creates an alignment allow gaps (using Smith-Waterman algorithm constrained by a band width). Gapped BLAST algorithm

19. The region of the path graph explored when seeded by the alignment of a pair of residues Ref: Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. NAR 25(17):3389, 1997

20. Evaluate the probability of an event taking place by chance. Very little is known about the random distribution of optimal global alignment scores Statistics for the scores of local alignments, are well understood (particularly true for local alignments lacking gaps) P-value Randomized data Distribution under the same setup Statistical tests

21. E-value is equivalent to standard P value (based on Karlin-Altschul theorem) Karlin-Alschul equation is probably one of the most recognized equation in bioinformatics. E-value represents the likelihood that the observed alignment is due to chance alone E-value: the number of alignments expected by chance (E) during a sequence database search. E = 1 indicates that an alignment this good would happen by chance with any random sequence searched against the same database. BLAST statistics

22. How to compute BLAST E-value Raw score (S) Bit score (S’) m: the effective query length n: the effective database length E P (the probability of getting at least one HSP)

23. Which BLAST? BLAST+ package has different user interface

24. Which database? • For blastp (blastx) • Non-redundant protein sequences (nr) • Reference proteins (refseq_protein) • Swissprot protein sequences (swissprot) • Patented protein sequences (pat) • Protein Data Bank proteins (pdb) • Environmental samples (env_nr) • For blastn (tblastx) • Human genomic plus transcript (not for tblastx) • Mouse genomic plus transcript (not available for tblastx) • Nucleotide collection (nr/nt) • etc

25. BLAST is approximate • BLAST makes similarity searches very quickly because it takes shortcuts. • looks for short, nearly identical “words” (11 bases) • It also makes errors • misses some important similarities • makes many incorrect matches • easily fooled by repeats or skewed composition

26. Interpretation of BLAST hits • very low E values (e-100) are homologs • moderate E values are related genes • long list of gradually declining of E values indicates a large gene family • long regions of moderate similarity are more significant than short regions of high identity

27. Is the hit with smallest E-value the closest sequence to the query? • Not necessarily • Some people argue that more strict phylogeny analysis is needed for further conclusion • Why? Different evolutionary rates, etc

28. Biological relevance • It is up to you, the biologist to scrutinize these alignments and determine if they are significant. • Were you looking for a short region of nearly identical sequence or a larger region of general similarity? • Are the mismatches conservative ones? • Are the matching regions important structural components of the genes or just introns and flanking regions?

29. 20 tips to improve your BLAST searches • “Don't Use the Default Parameters” (but what parameters?) • “Treat BLAST Searches as Scientific Experiments • “View BLAST Reports Graphically” (depends!) • “Be Skeptical of Hypothetical Proteins” • “Look for Stop Codons and Frame-Shifts to find Pseudo-Genes” • “Parse BLAST Reports with Bioperl” (your own python, or perl, program) • “How to Lie with BLAST Statistics” (“Lies, dammed lies, and statistics”) • … Ref: BLAST by Ian Korf; Mark Yandell; Joseph Bedell, 2003.

30. Borderline similarity • What to do with matches with E values that are not so impressive? (e.g., E values > 0.01, smaller numbers are more significant) • this is the “Twilight Zone” • retest these sequences and look for related hits (not just your original query sequence) • similarity is (is not) transitive: • if A~B and B~C, then A~C

31. Distant homology detection • Use more data (MSA) • Use advanced approaches (like HMM)

32. PSI-BLAST • Position Specific Iterated BLAST (PSI-BLAST) • Basic idea • Use results from BLAST query to construct a profile matrix (or PSSM) • Search database with profile instead of query sequence • Iterate • Position-specific scoring matrices are an extension of substitution scoring matrices

33. PSSM--Position Specific Scoring Matrix 3.0 Match A to position 3 Match D to position 3 -1.0 Match D to position 2 distribution of amino acids + pseudo-counts

34. Representing a profile as a logo • Logos are used to show the residue preferences or conservation at particular positions • Based on information theory helix-turn-helix motif from the CAP family of homodimeric DNA binding proteins http://weblogo.berkeley.edu/examples.html

35. Profile-profile alignment (alignment of alignments) • Inputs: two MSAs (profiles) • Alignment of two MSAs (similar to pairwise sequence alignment) • Many existing programs • Scoring functions vary • Dot product (FFAS)

36. FFAS scoring function Dot product of two amino acid “frequency” vectors

37. Sequence identity scoring zones • >25-30%: homology zone • 15-25%: twilight zone • <15%: midnight zone Weak sequence similarity detection is still not solved! And sequence similarity != structural similarity (for proteins) != functional similarity

38. Readings • Chapter 5 (Pairwise Sequence Alignment and Database Searching) • Box 5.1 – saving space by throwing away the intermediate calculations • “Time can be saved with a loss of rigor by not calculating the whole matrix” (Figure 5.18)