Øvelse: BLAST

Øvelse skrevet af: Rasmus Wernersson

Introduktion

I denne øvelse skal vi bruge BLAST (Basic Local Alignment Search Tool) til at søge i sekvensdatabaser såsom GenBank (DNA) og UniProt (protein). Et meget væsentligt aspekt mht. brug af BLAST er at kunne vurderer signifikansen af et givet resultat.

BLAST er en program-pakke der kan hentes kvit og frit, og som indeholder programmer til både at oprette nye søge-database (fortmatdb programmet) og til at søge i disse database (fx blastn til DNA-database og blastp til proteindatabaser). Vi skal her bruge det webinterface til BLAST som findes hos NCBI: NCBI BLAST - en af de store fordele er (ud over at vi ikke skal rode med en UNIX prompt) at alle NCBI's databaser (+ UniProt/SwissProt !) er tilgængelige til at søge i.

Links:

• NCBI BLAST hovedside: http://www.ncbi.nlm.nih.gov/BLAST/
• Bemærk: Der er links til "Nucleotide BLAST" (incl "blastn") og "Protein BLAST" (incl "blastp") herfra.
  
• NCBI "BLAST information GUIDE"
  

Del 1 - alignment af tilfældige sekvenser.

Som gennemgået i forelæsningen, vil er være risiko for at få falske positive resultater (altså et hit til en ikke-relateret sekvens) af rent stokastiske årsager. I første del af øvelsen skal vi kigge på hvad der sker hvis man søge med rent tilfældige sekvenser som input.

I har alle fået udleveret en fire-sidet og en tyve-sidet terning, som vi skal bruge til at generere DNA- og protein-sekvenser (hvis du laver denne øvelser derhjemme så kan du istedet bruge dette lille Java program "Remote Roller") .

TRIN 1 - DNA sekvenser og BLASTN:

1. Start med at slå 3 stk 25bp DNA sekvenser, ud fra følgende tabel.
   
   
   

1 : A
2 : C
3 : G
4 : T

2. Vi skal nu have fat i BLASTN algoritmen hos NCBI.
   
1. Vælg "Nucleotide BLAST" fra BLAST hovedsiden.
   
   
2. Under "Program Selection" skal I vælge "Somewhat similar sequences (blastn)"
   
   
3. Vælg "Nucletide Collection NR/NT" som database. NR står for "Non Redundancy" og er en samling af alle ikke-redundante sekvenser fra GenBank og genom-sekventerings projekterne.
   
   
4. VIGTIGT: Til denne kørsel med de små kunstige DNA sekvenser, skal vi slå noget af NCBI automatik mht. valg af options fra. Tryk på "Algorithm parameters"for at få adgang til de bagvedliggende options (se nedenstående screenshot). Fravælg "Automaticlly adjust parameters for short input sequences". Bemærk: her kan man også indstille E-value cutoff ("Expect threshold") og BLASTN "word size".
   
   
5. Indsæt dine tre sekvenser (hint: FASTA format) - og start BLAST søgningen.
   
   
   
   Svar: Husk at det er nemt at skrive sekvenser op i FASTA format i hånden - hver sekvens skal blot have en header der fortæller om dens navn:
   
   >seq_01
   CGCCCGACCGTGTAGGAGGTCCGGT
   >seq_02
   TTCCAAGTAGGTTATGTTAAAGCCG
   >seq_03
   CGGTTAATTCGCTGTCTATATAAGG
   
   
3. Lad os kigge lidt på resultaterne - bemærk at der både er en tabel med en oversigt i toppen, samt at selve de individuelle alignments vises længere nede på siden:
   
   
• Hvor stor (i basepar) er den database vi har BLAST'et imod?
  Svar: 23 milliarder bogstaver (dvs. baser i tilfældet med DNA).
  
  
• Er der nogen sekvenser de ligner vores tilfældigt genererede?
  Svar: Der er flere perfekte hits, fx:
  
  
    GENE ID: 5438907 BC1G_02961 | hypothetical protein
  [Botryotinia fuckeliana B05.10]
  
   Score = 33.7 bits (36),  Expect = 4.8
   Identities = 18/18 (100%), Gaps = 0/18 (0%)
   Strand=Plus/Minus
  
  Query  1    CGCCCGACCGTGTAGGAG  18
              ||||||||||||||||||
  Sbjct  574  CGCCCGACCGTGTAGGAG  557
  
  
• Hvor lange er de hits der findes? Svar: Typisk en 17-19 baser.
• Hvor mange % identity? Svar: 95% - 100%
• I hvilket område ligger bit-scores ("max score")? Svar: typisk på 33 - 35 bits.
• I hvilket område ligger E-value?
  Svar: 1.4 - 4.8 (for denne søgning, ligger typisk fra 0.1 - 10.0). BEMÆRK: at default cut-off for e-value er 10.0 - vi ser derfor ikke hits der har en dårligere (større) e-value end 10.0.
  
  
• Giver disse resultater nogen biologisk mening?
  Svar: NEJ! - hits'ne er rigtige nok, det er faktisk sekvens der findes i databasen. Men vi ved jo at vores sekvenser er helt tilfældige, og derfor ikke har noget slægtskab med dem i databasen. Vi har kun fundet vores hits, fordi databasen er så enorm stor, at vi af rent stokastiske årsager rammer nogle sekvenser der ligner.
  
  E-værdierne fortæller netop dette: Som jeg gennemgik i forelæsningenen, vil alignment score ("bit score" i BLAST) følge en ekstremværdi-fordeling for de sekvenser der ikke har noget med vores input-sekvenser at gøre (fordelt omkring en lav bit-score), og de hits der faktisk er relateret følger en normalfordeling (omkring en væsentlig højere bit-score).
  
  Ud fra ekstremværdi-fordelingen (som beregnes når BLAST databasen oprettes), er det mulig at estimere hvor mange tilfældige, ikke-relaterede hits, man må forvente ved en given bit-score. Dette udtrykkes som en e-værdi (EXPECT-value) som fortæller, hvor mange ikke-relaterede hits man kan forvente at finde ved en given bit-score. Hvis vi fx kigger på vores ovenstående hit med en bit-score på 33.7, giver dette en e-value på 4.8 (i denne database, af 23*10⁹ baser). Man må altså kunne forvente statistisk set af finde 4.8 hits med en score på 33.7 bits i denne database - af tilfældige årsager!
  
  
