1 / 70

Comparing biological sequences (3): Database searching and Multiple alignment

Comparing biological sequences (3): Database searching and Multiple alignment. Database searching. Goal: find similar (homologous) sequences of a query sequence in a sequence of database Input: query sequence & database Output: hits (pairwise alignments). Database searching.

carol-york
Download Presentation

Comparing biological sequences (3): Database searching and Multiple alignment

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Comparing biological sequences (3):Database searching and Multiple alignment

  2. 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. Database searching • Core: pair-wise 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

  4. What program to use for searching? 1) BLAST is fastest and easily accessed on the Web • limited sets of databases • nice translation tools (BLASTX, TBLASTN) 2) FASTA • precise choice of databases • more sensitive for DNA-DNA comparisons • FASTXand TFASTX can find similarities in sequences with frameshifts 3) Smith-Waterman is slower, but more sensitive • known as a “rigorous” or “exhaustive” search • SSEARCHin GCG and standalone FASTA

  5. FASTA 1) Derived from logic of the dot plot • compute best diagonals from all frames of alignment 2) 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

  6. FASTA Algorithm

  7. Makes Longest Diagonal 3) after all diagonals found, tries to join diagonals by adding gaps 4) computes alignments in regions of best diagonals

  8. FASTA Alignments

  9. FASTA Results - Histogram !!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:=============================================*

  10. FASTA Results - List 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

  11. FASTA Results - Alignment 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

  12. FASTAon 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.

  13. Institut de Génétique Humaine, Montpellier France, GeneStream server • http://www2.igh.cnrs.fr/bin/fasta-guess.cgi • Oak Ridge National Laboratory GenQuest serverhttp://avalon.epm.ornl.gov/ • European Bioinformatics Institute, Cambridge, UKhttp://www.ebi.ac.uk/htbin/fasta.py?request • EMBL, Heidelberg, Germanyhttp://www.embl-heidelberg.de/cgi/fasta-wrapper-free • Munich Information Center for Protein Sequences (MIPS)at Max-Planck-Institut, Germanyhttp://speedy.mips.biochem.mpg.de/mips/programs/fasta.html • Institute of Biology and Chemistry of Proteins Lyon, Francehttp://www.ibcp.fr/serv_main.html • Institute Pasteur, Francehttp://central.pasteur.fr/seqanal/interfaces/fasta.html • GenQuest at The Johns Hopkins University • http://www.bis.med.jhmi.edu/Dan/gq/gq.form.html • National Cancer Center of Japan • http://bioinfo.ncc.go.jp

  14. BLAST Searches GenBank [BLAST= Basic Local Alignment Search Tool] The NCBI BLASTweb server lets you compare your query sequence to various sections of GenBank: • nr = non-redundant (main sections) • month = new sequences from the past few weeks • ESTs • human, drososphila, yeast, or E.coli genomes • proteins (by automatic translation) • This is a VERY fast and powerful computer.

  15. BLAST • 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 • Does not handle gaps well • New “gapped BLAST” (BLAST 2) is better

  16. BLAST Algorithm

  17. BLAST Word Matching MEAAVKEEISVEDEAVDKNI MEA EAA AAV AVK VKE KEE EEI EIS ISV ... Break query into words: Break database sequences into words:

  18. Database Sequence Word Lists RTT AAQ SDG KSS SRW LLN QEL RWY VKI GKG DKI NIS LFC WDV AAV KVR PFR DEI … … Query Word List: MEA EAA AAV AVK VKL KEE EEI EIS ISV Compare Word Lists ? Compare word lists by Hashing (allow near matches)

  19. Find locations of matching words in database sequences ELEPRRPRYRVPDVLVADPPIARLSVSGRDENSVELTMEAT • MEA • EAA • AAV • AVK • KLV • KEE • EEI • EIS • ISV TDVRWMSETGIIDVFLLLGPSISDVFRQYASLTGTQALPPLFSLGYHQSRWNY IWLDIEEIHADGKRYFTWDPSRFPQPRTMLERLASKRRVKLVAIVDPH

  20. Extend hits one base at a time

  21. HVTGRSAF_FSYYGYGCYCGLGTGKGLPVDATDRCCWA Seq_XYZ: QSVFDYIYYGCYCGWGLG_GK__PRDA Query: E-val=10-13 • Use two word matches as anchors to build an alignment between the query and a database sequence. • Then score the alignment.

  22. HSPs are Aligned Regions • The results of the word matching and attempts to extend the alignment are segments - called HSPs (High-scoring Segment Pairs) • BLAST often produces several short HSPs rather than a single aligned region

  23. BLAST 2algorithm • The NCBI’s BLAST website now both use BLAST 2 (also known as “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

  24. Statistical tests • Evaluate the probability of an event taking place by chance (at random). • P-value • Randomized data • Distribution under the same setup • Z-score • Chebyshev Inequality

  25. BLASTStatistics • E value is equivalent to standard P value (based on Karlin-Altschul theorem) • Significant if E < 0.05 (smaller numbers are more significant) • The E-value represents the likelihood that the observed alignment is due to chance alone. A value of 1 indicates that an alignment this good would happen by chance with any random sequence searched against this database.

  26. BLAST variants for different searchesa (after S. Brenner, Trends Guide to Bioinformatics, 1998)

  27. 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

  28. Interpretation of output • very low E values (e-100) are homologs or identical genes • 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

  29. 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?

  30. Borderline similarity • What to do with matches with E() values in the 0.5 -1.0 range? • this is the “Twilight Zone” • retest these sequences and look for related hits (not just your original query sequence) • similarity is transitive: if A~B and B~C, then A~C

  31. Position Specific Iterated BLAST • Collect all database sequence segments that have been aligned with query sequence with E-value below set threshold (default 0.01) • Construct position specific scoring matrix for collected sequences. Rough idea: • Align all sequences to the query sequence as the template. • Assign weights to the sequences • Construct position specific scoring matrix • Iterate

  32. Motif finding • Observation : Some regions have been better conserved than others during evolution • Idea: By analyzing the constant and variable properties of such groups of similar sequences, it is possible to derive a signature for a protein family or domain (motifs)

  33. PROSITE patterns • PROSITE fingerprints are described by regular grammars • There is a number of programs that allow to search databases for PROSITE patterns (example GCG package) • Example [EDQH]-x-K-x-[DN]-G-x-R-[GACV] • Rules: • Each position is separated by a hyphen • One character denotes residuum at a given position • […] denoted a set of allowed residues • (n) denotes repeat of n • (n,m) denoted repeat between n and m inclusive Ex. ATP/GTP binding motive [SG]=X(4)-G-K-[DT]

  34. Multiple sequence alignment

  35. Generalizing the Notion of Pairwise Alignment • Alignment of 2 sequences is represented as a 2-row matrix • In a similar way, we represent alignment of 3 sequences as a 3-row matrix A T _ G C G _ A _ C G T _ A A T C A C _ A • Score: more conserved columns, better alignment

  36. Alignments = Paths • Align 3 sequences: ATGC, AATC,ATGC

  37. Alignment Paths x coordinate

  38. Alignment Paths • Align the 3 sequences: ATGC, AATC,ATGC x coordinate y coordinate

  39. Alignment Paths x coordinate y coordinate z coordinate • Resulting path in (x,y,z) space: • (0,0,0)(1,1,0)(1,2,1) (2,3,2) (3,3,3) (4,4,4)

  40. Aligning Three Sequences source • Same strategy as aligning two sequences • Use a 3-D “Manhattan Cube”, with each axis representing a sequence to align • For global alignments, go from source to sink sink

  41. 2-D vs 3-D Alignment Grid V W 2-D edit graph 3-D edit graph

  42. 2-D cell versus 2-D Alignment Cell In 2-D, 3 edges in each unit square In 3-D, 7 edges in each unit cube

  43. Architecture of 3-D Alignment Cell (i-1,j,k-1) (i-1,j-1,k-1) (i-1,j,k) (i-1,j-1,k) (i,j,k-1) (i,j-1,k-1) (i,j,k) (i,j-1,k)

  44. si-1,j-1,k-1 + (vi, wj, uk) si-1,j-1,k + (vi, wj, _ ) si-1,j,k-1 + (vi, _, uk) si,j-1,k-1 + (_, wj, uk) si-1,j,k + (vi, _ , _) si,j-1,k + (_, wj, _) si,j,k-1 + (_, _, uk) Multiple Alignment: Dynamic Programming cube diagonal: no indels • si,j,k = max • (x, y, z) is an entry in the 3-D scoring matrix face diagonal: one indel edge diagonal: two indels

  45. Multiple Alignment: Running Time • For 3 sequences of length n, the run time is 7n3; O(n3) • For ksequences, build a k-dimensional Manhattan, with run time (2k-1)(nk); O(2knk) • Conclusion: dynamic programming approach for alignment between two sequences is easily extended to k sequences but it is impractical due to exponential running time.

  46. Profile Representation of Multiple Alignment - A G G C T A T C A C C T G T A G – C T A C C A - - - G C A G – C T A C C A - - - G C A G – C T A T C A C – G G C A G – C T A T C G C – G G A 1 1 .8 C .6 1 .4 1 .6 .2 G 1 .2 .2 .4 1 T .2 1 .6 .2 - .2 .8 .4 .8 .4

  47. Profile Representation of Multiple Alignment - A G G C T A T C A C C T G T A G – C T A C C A - - - G C A G – C T A C C A - - - G C A G – C T A T C A C – G G C A G – C T A T C G C – G G A 1 1 .8 C .6 1 .4 1 .6 .2 G 1 .2 .2 .4 1 T .2 1 .6 .2 - .2 .8 .4 .8 .4 In the past we were aligning a sequence against a sequence Can we align a sequence against a profile? Can we align a profile against a profile?

  48. Aligning alignments • Given two alignments, can we align them? x GGGCACTGCAT y GGTTACGTC-- Alignment 1 z GGGAACTGCAG w GGACGTACC-- Alignment 2 v GGACCT-----

  49. Aligning alignments • Given two alignments, can we align them? • Hint: use alignment of corresponding profiles x GGGCACTGCAT y GGTTACGTC-- Combined Alignment z GGGAACTGCAG w GGACGTACC-- v GGACCT-----

  50. Multiple Alignment: Greedy Approach • Choose most similar pair of strings and combine into a profile , thereby reducing alignment of k sequences to an alignment of of k-1 sequences/profiles. Repeat • This is a heuristic greedy method u1= ACg/tTACg/tTACg/cT… u2 = TTAATTAATTAA… … uk = CCGGCCGGCCGG… u1= ACGTACGTACGT… u2 = TTAATTAATTAA… u3 = ACTACTACTACT… … uk = CCGGCCGGCCGG k-1 k

More Related