Events News Research CBS CBS Publications Bioinformatics
Staff Contact About Internal CBS CBS Other

Ø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:

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") .

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

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

  1. Åben BLASTN siden hos NCBI og indsæt dine tre sekvenser (hint: FASTA format). Vælg "NR" som database. Lad resten stå med default settings og tryk submit - bemærk at når kørslen først er sat igang, skal man trykke på "format" knappen for at se resultaterne.

    NR står for "Non Redundancy" og er en samling af alle ikke-redundante sekvenser fra GenBank og genom-sekventerings projekterne.

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

    • Er der nogen sekvenser de ligner vores tilfældigt genererede?
    • Hvor lange er de hits der findes?
    • Hvor mange % identity?
    • I hvilket område ligger bit-scores ("max score")?
    • I hvilket område ligger E-value?

    • Giver disse resultater nogen biologisk mening?

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

    • Hvor lange er de hits der findes?
    • Hvor mange % identity?
    • I hvilket område ligger bit-scores ("max score")?
    • I hvilket område ligger E-value?

    • De humane sekvenser findes også i NR - blev de samme sekvenser fundet i søgningen mod NR databasen? Hvorfor / hvorfor ikke?

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

  5. 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
 1 : A
 2 : R
 3 : N
 4 : D
 5 : C
 6 : Q
 7 : E
 8 : G
 9 : H
10 : I
11 : L
12 : K
13 : M
14 : F
15 : P
16 : S
17 : T
18 : W
19 : Y
20 : V
  1. Åbn BLASTP siden hos NCBI, og sæt sekvenserne ind. Vælg igen "NR" databasen (denne gang består den af oversætte CDS'er + SwissProt mm). Vælg BLOSUM62 som matrice, og sæt "Expect" cut-off til 1000 (default: 10). Submit derefter søgningen.

    • Hvor stor er databasen?

    • Hvor lange er hits'ne typisk og har de gaps?
    • Hvor mange % identity?
    • I hvilket område ligger E-value?
    • 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)?

    • Hvis vi havde brugt default cut-off'et for E-value på 10, ville der så være fundet nogen sekvenser?
    • Fanger vi mest junk ved at bruge DNA eller peptid i vores søgninger?

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

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

  1. Søg efter sekvensen i NR databasen (BLASTN) med default indstillinger.
    • Er der nogen signifikante hits?
    • Hvad koder disse gener for?
    • Bemærk farvekodningen i den grafiske oversigt over hits.
    • Hvordan er fordelingen af kvaliteten af hits: Er det en glidende overgang, eller er det fordelt i ekstremerne?

  2. 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.
    • Brug BLASTP mod NR.
      • Findes der nogen konserverede protein-domæner? (Viser på den side hvor I også skal trykke "format") - dette kan være en vigtig brik i puslespillet med at finde ud af hvad en given proteinsekvens har af funktion.

      • Er der nogen signifikante hits?
      • Hvilken type enzymer er de bedste hits?
      • Hvordan ser det ud med fordelingen af hits (høj -> lav kvalitet)?
      • Er protein-alignment'et bedre til at fange breden af similaritet?

  3. 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 peotiske 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).
LOCUS       KLON12.DNA    609 BP DS-DNA             UPDATED   06/14/98
DEFINITION  UWGCG file capture
ACCESSION   -
KEYWORDS    -
SOURCE      -
COMMENT     Non-sequence data from original file:
BASE COUNT      174 A    116 C    162 G    157 T      0 OTHER
ORIGIN      ?
    klon12.dna Length: 609   Jun 13, 1998 - 03:39 PM   Check: 6014 ..
        1 AACGGGCACG GGACGCATGT AGCTGGAACA GTGGCAGCCG TAAATAATAA TGGTATCGGA
       61 GTTGCCGGGG TTGCAGGAGG AAACGGCTCT ACCAATAGTG GAGCAAGGTT AATGTCCACA
      121 CAAATTTTTA ATAGTGATGG GGATTATACA AATAGCGAAA CTCTTGTGTA CAGAGCCATT
      181 GTTTATGGTG CAGATAACGG AGCTGTGATC TCGCAAAATA GCTGGGGTAG TCAGTCTCTG
      241 ACTATTAAGG AGTTGCAGAA AGCTGCGATC GACTATTTCA TTGATTATGC AGGAATGGAC
      301 GAAACAGGAG AAATACAGAC AGGCCCTATG AGGGGAGGTA TATTTATAGC TGCCGCCGGA
      361 AACGATAACG TTTCCACTCC AAATATGCCT TCAGCTTATG AACGGGTTTT AGCTGTGGCC
      421 TCAATGGGAC CAGATTTTAC TAAGGCAAGC TATAGCACTT TTGGAACATG GACTGATATT
      481 ACTGCTCCTG GCGGAGATAT TGACAAATTT GATTTGTCAG AATACGGAGT TCTCAGCACT
      541 TATGCCGATA ATTATTATGC TTATGGAGAG GGAACATCCA TGGCTTGTCC ACATGTCGCC
      601 GGCGCCGCC
//


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

    Tænk over følgende inden I går i gang med søgningen:
    • Er der tale om proteinkodende sekvens?
    • Hvis ja, kan man så regne med at START og STOP codon er med?

    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 er relevante at bruge?






    • Hvad er resultaterne af de forskellige typer søgning - hvad virker og hvad virker ikke (og hvorfor)?







    • Endelig konklusion: Hvad er KLON12?

Del 3 - BLAST og genomer

Nu har vi indtil videre brugt BLAST til at søge i de store brede databaser (GenBank + andre DNA-databaser, UniProt + andre proteindatabaser). Lad os runde øvelsen af med at kigge på hvordan man bruger BLAST på nogle mere smalle databaser, nemlig til at søge i specifikke genomdatabaser.

Dette vil typisk være nyttigt hvis man har en sekvens med kendt funktion fra en bestemt organisme (fx. et gen der styrer cell-deling hos gær, Saccharomyces cerevisiae), og gerne vil finde de tilsvarende gener hos en anden organisme (fx. menneske - gener der styrer cell-deling er ofte indvolverede i kræft).

Som I sikkert har lagt mærke til har det hos NCBI været muligt at søge i både musen og menneskets genom med BLAST. Der har også været mulighed for kun at ville se resultater fra en bestemt organisme.

Der findes en række andre servere der specifikt handler om en bestemt organismes genom. For eksempel har gær-genomet "hjemme" i Saccharomyces Genome Database (SGD - www.yeastgenome.org) - bemærk iøvrigt at også denne side tilbyder BLAST som en måde at søge i databasen på. BLAST er ganske enkelt standard-værktøjet til at finde sekvenser ud fra andre sekvenser.

Vi skal her nøjes med at kigge på de genom-ressourcer der ligger hos NCBI (med en enkelt afstikker til SGD):

http://www.ncbi.nlm.nih.gov/Genomes/
  1. 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"?

    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.

  2. 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?
    • Giver det mening ud fra det der stod under beskrivelsen af HTA2?

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

    • Hvor mange pålidelige hit findes der?
    • Fra hvor mange gener stammer disse?
    • Er der nogen af disse gener der har iso-former (altså forskellige protein-produkter fra samme gen, typisk pga. alternativ splejsning)?

    • 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!).

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

  5. 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).