4. Lad os prøve at vælgde en anden database. Vælg "Human genomic plus transcripts" som database og kør søgningen en gang til.
   
   
• Hvor stor er denne database? Svar: 6 milliader baser.
  
  
• Hvor lange er de hits der findes? Svar: 16-21 baser.
• Hvor mange % identity? Svar: 95% - 100%
• I hvilket område ligger bit-scores ("max score")? Svar: 31.9-35.6 (for mine sekvenser!)
  
• I hvilket område ligger E-value? Svar: 0.36 - 4.3
  
  
• De humane sekvenser findes også i NR - blev de samme sekvenser fundet i søgningen mod NR databasen? Hvorfor / hvorfor ikke?
  VIGTIGT (marts 2008) - NCBI har åbenbart valgt at pille de komplette humane kromosomer ud af NR/NT databasen efter jeg har skrevet øvelsen (de enkelte gener er der stadig) - det er derfor ikke lige til at finde det samme entry i begge søgninger. MEN - det er muligt at finde entries der har samme alignment score (bit-score) i begge søgninger, og så se hvordan E-value varierer fra den humane database (6 mia. baser) til NR/NT databasen (23 mia. baser). hver gang database-størrelsen fordobles, fordobles E-værdien ligeledes!
  
  
  Svar (2007): For sekvens nummer to (i mit eksempel) dukker der nu en lang række hits op - med en e-værdi på 8.4. Disse hits ser vi ikke i søgningen mod NR databasen. Forklaringen på dette er forskellen i størrelsen (NR er 3.5 gange støre end human+transcripts databasen). Et hit der har en e-værdi på 8.4 i den mindre database vil være 3.5 gange så sandsynlig at finde af tilfældige årsager i en database der er 3.5 gange større. Dermed vil det få en e-værdig på omkring 30 i den store database - disse hits ser vi ikke pga. en default cut-off e-value på 10.0.
  
  
  
5. En kort kommentar om word-size og DNA sekvenser: Som default er word-size 11 for DNA sekvenser - dette betyder at algorithmen først kigger efter et perfect match på 11 baser, og hvis det findes, så kigges der på hvad der ligger af sekvens omkring dette stykke - man kan ændre den til 7 (langsommere, men mere grundigt) eller 15 (hurtigere - finder ikke helt så meget). Det gør ikke den store forskel i dette eksempel med tilfældige sekvenser, da de samle hits altid er korte og næsten altid uden gaps - men det kan betyde en forskel, hvis man aligner rigtigt DNA, hvor resultatet er et længe alignment med både mismatches og gaps.
   Kommentar:  Som gennemgået i forelæsningen - og i PodCast'en.
   

TRIN 2 - Protein sekvenser og BLASTP:

6. Lad os nu istedet kaste os over nogle protein-sekvenser. Start med at slå 3 peptider med en længe på 25 aa, ud fra nedenstående tabel.
   
   
   

