I en database er det 1000 000 sekvenser. Vi har en query-sekvens, og på en eller annen måte får vi kartlagt at databasen inneholder i alt 1500 homologer til vår query. Vi gjør et søk med den nye super-metoden PLAST og får 1000 treff. Av disse treffene viser det seg etterhvert at 950 er homologer til vår query-sekvens.
Hva er sensitiviteten (sensitivity eller recall) til PLAST i dette søket?
Hva er spesifisiteten (specificity) til PLAST her?
Hva blir presisjonen (precision)?
HINT: Det er lurt å sette opp en 2x2 tabell som lister hvor mange vi har av Sanne Positive, Falske Positive, Sanne Negative og Falske Negative.
Hvis du er i tvil om begrepene, les litt på denne wikipedia-siden
Bruk den enkle scoringsregelen 1 for likhet, 0 for ulikhet og 0 for indel. Regn ut SP-score for sammenstillingen
CCGAT
GCCAT
-CGAT
CAGTTNotér score for hver enkelt kolonne/posisjon. Hva er den maksimale score-verdien vi kan få i en kolonne med denne scorings-regelen? Hva er nest største mulige score? Hvilke score-verdier kan en kolonne egentlig få med denne regelen?
Finn konsensus-sekvensen, det vil si den sekvensen som framkommer ved majoritets-avstemming i hver posisjon. Graden av konservering i en kolonne angir vi ofte med andelen sekvenser som har konsensus-symbolet i den posisjonen. Finn disse verdiene for sammenstillingen over.
Hvordan henger dette sammen med SP-scoren for en posisjon?
Det sies at like oppstrøms av kodende gener hos bakterier finner vi et konservert motiv som kalles Ribosomalt Binde Sete (RBS). Dersom vi ser på del-sekvensen som består av de 30 basene oppstrøms start-kodonet så bør vi finne RBS motivet der.
I filen (høyre-klikk, lagre som) upstream30_Ecoli.fasta finner du 8 slike sekvenser fra et genom av Escherichia coli. Hva slags motiv forventer vi å finne i disse sekvensene (RBS for E. coli, søk på Wikipedia hvis du ikke vet det fra før)?
Vi vil nå prøve fire ulike MSA-verktøy for å se om de lager sammenstillinger der dette motivet framtrer.
Vi bruker EMBL-EBI sine websider denne gangen. Gå til siden for T-Coffee hos EBI, lim inn sekvensene (med header-tekster) og kjør i vei med default parametre. Husk å velge DNA i første pop-up meny. Velg output på ClustalW format.
Gjør deretter tilsvarende for Muscle hos EBI
Gjør deretter tilsvarende for Mafft hos EBI
Gjør deretter tilsvarende for Clustal hos EBI
Er det noen som lager en MSA der motivet trer tydelig fram?
Hva kan være grunnen til at dette viste seg litt vanskelig?
Vi skal sammenstille fire sekvenser A, B, C og D. Utfra noen parvise sammenligninger av disse sekvensene har vi en litt grov avstandstabell som ser slik ut:
A B C D
A 0 2 6 1
B 2 0 7 3
C 6 7 0 4
D 1 3 4 0En progressiv sammenstillings-metode vil bygge den multiple sammenstillingen ved å legge inn sekvensene i en bestemt rekkefølge. Angi rekkefølgen på hvordan disse sekvensene sammenstilles ved å lage et guide-tre med UPGMA-metoden (se f.eks. Figure 8.3 på side 279 i Zvelebil&Baum).
I bioinformatikk kommer vi borti problemer der selv de kraftigste datamaskiner sliter med å komme i mål. La oss se litt mer på hvor lett tilsynelatende små problemer fort kan vise seg å være veldig store. Vi vil også her gjøre bruk av enkel lineær regressjon i R, som vi kan bruke til å lage prognoser.
Følgende historie er sann: I forbindelse med et prosjekt (utenfor NMBU) for en tid tilbake ble det snakk om å samle inn alle del-sekvenser av lengde 25 fra mange tusen bakterie-genomer. Det ble da foreslått (av noen uten BIN210) at man istedenfor å søke gjennom alle genomene skulle lage seg en database bestående av alle mulige DNA-sekvenser av lengde 25, og så kjøre et søk mot denne databasen, f.eks. med BLAST.
OK, lag et program som genererer alle mulige DNA-sekvenser av lengde 25, og lagre dette som en FASTA-fil. Hvor vanskelig kan det være?
Her er en funksjon i R som genererer alle DNA-sekvenser av lengde K. Lag deg et nytt R-script og kopier inn følgende kode:
Kmers <- function(K, iter = 1, alphabet = c("A","C","G","T")){
require(stringr)
if(iter < K){
w <- Kmers(K, iter + 1, alphabet)
w <- str_c(rep(w, each = length(alphabet)),
rep(alphabet, times = length(alphabet)))
} else {
w <- alphabet
}
return(sort(w))
}Etter å ha kjørt scriptet, skal funksjonen Kmers() nå være klar til å brukes. Kjør følgende kode i Console-vinduet, for å teste at det virker:
Kmers(1) # skal lage alle DNA-sekvenser av lengde 1
Kmers(2) # skal lage alle DNA-sekvenser av lengde 2
Kmers(3) # skal lage alle DNA-sekvenser av lengde 3NB! Ikke prøv å kjøre Kmers(25)! Vi skal straks se hvorfor…
La oss nå beregne
vi bruker på å generere alle DNA-sekvenser av lengde 1 opp til 12. Kopier inn følgende kode i scriptet ditt:
tbl <- tibble(Lengde = 1:12, # DNA-sekvenser av lengder 1,...,12
Sekunder = 0,
Bytes = 0)
for(k in 1:nrow(tbl)){
cat("Finner alle DNA sekvenser av lengde", tbl$Lengde[k], "\n")
tbl$Sekunder[k] <- system.time(v <- Kmers(K = tbl$Lengde[k]))[3]
tbl$Bytes[k] <- as.numeric(object.size(v))
}Her måler vi tiden i sekunder, og mengden lagringsplass i bytes. Kjør scriptet. Spesielt tiden det tar vil variere noe fra PC til PC, og også fra gang til gang. Mengden lagringsplass blir ganske stor for lengde 12, så hvis dette er for drøyt for din PC kan du stoppe ved 11 eller 10.
Åpne tabellen i Viewer i RStudio, og ta en titt. Det vil fort være slik at for de minste lengdene tar det 0 tid. Legg merke til for hvilke lengder tidsforbruket er stabilt større enn 0, hos meg var det fra lengde 8 og oppover.
La oss tilpasse en regressjons-modell til disse dataene, slik at den beskriver hvordan tiden avhenger av lengden. Vi log-transformerer tiden her, og derfor kan vi ikke bruke data der tiden er 0. Vi bruker subset til å velge kun de radene fra data-tabellen som har tid større enn 0:
modell.tid <- lm(log10(Sekunder)~Lengde, data = tbl, subset = 8:12)Her brukte jeg subset = 8:12 fordi radene 8 til 12 har positive verdier i min tabell, du tilpasser dette til dine tall.
Vi kan så predikere antall sekunder det tar å generere alle DNA-sekvenser av ymse lengder. La oss først prøve de lengdene vi allerede har sett:
sekunder.predikert <- 10^(predict(modell.tid, newdata = tibble(Lengde = 1:12)))Legg merke til at siden vår modell har log10 av sekunder som respons, så må vi ta 10 opphøyd i den predikerte verdien for å få predikert antall sekunder (ikke log-sekunder).
Sjekk at det predikerte antall sekunder er noenlunde i samsvar de observerte i tbl, ihvertfall av samme størrelse-orden.
Prediker så antall sekunder for lengde 25. Omtrent hvor lang tid vil det ta å generere alle DNA-sekvenser av lengde 25?
Gjør så eksakt samme analyse for lagringsplass (Bytes), men her trenger du ikke bruke subset, for ingen verdier er 0. Hvor mye lagringsplass trenger vi? Hvor mye lagringsplass har du på din PC?
Sist uke snakket vi litt om store-O notasjonen, som typisk brukes for å beskrive hvordan tidsforbruket til en algoritme øker når størrelsen på problemet vokser. Beskriv med store-O notasjon hvordan både tidsforbruk og lagringsplass øker med lengden på DNA-sekvensene.
Google begrepet “NP-komplett”. Hva er dette? Problemet over er av typen NP-komplette.
Alle ‘kommandoer’ vi bruker i R er funksjoner. Når vi installerer og laster inn en R-pakke, så får vi tilgjengelig en rekke nye funksjoner. Vi kan også lage oss egne funksjoner, noe vi så et eksempel på over med funksjonen Kmers(). Vi skal ikke dvele veldig mye ved funksjoner her i BIN210, men ved å skjønne det mest grunnleggende blir det enklere å bruke funksjoner, og å lese Help-filene.
En funksjon i R er alltid bygget på følgende skjelett:
funksjonsnavn <- function(input1, input2, ...){
# her kommer koden, som er funksjonens innmat
return(et.eller.annet)
}Funksjoner har et navn som vi velger selv. Nøkkelordet function viser at vi nå spesifiserer en funksjon. Deretter kommer parenteser, der vi lister alle input (argumenter) til funksjonen. Så følger et par krøll-parenteser, og mellom disse står selve funksjonens kode, det som utføres av denne funksjonen. Dette avsluttes alltid med nøkkelordet return og parenteser, og inni disse angir vi objektet som funksjonen gir som output. Dette objektet må da ha blitt opprettet og fått en fornuftig innhold i koden foran.
Dersom vi ser på funksjonen Kmers vi laget over, så har den 3 inputs. Den første, K, er det vi kaller en essensiell input. Det betyr at denne må vi gi en verdi når vi bruker funksjonen. Når vi kjørte funksjonen med koden Kmers(1) så gir vi altså verdien 1 til input K. Vi kunne (og burde?) ha skrevet Kmers(K = 1) for å understreke dette.
Men funksjonen Kmers() har også to inputs til, iter og alphabet. Begge disse er det vi kaller opsjoner eller valgfri input fordi de blir gitt et innhold i definisjonen av funksjonen. Dette gjør at vi ikke trenger å gi disse verdier når vi bruker funksjonen, siden de allerede har fått default innhold. Men, vi kan gi disse nye verdier hvis vi vil. Prøv å kjøre Kmers(3, alphabet = c("a", "b")) og se at funksjonen nå endrer oppførsel i forhold til Kmers(3).
Koden inne i funksjonen kan være alt mulig, avhengig av hva denne ‘maskinen’ skal gjøre. Uansett, så vil denne koden produsere et eller annet resultat. Dette lagres da typisk i et eller annet objekt, og enhver funksjon avsluttes med return() der det som angis inni parentesene er det objektet som inneholder funksjonens output (vel, det hender return() utelates, da er det siste beregnede objekt som retuneres).
Poenget med funksjoner er å ‘pakke inn’ biter av kode som utfører noe vi gjerne vil ha gjort flere ganger i litt ulike sammenhenger. Tenk på funksjonen mean() som beregner gjennomsnitt av en vektor med tall. Uten denne funksjonen måtte vi ha kodet formelen for gjennomsnitt hver gang vi skal beregne noe slikt. Nå er dette kodet en gang for alle i denne funksjonen. Vi pakker derfor kode inn i funksjoner når vi innser at dette er noe vi vil komme til å bruke om igjen flere ganger, og gjerne i litt ulike sammenhenger. Det er derfor vanlig at vi lagrer funksjoner på egne filer, og ikke som en del av en script. Vi så dette i uke 5, da vi lastet inn funksjoner fra fila localign.R. Ta en titt i denne fila, og se at funksjonene følger samme oppskrift som nevnt over.
Når du leser en Help-fil for en funksjon i R, så følger disse en standard oppskrift for utseende. Ta en titt på Help-fila for mean(), skriv ?mean og return i Console-vinduet. Under Usage angis hvordan man bruker funksjonen. Vi ser at input x er essensiell, mens trim og na.rm er opsjoner. Under Arguments finner du en beskrivelse av hva de ulike input skal være. Under Value beskrives hva funksjonen gir som output. Help-fila kan også ha andre avsnitt, det varierer noe fra funksjon til funksjon.
Forrige uke sanket vi en del DNA-sekvenser ved å ta utgangspunkt i en sekvens fra huggorm, og brukte BLAST til å finne tilsvarende sekvenser hos andre arter, en fra hver art. Dersom du mangler en slik fil kan du bruke denne (høyre-klikk og lagre hos deg).
Query-sekvensen vi søkte i BLAST med er en bit av et kodende gen (‘partial cds’). Siden BLAST gir oss lokale sammenstillinger må vi forvente at de sekvensene vi lastet ned også er slike biter av et kodende gen. Ikke så sjelden vil vi gjerne translatere slike kodende gener til proteinet, dvs. oversette ett og ett kodon (triplett av baser) til en aminosyre. Vi vil nå se litt på dette, og se hvordan vi kan bruke litt koding for å komme rundt problemene som oppstår.
translate()Lag et script, og les inn fila med huggorm-sekvensene med readFasta() fra microseq-pakken, og lagre dette i en tabell kalt vipera.tbl.
I microseq-pakken er det en funksjon som heter translate(), og som gjør en translasjon fra DNA til protein. Lag en ny kolonne som heter Protein, og bruk mutate() og tranlsate() til å fylle denne med translaterte sekvenser.
Ta en titt på de translaterte sekvensene. Du vil trolig finne at de inneholder symbolet *, gjerne flere ganger. Eller gjør et kjapt tekst-søk etter * i sekvensene i Console-vinduet: str_detect(vipera.tbl$Protein, "\\*"). Hos meg fant str_detect dette symbolet i alle sekvensene!
Symbolet * betyr (translatert) stopp-kodon. Det er vanligvis kodonene TAA, TGA eller TAG. Disse skal jo ikke forekomme inne i en kodende sekvens! Vel, det er faktisk to grunner til at vi ser slike her. Den første er at dette er mitokondrisk DNA. DNA fra mitokondriene følger en litt annen genetisk kode, nemlig at kodonet TGA oversettes til W som er aminosyren tryptofan, ikke til stopp eller *. Heldigvis kan funksjonen translate() ordne dette. Den kan også bruke translasjons-tabell nummer 4 (default er 11). Dette er ikke eksakt det samme som brukes ved mitokondrisk DNA, men gjør samme oversettelse av stopp-kodonene.
Kjør koden på nytt, men endre translasjons-tabell til 4. Se Help-filen for translate() og prøv å få til dette på egenhånd (et av argumentene er en opsjon for dette).
Kjør et nytt søk med str_detect() etter * blant sekvensene. Forhåpentligvis ble det færre sekvenser som nå har dette symbolet. Men, noen har det kanskje fortsatt?
Vel, stopp-kodon opptrer lett dersom vi translaterer i feil leseramme. Dersom vi har et helt gen, så starter vi på første kodon, start-kodonet, og translaterer kodon for kodon. Men, her har vi bare biter av genet, og da starter vi kanskje inni et kodon. Det betyr at 1 eller 2 baser i starten av sekvensen egentlig tilhører forrige kodon, og at vi må starte translasjonen på base 2 eller 3 for at leserammen skal bli riktig. Hvordan kan vi fikse dette?
Her er en enkel kort sekvens for å demonstrere: x <- "AGATAGGTAACGAT". Dersom vi starter på første base, så er kodon 2 et stopp-kodon (TAG). Dersom vi starter på base 2, så blir kodon nummer 3 et stopp-kodon (TAA). Dersom vi starter på base 3 er det ingen stopp-kodoner.
Når vi har en bit av et gen, og ikke vet leserammen, så vil vi typisk prøve oss fram, og se om en av leserammene blir uten stopp-kodon, og da gå for den.
Hvordan gjør vi dette i R? Hvordan trimmer vi vekk symboler i en sekvens (=tekst)? Med funksjonen str_sub(). Vi har tidligere sett flere funksjoner i str_ familien, dette er enda en. Kikk i Help-filen for str_sub, den har følgende 3 input:
string)start)end)Posisjonene vi angir i start og end kan også telles bakfra! Hvis vi skriver str_sub(x, 1, -1) så betyr det at vi fra teksten x vil ha del-teksten som starter på posisjon 1 og ender på siste posisjon. Dette gir oss hele x, men hvis vi bruker str_sub(x, 2, -1) så får vi x minus første symbol.
Prøv å kjøre følgende koder direkte i Console-vinduet:
translate(str_sub(x, start = 1, end = -1)) # har stopp-kodon
translate(str_sub(x, start = 2, end = -1)) # har også stopp
translate(str_sub(x, start = 3, end = -1)) # har ingen stoppVi kan selvsagt også trimme vekk symboler i slutten av sekvensen, men det trenger vi ikke foreløpig. Leserammen styres kun av hvor i sekvensen vi starter translasjonen.
if_else()Vi har sett hvordan vi kan trimme starten av sekvenser for å skifte leseramme. Men, vi kan jo ikke gjøre dette på alle sekvenser, da ødelegger vi jo de som allerede har riktig leseramme. Vi trenger en betingelse.
Betingelser forekommer i alle programmerings-språk, og her skal vi se på en variant i R, nemlig funksjonen if_else(). En betingelse betyr at noe utføres dersom en eller annen test er sann, og noe annet utføres hvis testen er usann. Funksjonen if_else() tar typisk 3 input:
TRUEeller FALSE (logical)TRUE så tilordnes det som kommer i input 2FALSE, så tilordnes det som kommer i input 3Et lite eksempel:
x <- 1:5
y <- if_else(x > 3, "Tittei", "Heisann")Først fylles x med heltallene fra 1 til 5. Så fylles y med tekster, men to ulike avhengig av hva slags verdi som ligger i de ulike elementene av x. De 3 første elementene av x er ikke større enn 3, og dermed er testen x > 3 usann eller FALSE i disse tilfellene. Ergo vil de tre første elementene av y fylles med "Heisann". De siste to elementene i x gir derimot TRUE i testen, og de tilsvarende elementene i y fylles med "Tittei". Kjør og verifiser.
Hvordan skal vi bruke dette? Testen vår blir nå str_detect(Protein, "\\*"). Vi sjekker om translatert sekvens inneholder stopp-kodon. Dette er enten TRUE eller FALSE. I de tilfellene der det er FALSE beholder vi bare Protein som den er, men hvis TRUE så trimmer vi av 1 base, og gjør en ny translasjon. Koden kan se noe slikt ut:
vipera.tbl %>%
mutate(Protein = translate(Sequence, trans.tab = 4)) %>%
mutate(Protein = if_else(str_detect(Protein, "\\*"),
translate(str_sub(Sequence, 2), trans.tab = 4),
Protein)) -> vipera.tblLegg dette inn i din kode, og legg på enda en mutate der du igjen tester om det finnes stopp-kodoner i Protein, og isåfall starter på base 3 før en ny translasjon.
Dette skal nå ha fanget opp alle muligheter for feil leseramme. Dersom sekvensene fortsatt inneholder stopp-kodoner, så må det skyldes at sekvensen ikke er fra en kodende sekvens i det hele tatt. Da må vi nesten se det som en falsk positiv fra BLAST-søket, og derfor filtrere vekk hele den raden fra tabellen.
Til slutt kan du jo lage en multippel sammenstilling av disse protein-sekvensene. Dette virker lurt, som en sjekk på at disse sekvensene virkelig fra samme kodende gen, hentet fra ulike slanger. De bør alle være ganske like over lengre partier.
Skriv vipera.tbl til en FASTA-fil med writeFasta(). NB! Kun en kolonne med navn Sequence blir skrevet ut av denne funksjonen, så kolonnen Protein må muteres inn i Sequence for å få skrevet ut proteinene.
Bruke MUSCLE hos EBI, slik vi gjorde i en oppgave over. Last opp FASTA-filen med slange-sekvensene. Hvis du ser litt over sammenstillingen så kan du muligens se om noen av sekvensene er falske positive, dvs at de skiller seg helt fra de andre. Vi skal neste uke bruke disse dataene videre.