7. Åbn "Protein BLAST" siden hos NCBI og vælg BLASTP som metode. Sæt sekvenserne ind og vælg igen "NR" databasen (denne gang består den af oversætte CDS'er + SwissProt mm).
   
   Under "Algorithm Parameter"vælges BLOSUM62 som matrice, og "Expect treshold" sættes til 1000 (default: 10) VIGTIGT. Husk også at slå "Short sequence" autoparameterne fra (ligesom med DNA søgningen) - ellers ignoreres vores omhyggeligt indstillede parametere.
   
   Submit derefter søgningen.
   
   
   HUSK E-VALUE CUT-OFF PÅ 1000!
   
   Mine sekvenser:
   
   >seq_01
   aymgsplsflshdhncirdakyvvs
   >seq_02
   danshrsqwnrlgmanwfqhwsewv
   >seq_03
   vsltdgpdkhesqhirqdkyfddia
   
   
   
• Hvor stor er databasen? Svar: 2.1 milliader bogstaver (aminosyrere).
  
  
• Hvor lange er hits'ne typisk og har de gaps? Svar: Typisk 17-21. Ingen gaps, men med en del mismatches.
• Hvor mange % identity? Svar: ca. 40-60
• I hvilket område ligger E-value? Svar: 229 - 986 (!).
• Prøve at kigge grundigt på et par af selve alignments'ne - ("+" betyder similar) - er er nogle af dem der som sådan ser rimeligt plausible ud (hvis vi for et øjeblik ser bort fra længden)?
  Svar: Ja se fx. følgende alignment - det ser da overbevisnede ud (men det er alt, alt for kort til at være signifikant):
  
  
   Score = 27.7 bits (60),  Expect =   174, Method: Composition-based stats.
   Identities = 10/17 (58%), Positives = 12/17 (70%), Gaps = 0/17 (0%)
  
  Query  1   AYMGSPLSFLSHDHNCI  17
             AY GSP  +  HDHNC+
  Sbjct  49  AYPGSPHGYDIHDHNCL  65
  
  
  
• Hvis vi havde brugt default cut-off'et for E-value på 10, ville der så være fundet nogen sekvenser?
  Svar: NEJ. Hvis vi også havde kørt søgningerne mod DNA databasen med en e-value cut-off på 1000, ville vi har set sidevis af hits for hver sekvens.
  
  
• Fanger vi mest junk ved at bruge DNA eller peptid i vores søgninger?
  Svar: Risikoen for at samle et "falsk" hit op, er større med DNA - bemærk at vi var nødt til at "snyde" med cut-off for i det hele taget at se nogle hits mod proteinsekvenserne.
  
  
8. Kommentar ang. word-size og proteinsekvenser. blastp fungerer noget anderledes end blastn mht. at finde et sted et sekvenserne med høj-similariet: istedet for at kræve et længede stykke der er 100% det samme (det vil ikke give mening for protein-sekvens), prøver algoritmen at finde to sekvens-stumper (af længde "word-size", default =3) med høj similaritet inden for en kort afstand (default: 40 aa). Så snart et sådan sekvens par er fundet, fylder BLAST ind med almindeligt lokalt alignment.
   Kommentar:  Som gennemgået i forelæsningen - og i PodCast'en.
   

Del 2 - overførelse af viden om funktion

En af de allermest brugte anvendelser af BLAST er til, at overføre viden om funktion fra en sekvens til en anden. En vigtig forudsætning for dette er erkendelsen af alle levende organismer er evolutionært beslægtede. Hvis gen "X" i en organisme fx. vides at kode for et bestemt enzym (e.g. en protease) kan man med rimelig sikkerhed antage at gen "Y" fra en anden organisme har sammen funktion, hvis de to gener er tilstrækkeligt ens (hvilket betyder at de ikke er gået alt for meget deres egne veje, siden sidste fælles stamform for genet).

Vi undersøgte præcis dette fænomen i "Parvis alignment" øvelsen (med en serin-protease fra to Bacills arter). Denne gang har vi ikke to sekvenser vi iforvejen har en formodning om er beslægtede - vi skal derimod bruge BLAST til at finde de sekvenser der ligner i de samlede database med millioner af sekvenser.

Lad os starte med en situation der giver nogle ganske fornuftige hits i databasen. Nedenstående sekvens er et komplet transcript (mRNA) fra en bakterie. Lad os finde ud af hvad det er.

>Unknown_transcript01
CCACTTGAAACCGTTTTAATCAAAAACGAAGTTGAGAAGATTCAGTCAACTTAACGTTAATATTTGTTTC
CCAATAGGCAAATCTTTCTAACTTTGATACGTTTAAACTACCAGCTTGGACAAGTTGGTATAAAAATGAG
GAGGGAACCGAATGAAGAAACCGTTGGGGAAAATTGTCGCAAGCACCGCACTACTCATTTCTGTTGCTTT
TAGTTCATCGATCGCATCGGCTGCTGAAGAAGCAAAAGAAAAATATTTAATTGGCTTTAATGAGCAGGAA
GCTGTTAGTGAGTTTGTAGAACAAGTAGAGGCAAATGACGAGGTCGCCATTCTCTCTGAGGAAGAGGAAG
TCGAAATTGAATTGCTTCATGAATTTGAAACGATTCCTGTTTTATCCGTTGAGTTAAGCCCAGAAGATGT
GGACGCGCTTGAACTCGATCCAGCGATTTCTTATATTGAAGAGGATGCAGAAGTAACGACAATGGCGCAA
TCAGTGCCATGGGGAATTAGCCGTGTGCAAGCCCCAGCTGCCCATAACCGTGGATTGACAGGTTCTGGTG
TAAAAGTTGCTGTCCTCGATACAGGTATTTCCACTCATCCAGACTTAAATATTCGTGGTGGCGCTAGCTT
TGTACCAGGGGAACCATCCACTCAAGATGGGAATGGGCATGGCACGCATGTGGCCGGGACGATTGCTGCT
TTAAACAATTCGATTGGCGTTCTTGGCGTAGCGCCGAGCGCGGAACTATACGCTGTTAAAGTATTAGGGG
CGAGCGGTTCAGGTTCGGTCAGCTCGATTGCCCAAGGATTGGAATGGGCAGGGAACAATGGCATGCACGT
TGCTAATTTGAGTTTAGGAAGCCCTTCGCCAAGTGCCACACTTGAGCAAGCTGTTAATAGCGCGACTTCT
AGAGGGGTTCTTGTTGTAGCGGCATCTGGGAATTCAGGTGCAGGCTCAATCAGCTATCCGGCCCGTTATG
CGAACGCAATGGCAGTCGGAGCGACTGACCAAAACAACAACCGCGCCAGCTTTTCACAGTATGGCGCAGG
GCTTGACATTGTCGCACCAGGTGTAAACGTGCAGAGCACATACCCAGGTTCAACGTATGCCAGCTTAAAC
GGTACATCGATGGCTACTCCTCATGTTGCAGGTGCAGCAGCCCTTGTTAAACAAAAGAACCCATCTTGGT
CCAATGTACAAATCCGCAATCATCTAAAGAATACGGCAACGAGCTTAGGAAGCACGAACTTGTATGGAAG
CGGACTTGTCAATGCAGAAGCGGCAACACGCTAATCAATAATAATAGGAGCTGTCCCAAAAGGTCATAGA
TAAATGACCTTTTGGGGTGGCTTTTTTACATTTGGATAAAAAAGCACAAAAAAATCGCCTCATCGTTTAA
AATGAAGGTACC

9. Søg efter sekvensen i NR/NT databasen (BLASTN) med default indstillinger.
   
• Er der nogen signifikante hits?
  Svar: ja, der er 6 hits med en e-value på "0.0" (dvs så lille at den bliver afrundet til nul) - og et på 6e-131.
  Det første hit (S48754) har også en coverage på 100% og en identity på 100%.
  
  
• Hvad koder disse gener for?
  Svar: De er alle basiske ("alkaline") serin-protease fra Bacillus -> vores DNA stammer faktisk fra det første entry ("S48754").
  
  
• Bemærk farvekodningen i den grafiske oversigt over hits.
  Svar: 4 (næsten) perfekte hits (bit-score: 2100-2700), et par stykker der er signifikante men korte (bit-score: 75-100, e-value: 1e-10 til 1e-18) og en helt masse junk, der består af korte fragmenter med en dårlig e-value. Altså, en rimelige brat ovegang fra de gode hits til de virkeligt dårlige.
  
  
• Hvordan er fordelingen af kvaliteten af hits: Er det en glidende overgang, eller er det fordelt i ekstremerne?
  Som gennemgået i forelæsningen, er DNA alignment-matricen optimeret til at finde hits der matcher godt - bemærk at mismatches straffes hårdt (Bemærk at NCBI også har en 2/-3 matrice - men effekten er basalt set den samme):
  
  
    a  c  g  t
  a  1
  c -3  1
  g -3 -3  1
  t -3 -3 -3  1
  
  
  
  
  
10. Lad os prøve det samme med sekvensen på protein-niveau.
• Oversæt den længste ORF med VirtualRibosome (tip: husk at vælge alle læserammer), og gem/kopier proteinsekvesen i FASTA format.
  
  Svar: HUSK at bruger OFR finderen!
  
  Eftersom vi ved at det er et fuld-længde transcript, kan vi godt regne med at både START og STOP codon er med. Vi kan derfor vælge "START-codon = strict".
  
  Bemærk: man kan oversætte med både den standard-genetiske kode (tabel 1) eller alternativt med tabel 11 ("The Bacterial and Plant Plastid Code"). Den eneste forskel er at tabel 11 tillader nogle ekstra yderst sjældne START codons, bla. fra bacterielle phager.
  
  
• Brug BLASTP mod NR.
• Findes der nogen konserverede protein-domæner? Dette kan være en vigtig brik i puslespillet med at finde ud af hvad en given proteinsekvens har af funktion.
  
  Svar: Ja, der findes et "Peptidase S8" domæne. Dette kan se på både "pause-siden" mens der søges, og ved at klikke på "Show conserved domains" øverst på resultat siden.
  
  
• Er der nogen signifikante hits? Svar: Ja, ganske mange.
  
  
• Hvilken type enzymer er de bedste hits?
  Svar: Proteaser og Peptidaser (hvilket dækker over det samme). Bemærk at man kan klikke på hits'ne og gå direkte til sekvensen, hvor man kan læse yderligere. Vi kan slå følgende fast: 1) Der er en serin protease, 2) af S8 familien, 3) den er basisk, 4) den kommer fra Bacillus.
  
  
• Hvordan ser det ud med fordelingen af hits (høj -> lav kvalitet)?
  Svar: Fra fuld-længde hits med en e-value på mindre end 1e-169 og gradvis dernedaf (per default ser vi kun de 100 bedste hits).
  
  
• Er protein-alignment'et bedre til at fange breden af similaritet? Svar: Ja.
  
  
11. OK - det var jo lidt snyd at bruge en sekvens, der allerede var tilstede i databasen, så lad os istedet kaste os over følgende sekvens.
    
    Det er her tale om et DNA fragment sekventeret direkte fra DNA ekstraheret fra en jordprøve. Fragmentet der går under det poetiske navn "KLON12" er amplificeret med PCR primere der er særligt designet til at amplificere kernen af gener der koder for en bestem klasse af enzymer. (I parantes bemærket er det er gen-fragment jeg selv har sekventeret - mit speciale handlede om ikke-kultiverbare mikroorganismer fra jord).
    

12. Opgaven er nu at finde ud af hvad denne stump sekvens koder for, ved at bruge de metoder vi netop har gennemgået.
    
    Bemærk: Sekvensen er (mere eller mindre) i GenBank format og NCBI BLAST forventer enten FASTA eller rå sekvens. De fleste webservere (ikke alle!), der også kan acceptere en sekvens i rå format, sørger automatisk for at smide alle mellemrum og tal væk. Dette gælder også for både NCBI BLAST og de fleste af de servere jeg selv har skrevet, fx VirtualRibosome. Man kan altså kopiere de linier ind hvor selve sekvensen står. NB NB NB: Hvis man skal submitte flere sekvenser på een gang, virker dette trick ikke - brug korrekt formaterede FASTA sekvenser istedet.
    
    Tænk over følgende inden I går i gang med søgningen:
    
    • Er der tale om proteinkodende sekvens? Svar: Ja.
    • Hvis ja, kan man så regne med at START og STOP codon er med?
      Svar: Nej - det er et fragment af en gen - det ved vi kun fordi jeg allerede have giver jer en hint om det. Hvis man ikke har den oplysning kan man prøve sig frem.
      
    
    Lad os som udgangspunkt sige at vi ikke tror på noget med en E-værdi på større end ("dårligere end") 10^-10.
    
    Dette kunne være et typisk eksempel på en eksamensopgave. Det er vigtigt at I i jeres svar argumenterer (bare i stikords form, for jeres egen skyld, siden I ikke skal aflevere denne opgave) for jeres strategi for at nå frem til et resultat:
    
    
• Hvilke værktøjer / databaser er relevante at bruge?
  
  Svar: Vi vil i sidste ende gerne bruge BLAST til at søge i de store samlede databaser.
  
  Lad os derfor prøve følgende:
  
  1) BLASTN
  2) Oversætte til protein (fx. VirtualRibosome).
  3) BLASTP
  
  Både med BLASTN og BLASTP ønsker vi at bruge NR databasen - det er den der giver den største bredde. Der er godt at started med en stor, bred database - så kan man evt. gå i dybten i en af de mere specifikke, når man har skudt sig ind på hvad for en type sekvens man har.
  
  
• Hvad er resultaterne af de forskellige typer søgning - hvad virker og hvad virker ikke (og hvorfor)?
  
  Svar:
  
  1) BLASTN - når vi prøver at BLAST'e mod NR får vi slet ikke nogen signifikante hits:
  
  
  gb|CP000159.1|  Salinibacter ruber DSM 13855, complete genome      51.8    0.003
  gb|CP000813.1|  Bacillus pumilus SAFR-032, complete genome         50.0    0.012
  gb|CP000485.1|  Bacillus thuringiensis str. Al Hakam, complete...  48.2    0.042
  gb|CP000001.1|  Bacillus cereus E33L, complete genome              48.2    0.042
  dbj|AB085752.1|  Bacillus sp. KSM-LD1 gene for protease, compl...  48.2    0.042
  ref|XM_001193501.1|  PREDICTED: Strongylocentrotus purpuratus ...  44.6    0.51
  ref|XM_001196149.1|  PREDICTED: Strongylocentrotus purpuratus ...  44.6    0.51
  gb|BT014524.1|  Lycopersicon esculentum clone 133938F, mRNA se...  44.6    0.51
  emb|AM486171.2|  Vitis vinifera contig VV78X025776.9, whole ge...  42.8    1.8 
  emb|AM478115.2|  Vitis vinifera contig VV78X115679.5, whole ge...  42.8    1.8 
  emb|CR392038.9|  Zebrafish DNA sequence from clone DKEY-172J4 ...  42.8    1.8 
  ref|XM_001191953.1|  PREDICTED: Strongylocentrotus purpuratus ...  42.8    1.8 
  dbj|AB214317.1|  Bacillus licheniformis prz gene for keratinaz...  42.8    1.8 
  gb|AY136509.1|  Glycine max subtilisin-like protease C1 gene, ...  42.8    1.8 
  gb|AC133779.21|  Medicago truncatula clone mth2-27n4, complete...  42.8    1.8 
  gb|AF036960.3|  Glycine max subtilisin-like protease C1 mRNA, ...  42.8    1.8 
  ref|XM_779662.1|  PREDICTED: Strongylocentrotus purpuratus sim...  42.8    1.8 
  emb|CU329670.1|  Schizosaccharomyces pombe chromosome I            41.0    6.2 
  dbj|AB305158.1|  Alkalimonas collagenimarina aprI gene for sub...  41.0    6.2 
  emb|CT030005.14|  Zebrafish DNA sequence from clone DKEY-206D1...  41.0    6.2 
  ref|XM_001583403.1|  Trichomonas vaginalis G3 Dynein heavy cha...  41.0    6.2 
  ref|XM_001180101.1|  PREDICTED: Strongylocentrotus purpuratus ...  41.0    6.2 
  ref|XM_785195.2|  PREDICTED: Strongylocentrotus purpuratus sim...  41.0    6.2 
  ref|XR_023110.1|  PREDICTED: Pan troglodytes kinase suppressor...  41.0    6.2 
  ref|XM_001105629.1|  PREDICTED: Macaca mulatta kinase suppress...  41.0    6.2 
  ref|NM_115638.1|  Arabidopsis thaliana unknown protein (AT3G57...  41.0    6.2 
  dbj|AK124344.1|  Homo sapiens cDNA FLJ42353 fis, clone UTERU20...  41.0    6.2 
  ref|XM_388640.1|  Gibberella zeae PH-1 strain PH-1; NRRL 31084...  41.0    6.2 
  gb|AC108679.6|  Homo sapiens 3 BAC RP11-485G4 (Roswell Park Ca...  41.0    6.2 
  gb|AC148994.3|  Medicago truncatula chromosome 7 clone mth2-3c...  41.0    6.2 
  gb|AY126451.1|  Dactylaria parvispora serine protease mRNA, pa...  41.0    6.2 
  
  
  (TIP: Man kan få dette tekst-baserede output ved at klikke på "Reformat these result" i toppen af siden, og så vælge "Plain text").
  
  Der er simpelt hen ikke noget i databasen, der er ligner vores input tilstrækkeligt - derfor finder vi kun de tilfældige resultater. En søgning på DNA niveau er kun egnet til at finde de sekvenser, der ligner virkeligt meget.
  
  2) Virtual ribosome - resultat fra ORF finderen (NB: tabel 1, "Standard Genetic Code" kan også sagtens bruges).
  
  VIRTUAL RIBOSOME
  ----------------
  Translation table: Bacterial and Plant Plastid
  
  >Seq1_rframe1_ORF
  Reading frame: 1
  
      N  G  H  G  T  H  V  A  G  T  V  A  A  V  N  N  N  G  I  G  V  A  G  V  A  G  G  N  G  S 
  5' AACGGGCACGGGACGCATGTAGCTGGAACAGTGGCAGCCGTAAATAATAATGGTATCGGAGTTGCCGGGGTTGCAGGAGGAAACGGCTCT 90
     ..............................))).....................))).................................
  
      T  N  S  G  A  R  L  M  S  T  Q  I  F  N  S  D  G  D  Y  T  N  S  E  T  L  V  Y  R  A  I 
  5' ACCAATAGTGGAGCAAGGTTAATGTCCACACAAATTTTTAATAGTGATGGGGATTATACAAATAGCGAAACTCTTGTGTACAGAGCCATT 180
     .....................>>>.........))).......................................))).........)))
  
      V  Y  G  A  D  N  G  A  V  I  S  Q  N  S  W  G  S  Q  S  L  T  I  K  E  L  Q  K  A  A  I 
  5' GTTTATGGTGCAGATAACGGAGCTGTGATCTCGCAAAATAGCTGGGGTAGTCAGTCTCTGACTATTAAGGAGTTGCAGAAAGCTGCGATC 270
     ........................))))))...........................)))...)))......)))............)))
  
      D  Y  F  I  D  Y  A  G  M  D  E  T  G  E  I  Q  T  G  P  M  R  G  G  I  F  I  A  A  A  G 
  5' GACTATTTCATTGATTATGCAGGAATGGACGAAACAGGAGAAATACAGACAGGCCCTATGAGGGGAGGTATATTTATAGCTGCCGCCGGA 360
     .........)))............>>>...............)))............>>>.........)))...)))............
  
      N  D  N  V  S  T  P  N  M  P  S  A  Y  E  R  V  L  A  V  A  S  M  G  P  D  F  T  K  A  S 
  5' AACGATAACGTTTCCACTCCAAATATGCCTTCAGCTTATGAACGGGTTTTAGCTGTGGCCTCAATGGGACCAGATTTTACTAAGGCAAGC 450
     ........................>>>...........................)))......>>>........................
  
      Y  S  T  F  G  T  W  T  D  I  T  A  P  G  G  D  I  D  K  F  D  L  S  E  Y  G  V  L  S  T 
  5' TATAGCACTTTTGGAACATGGACTGATATTACTGCTCCTGGCGGAGATATTGACAAATTTGATTTGTCAGAATACGGAGTTCTCAGCACT 540
     ...........................)))..................)))............)))........................
  
      Y  A  D  N  Y  Y  A  Y  G  E  G  T  S  M  A  C  P  H  V  A  G  A  A 
  5' TATGCCGATAATTATTATGCTTATGGAGAGGGAACATCCATGGCTTGTCCACATGTCGCCGGCGCCGCC 609
     .......................................>>>...........................
  
  
  
  (TIP: Husk at du kan få sekvensen i FASTA format ved at klikke på "FASTA" linket på resultat-siden).
  
  3) BLASTP
  
  Det første vi ser, er at at der finds et "Peptidase_S8" domæne. Dette er en meget kraftig indikator på funktionen.
  
  Vi får en række signifinkante hits. Hvis vi kigger på de bedste, ser at de alle er serin-proteaser (af S8 familien).
  
  
  ref|NP_809125.1|  protease [Bacteroides thetaiotaomicron VPI-5...   154    4e-36
  ref|ZP_02425799.1|  hypothetical protein ALIPUT_01953 [Alistip...   147    5e-34
  ref|ZP_02437453.1|  hypothetical protein BACSTE_03728 [Bactero...   142    2e-32
  ref|ZP_02424566.1|  hypothetical protein ALIPUT_00684 [Alistip...   137    3e-31
  ref|YP_446403.1|  peptidase families S8 and S53 domain protein...   129    1e-28
  ref|ZP_02063425.1|  hypothetical protein BACOVA_00373 [Bactero...   126    7e-28
  gb|AAB49694.1|  protease [Bacillus sp. PD498]                       117    4e-25
  ref|YP_475789.1|  peptidase, S8A (subtilisin) family [Synechoc...   111    3e-23
  ref|YP_477062.1|  peptidase, S8A (subtilisin) family [Synechoc...   110    6e-23
  ref|NP_693519.1|  thermophilic serine proteinase [Oceanobacill...   107    5e-22
  gb|AAB36499.1|  thermostable alkaline protease [Thermoactinomy...   103    5e-21
  sp|Q45670|THES_BACSJ  Thermophilic serine proteinase precursor...   103    7e-21
  gb|AAK29176.1|  thermophilic alkaline protease [Geobacillus st...   103    8e-21
  ref|ZP_01860436.1|  thermophilic serine proteinase [Bacillus s...   103    9e-21
  pdb|1DBI|A  Chain A, Crystal Structure Of A Thermostable Serin...   102    2e-20
  ref|ZP_01170337.1|  subtilisin-type proteinase [Bacillus sp. N...   101    3e-20
  ref|ZP_02328516.1|  subtilisin-type proteinase [Paenibacillus ...   100    5e-20
  ref|NP_487583.1|  protease [Nostoc sp. PCC 7120] >dbj|BAB75242...   100    5e-20
  ref|ZP_00106975.2|  COG1404: Subtilisin-like serine proteases ...   100    5e-20
  ref|YP_324024.1|  Peptidase S8 and S53, subtilisin, kexin, sed...   100    6e-20
  ref|ZP_01630730.1|  protease [Nodularia spumigena CCY9414] >gb...   100    8e-20
  ref|YP_175254.1|  serine protease [Bacillus clausii KSM-K16] >...  99.8    9e-20
  ref|NP_242357.1|  protease [Bacillus halodurans C-125] >dbj|BA...  99.8    9e-20
  
  
  (BEMÆRK: "Hypotetical" proteiner er ikke nær så stærk evidens som et velbeskrevet protein).
  
  Lad os kigge lidt nærme på det bedste hit (NP_809125):
  
  
  >ref|NP_809125.1| protease [Bacteroides thetaiotaomicron VPI-5482]
   gb|AAO75319.1| protease [Bacteroides thetaiotaomicron VPI-5482]
  Length=695
  
   Score =  154 bits (389),  Expect = 4e-36, Method: Compositional matrix adjust.
   Identities = 98/238 (41%), Positives = 134/238 (56%), Gaps = 45/238 (18%)
  
  Query  2    GHGTHVAGTVAAVNNNGIGVAGVAGGNGSTNSGARLMSTQIFNSDGDYT-NSETLVYRAI  60
              GHGTHVAGT+AA+NNNG GV G+AGG+GS NSG ++MS Q+F  +   T ++E    RAI
  Sbjct  295  GHGTHVAGTIAAMNNNGEGVCGIAGGDGSKNSGVKIMSCQVFAGEAGVTLDAEA---RAI  351
  
  Query  61   VYGADNGAVISQNSWGSQS-------------------LTIKELQKAAIDYFIDYAGMDE  101
               Y ADNGAVI Q SWG  S                     +  L+K A+DYFI+ AG  
  Sbjct  352  KYAADNGAVILQCSWGYNSSLANLIEGYTPGPGSEEEWEKLYPLEKDALDYFINNAGS--  409
  
  Query  102  TGEIQTGPMRGGIFIAAAGNDNVSTPNMPSAYERVLAVASMGPDFTKASYSTFGTWTDIT  161
                    G + GG+ I A+GN+       P+AY + ++V+++  DFT ASYS +G    I+
  Sbjct  410  ----PNGVIDGGLAIFASGNEYAGMAAFPAAYSKCISVSAVAADFTPASYSNYGKEVTIS  465
  
  Query  162  APGGDIDKFD-------------LSEYGVLSTYADN---YYAYGEGTSMACPHVAGAA  203
              APGGD + ++             +    +LST+  N    Y + +GTSMACPHV+G A
  Sbjct  466  APGGDTEYYNPVGQDDPEGWEEGIHSGSILSTWIQNGNATYGFMDGTSMACPHVSGVA  523
  
  
  Bemærk at selv om det ikke er et perfekt hit (vores input-sekvens findes jo netop ikke i databasen), så ser det ganske fornuftigt ud: alignment'et er langt med en Identity på 41% og en Similarity på 56%. Sammenholdt med det faktum at alle de bedste hits er serin-proteaser, har vi en meget stærk indikation på hvilken type af enzym vi har med at gøre.
  
  
• Endelig konklusion: Hvad er KLON12?
  Svar: En serin-protease, fra S8 familien.
  
  

Del 3 - BLAST og genomer

13. Lad os undersøge hvordan histonerne (ansvarlige for at pakke DNA) er beslægtede mellem gær og menneske (evolutionær afstand: 1-1,5 milliard år).
    
    Start med at slå HTA2 genet op i SGD (hint: der er en Quick Search box øverst på siden). Bemærk at der er en kort forklaring om hvad genet og dets protein-produkt gør. Hvad står der her om dets forhold til genet "HTA1"?
    Svar: HTA2 og HTA1 er næsten ens ("nearly identical").
    
    Slå nu protein-sekvensen op (se under "Protein Info & Structure" ude til højre). Hold vinduet åbent eller gem sekvensen, vi skal bruge den om lidt.
    
    
14. Gå nu til NCBI's genom side. Læg mærke til at der ude til højre findes en liste af organismer (under "Organism-specific") hvor de informationer der findes er vist med bogstavkode ("B" = BLAST). Hvis man her vælger BLAST for en specifik organisme, får man (ikke overraskende) mulighed for at søge i dette specifikke genom (og de oversatte proteinsekvenser) med BLAST.
    
    Inden vi søger i det humane genom, så lad os undersøge hvad vi kan finde hvis vi søget i gærs eget genom. Vælg derfor "B" under "Yeast".
• Vælg "RefSeq protein" som database.
  
• Vælg "BLASTP" som program.
• Søg efter HTA2 protein-sekvensen.
  
  
• Hvor mange pålidelige hits findes der? Svar: 3
• Giver det mening ud fra det der stod under beskrivelsen af HTA2?
  Svar: Ja - HTA1 og HTA2 er meget ens (kun 2 aminosyrer er forskellige). Som bonus finder vi også "H2AZ" - en anden variant af HTA2 der noget mere divergeret i sekvens.
  
  
15. Næste trin er at søge i det (oversatte) humane genom. Gå ind på BLAST siden for det humane genom.
• Vælg "RefSeq protein" som database.
• Bemærk at der er flere muligheder for valg af database end hos gær - det skyldes simpelthen at man har dårligere styr på det humane genom, og at man som bruger kan have brug at at vælge aktivt mellem forskellige fortolkninger.
• Vælg "BLASTP" som program og søg efter homologer til HTA2 hos mennesket.
  
  HITS:
  
  
  ref|NP_003503.1|  H2A histone family, member L [Homo sapiens]       191    1e-49
  ref|NP_003507.1|  H2A histone family, member O [Homo sapiens] ...   190    2e-49
  ref|NP_254280.1|  histone H2a [Homo sapiens]                        190    3e-49
  ref|NP_003500.1|  H2A histone family, member C [Homo sapiens] ...   190    3e-49
  ref|NP_542163.1|  H2A histone family member [Homo sapiens]          190    3e-49
  ref|NP_734466.1|  histone H2A [Homo sapiens]                        189    4e-49
  ref|NP_778235.1|  histone H2A [Homo sapiens]                        189    5e-49
  ref|NP_002096.1|  H2A histone family, member X [Homo sapiens]       189    6e-49
  ref|NP_003508.1|  H2A histone family, member Q [Homo sapiens]       189    6e-49
  ref|NP_066409.1|  histone 1, H2ad [Homo sapiens]                    189    6e-49
  ref|NP_066390.1|  H2A histone family, member A [Homo sapiens] ...   189    6e-49
  ref|NP_066544.1|  H2A histone family, member E [Homo sapiens]       188    8e-49
  ref|NP_808760.1|  H2A histone family, member J isoform 2 [Homo sa   188    1e-48
  ref|NP_060737.1|  H2A histone family, member J isoform 1 [Homo sa   184    2e-47
  ref|NP_061119.1|  core histone macroH2A2.2 [Homo sapiens]           150    2e-37
  ref|NP_613258.2|  H2A histone family, member Y isoform 3 [Homo sa   149    7e-37
  ref|NP_613075.1|  H2A histone family, member Y isoform 1 [Homo sa   149    7e-37
  ref|NP_004884.1|  H2A histone family, member Y isoform 2 [Homo...   149    7e-37
  ref|NP_036544.1|  H2A histone family, member V isoform 1 [Homo sa   130    3e-31
  ref|NP_002097.1|  H2A histone family, member Z [Homo sapiens]       129    6e-31
  ref|NP_619541.1|  H2A histone family, member V isoform 2 [Homo sa   111    1e-25
  ref|NP_958844.1|  H2A histone family, member V isoform 3 [Homo sa   108    1e-24
  ref|XP_946222.1|  PREDICTED: similar to Histone H2a [Homo sapiens   104    2e-23
  ref|NP_542451.1|  H2A histone family, member B3 [Homo sapiens]...  87.4    3e-18
  ref|NP_001017990.1|  H2A histone family, member B1 [Homo sapiens]  84.7    2e-17
  ref|XP_934613.2|  PREDICTED: similar to testis-specific histon...  77.4    3e-15
  ref|NP_958925.1|  H2A histone family, member V isoform 5 [Homo sa  77.4    3e-15
  ref|XP_944049.2|  PREDICTED: similar to H2A histone family, me...  65.5    1e-11
  ref|XP_370909.1|  PREDICTED: similar to H2A histone family, me...  65.5    1e-11
  ref|NP_958924.1|  H2A histone family, member V isoform 4 [Homo sa  63.9    3e-11
  ref|NP_005624.2|  son of sevenless homolog 1 [Homo sapiens]        47.8    2e-06
  ref|NP_001018082.1|  BTB (POZ) domain containing 11 isoform 3 [Ho  42.7    7e-05
  ref|NP_665803.1|  ankyrin repeat and BTB (POZ) domain containing   37.0    0.004
  
  
  
  
• Hvor mange pålidelige hit findes der?
  Svar: Lad os ligesom før kræve at hits skal have en e-value på 1e-10 eller bedre (mindre) for at vi regner det for pålideligt. Dette passer også fint med, at disse hits er enige om at der er tale om en histon.
  
  Som udgangspunkt er der derfor 30 gode hits. Man kan dog godt argumentere for at man skal passe på hits der kun har en forudsagt funktion ("PREDICTED") - dem har vi fire af. Der er også tre hits der hedder helt præcist det samme ("histone H2A [Homo sapiens]"), om her kunne men evt. gå ind på link'et til sekvensen og læse hvad der står - her kan man se, at der faktisk er mindre forskelle i sekvens, og at de er stammer fra forskellige gener (/gene="HIST3H2A", /gene="HIST1H2AA", /gene="HIST2H2AB).
  
  Men lad os arbejde videre med alle 30 gode hits i de næste analyser.
  
  
• Fra hvor mange gener stammer disse?
  Svar: Trick'et er her af en del af hits'ne er fra isoformer af det sammen protein (se næste spørgsmål). Et protein med to isofomer stammer stadig kun fra en enkelt gen.
  
  Lad os fx kigge på følgende hits:
  
  ref|NP_613075.1|  H2A histone family, member Y isoform 1 [Homo sa   149    7e-37
  ref|NP_004884.1|  H2A histone family, member Y isoform 2 [Homo...   149    7e-37
  
  Hvis man gå ind og læser sekvens-entry'et for disse hits kan man faktisk direkte se, at de stammer fra sammen gen (/gene="H2AFY", /gene="H2AFY").
  
  Vi skal altså slå hits sammen, der dækker over isoformer kodet fra sammen gen. Der er 5 "V" isoformer, 3 "Y" isoformer og 2 "J" isoformer, så det bliver til 30 - 4 -2 -1 = 23 gener.
  
  
• Er der nogen af disse gener der har iso-former (altså forskellige protein-produkter fra samme gen, typisk pga. alternativ splejsning)?
  Svar: Se ovenstående.
  
  
• Bemærk: Noget af sekvensen er markeret i gråt - det er lav-komplekse områder (dvs. at de indeholder brudstykker af ofte forekommende sekvens, som ikke er så meget værd for BLAST algorithmen). Man kan vælge hvordan BLAST skal tage sig af disse område på den første side. Hvis man vælger "low complexity" under "filter", bliver disse områder slet ikke taget med når de to indledende hits skal findes (default i almindelig BLASTP!).
  
  
16. Prøv at kører søgning igen (i et nyt vindue eller tab) og denne gang slå lav-kompleksitets filteret til.
• Hvordan ser det nu ud med længderne af de alignments man får?
  Svar: De bliver kortere i de tilfælde, hvor det lav-komplekse område ligger i starten eller i slutningen af hit'et.
  
  
17. Afsluttende bemærkning: Vi har nu med BLAST fundet en række homologe gener (eller rettere protein-produkter). Hvis man ønsker at grave dybere, vil det næste logiske trin være at trække den komplette sekvens for kandidat proteinerne ud af databasen (og evt. ignorer isoformerne i første omgang), for derefter at analysere dem med enten parvis alignment (for at finde specifikke ligheder/forskelle mellem HTA2 og kandidaten) eller multiple alignment (kommende emne) for at kunne analysere kandidaternes indbyrdes slægskabsforhold.
    
    BLAST kan altså også være en god måde at begynde at bygge et datasæt af sekvenser ud fra en kendt sekvens (enten fra en anden organisme eller simpelthen bare et andet gen fra samme organisme).