For å få godkjent denne øvingen må du ha svart på alle oppgavene markert med rosa bakgrunn, dokumentert svarene (for eksempel med kode, utregning, kommando, forklaring og/eller screenshot av resultatet), ha riktig på 75% og levert riktig type fil i Canvas innen fristen.
Besvarelser leveres som RMarkdown-genererte html-filer i Canvas og dere leverer i grupper på to-og-to. Husk å skrive inn navnene deres i rapporten og i tillegg gruppe-nummeret. Fristen er søndag kl. 21.00 den gjeldende uka. Dere får bare en sjanse per uke siden fasit legges ut på mandag.
I BIN210 skal du levere rapporter laget med RMarkdown. Dette gjør du enkelt i RStudio. I slike rapporter kan du kombinere vanlig tekst med kode som gjør all håndtering av data, utfører beregninger og lager figurer eller tabeller.
Du genererer (“strikker”) en HTML-fil fra RMarkdown-filen og det er denne du leverer i Canvas. Vi vil at kode ikke skal skjules i disse rapportene. Bruk de samme overskriftene som på oppgavene.
Her er en kort film med litt introduksjon til RMarkdown
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 et sted inni 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. Spør AI om dette 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 igjen.
Gå til siden for Multiple Sequence Alignment hos EMBL-EBI
Velg Clustal Omega og lim inn sekvensene eller last opp FASTA-fila og kjør i vei med default parametre. Velg output på ClustalW format. Konserverte biter av sammenstillingen amrkeres med en stjerne under sekvensene.
Gjør deretter tilsvarende for Muscle.
Gjør deretter tilsvarende for Mafft.
Gjør deretter tilsvarende for T-coffee
Lag en liten tabell som viser resultatene. Første kolonne lister verktøyene du prøve og andre kolonne angir den konserverte del-sekvensen som ble funnet i sammenstillingen.
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 sammenstillinger 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 av disse ved å legge inn sekvensene i en bestemt rekkefølge.
Angi rekkefølgen på hvordan disse sekvensene sammenstilles ved å manuelt lage et guide-tre med UPGMA-metoden.
HINT: Se f.eks. Figure 8.3 på side 279 i boka til Zvelebil&Baum.
La oss se hvordan vi kan lage multiple sammenstillinger direkte i R.
Innstaller R pakken msa fra Bioconductor.
Med denne kan vi bruke ClustalW, Clustal Omega eller Muscle til å lage slike sammenstillinger. I en tidligere oppgave denne uka brukte vi sekvenser fra E. coli:
Her lager vi litt kode for å komme i gang:
library(msa)
ecoli.ss <- readDNAStringSet(filepath = "data/upstream30_Ecoli.fasta") # rediger sti
ecoli.msa <- msa(ecoli.ss, method = "Muscle")Legg merke til at
readDNAStringSet() istedenfor
med readFasta(). Dette fordi msa-pakken
forventer som input et objekt av typen XStringSet (fra
pakken Biostrings). Ved å lese inn data på denne måten får
vi et slikt objekt, og ikke som en tabell slik vi har gjort før.msa() og angir at vi vil
bruke Muscle som metode. Det objektet vi får tilbake er av typen
MsaDNAMultipleAlignment. Dette er å se på som en slags
liste, men vi kan ikke uten videre hente ut data fra dette
objektet. Dette er (dessverre) ganske vanlig for R pakker fra
Bioconductor.Vi kan se på sammenstillingen ved å bare skrive
ecoli.msa i Console-vinduet. Dette gir oss en rask, men
ikke så pen, oversikt over sammenstillingen:
ecoli.msa## MUSCLE 3.8.31
##
## Call:
## msa(ecoli.ss, method = "Muscle")
##
## MsaDNAMultipleAlignment with 8 rows and 37 columns
## aln names
## [1] ----AACGTCTGCTGGAATGG-CAGGAGGCCCATC-- GenomeSequence=NC...
## [2] ---TCGCGCTTCGTGCAATAA-AAGGAGGCGCAG--- GenomeSequence=NC...
## [3] -GTTATCGCCAGGCTT------TAGGAGGTTAATAAC GenomeSequence=NC...
## [4] -----TAGACTGATTTTCGGCTAAGGAGGAAGGCG-- GenomeSequence=NC...
## [5] GATTAAGGAAGGATTTTC----CAGGAGGAACAC--- GenomeSequence=NC...
## [6] -----ACGCTCCGGAAAATCGCGAGGAGGCCGTCG-- GenomeSequence=NC...
## [7] GATAATGGTTGCATGTACT---AAGGAGGTTGT---- GenomeSequence=NC...
## [8] AACCAGCGCTGCACT-------AAGGAGGGCCGCACA GenomeSequence=NC...
## Con ---?A?CGCT??ATTTA?T---AAGGAGG??CAC?-- ConsensusLegg merke til at nederst skrives ut en
konsensus-sekvens. Dette er en slags
‘flertalls-sekvens’ som viser hvilken base som er mest vanlig i hver
posisjon i sammenstillingen. Her ser det ut til at ? betyr
at flertallet er uavklart.
Men, vi vil selvsagt ha tak i selve sammenstillingen, og gjerne da
som en tabell med kolonnene Headerog Sequence
slik vi er vant til. Vi legger til følgende kodelinjer:
seq.lst <- ecoli.msa |>
msaConvert(type = "seqinr")
msa.tbl <- data.frame(Header = seq.lst$nam,
Sequence = seq.lst$seq)Nå skal msa.tbl være en tabell med Header
og Sequence som vi kan bruke videre.
I den siste koden over brukte vi først msaConvert() til
å konvertere fra det komplekse msa-objektet til en vanlig
liste, og så laget i en tabell ved å hente ut to elementer fra den
lista.
Kjør koden og se på msa.tbl i Console-vinduet.
Med R pakken DECIPHER kan vi enkelt visualisere en
multippel sammenstilling med funksjonen BrowseSeqs(). Den
lager resultatet som en HTML-fil som vi kan åpne i en nettleser eller
innkludere i et RMarkdown-dokument.
Installer pakken DECIPHER fra Bioconductor.
Funksjonen BrowseSeqs() krever at vi gir som input
sammenstillingen som et XStringSet objekt, slik vi hadde
for rådataene over. Vi kan lage en slik fra tabellen
msa.tbl ved å bruke funksjonen DNAStringSet()
og gi den sekvensene som input. Vi bruker deretter names()
til å bruke Header tekstene som navn:
library(DECIPHER)
msa.ss <- DNAStringSet(msa.tbl$Sequence)
names(msa.ss) <- msa.tbl$HeaderNå kan vi bruke funksjonen BrowseSeqs() fra
DECIPHER til å lage en HTML-fil som vi igjen innkluderer i
vårt dokument:
msa.ss |>
BrowseSeqs(openURL = F) |>
htmltools::includeHTML()
20
'''''''''|'''''''''|'''''''''|'''''''
GenomeSequence=NC_002695.1;Strand=1;Left=52183;Right=54045 ----AACGTCTGCTGGAATGG-CAGGAGGCCCATC-- 30
GenomeSequence=NC_002695.1;Strand=1;Left=447209;Right=447976 ---TCGCGCTTCGTGCAATAA-AAGGAGGCGCAG--- 30
GenomeSequence=NC_002695.1;Strand=1;Left=278411;Right=278884 -GTTATCGCCAGGCTT------TAGGAGGTTAATAAC 30
GenomeSequence=NC_002695.1;Strand=1;Left=148949;Right=149389 -----TAGACTGATTTTCGGCTAAGGAGGAAGGCG-- 30
GenomeSequence=NC_002695.1;Strand=1;Left=213022;Right=215163 GATTAAGGAAGGATTTTC----CAGGAGGAACAC--- 30
GenomeSequence=NC_002695.1;Strand=1;Left=536699;Right=538480 -----ACGCTCCGGAAAATCGCGAGGAGGCCGTCG-- 30
GenomeSequence=NC_002695.1;Strand=1;Left=302753;Right=302953 GATAATGGTTGCATGTACT---AAGGAGGTTGT---- 30
GenomeSequence=NC_002695.1;Strand=1;Left=319558;Right=321891 AACCAGCGCTGCACT-------AAGGAGGGCCGCACA 30
Consensus RRYHMDVGHHNSVBD+++++++NAGGAGGNNVDBVMM 30
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). Du laster tabellen inn i R med
load(<filnavnet>) der du erstatter
<filnavnet> med navnet (og stien) til fila 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()Les inn tabellen med sekvensene fra huggormene. I
microseq-pakken er det en funksjon som heter
translate(), og som gjør en translasjon fra DNA til
protein. Denne funksjonen tar som input en vektor av DNA-sekvenser, og
returnerer en tilsvarende vektor med protein-sekvenser. Vi kan bruke
denne funksjonen direkte på en kolonne med DNA-sekvenser i en
tabell.
Det genet vi her ser på forekommer i mitokondrisk DNA hos
disse slangene. DNA i mitokondirene oversettes litt annerledes enn annet
DNA, ved at kodonet TGA oversettes til W som
er aminosyren tryptofan, ikke til stopp som ellers. Vi kan angi dette
ved å angi trans.tab = 4 som et argument til
translate(). Dette er en av de mange translasjons-tabellene
som brukes i biologien.
Lag en ny kolonne i vipers.tbl som
heter Protein, og bruk mutate() og
microseq::translate() til å fylle denne med translaterte
sekvenser. Bruk translation table 4 som angitt over.
NB! Merk at både R pakkene microseq og
Biostrings (lastet inn med msa) har funksjoner
som heter translate! Derfor er det viktig at du angir
microseq::translate() for å bruke den fra
microseq-pakken. Hvis du bare skriver
translate() så vil den siste som ble lastet inn brukes, og
det kan være “feil” versjon av funksjonen. Dette forekommer ikke så
sjelden, og kan gi feilmeldinger som er vanskelige å forstå.
Se litt på Protein sekvensene du fikk. Trolig finner du
symbolet * i noen av dem. Dette er tegnet for stopp-kodon,
og betyr at translasjonen har støtt på et stopp-kodon i DNA-sekvensen.
Dette er ikke bra, og betyr sannsynligvis at leserammen er feil, for vi
skal ikke ha stopp-kodon midt inne i en kodende sekvens. Her er hvordan
du raskt sjekker hvilke av de translaterte sekvensene som har
stopp-kodon:
vipers.tbl %>% # antar tabellen heter vipers.tbl
filter(str_detect(Protein, "\\*")) |>
select(Species, Length) |>
knitr::kable()| Species | Length |
|---|---|
| Atheris squamigera | 655 |
| Bitis nasicornis | 655 |
| Cerastes gasperettii | 659 |
| Daboia russelii | 704 |
| Echis carinatus | 658 |
| Echis pyramidum | 640 |
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.
Med funksjonen str_sub() fra tidyverse kan
vi trimme vekk biter av tekster. 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:
x <- "AGATAGGTAACGAT"
microseq::translate(str_sub(x, start = 1, end = -1)) # har stopp-kodon
microseq::translate(str_sub(x, start = 2, end = -1)) # har også stopp
microseq::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 betingelse er
sann, og noe annet utføres hvis den er usann. Funksjonen
if_else() tar typisk 3 input:
TRUEeller
FALSE (logical).TRUE så tilordnes det som kommer i
input 2.FALSE, så tilordnes det som kommer i
input 3.Et 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? Betingelsen 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:
vipers.tbl <- vipers.tbl %>%
mutate(Protein = microseq::translate(Sequence, trans.tab = 4)) %>%
mutate(Protein = if_else(str_detect(Protein, "\\*"),
microseq::translate(str_sub(Sequence, 2), trans.tab = 4),
Protein))Legg dette inn i din kode, og legg så 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 lage en multippel sammenstilling av disse protein-sekvensene. Dette virker lurt, som en sjekk på at disse sekvensene virkelig er fra samme kodende gen, hentet fra ulike slanger. De bør alle være ganske like over lengre partier.
Bruk koden vi laget tidligere i dette dokumentet.
Lag en multippel sammenstilling at huggorm-proteinene. Bruk artsnavnet
som Header tekst for hver sekvens, og husk at du nå skal
bruke protein-sekvensen, ikke DNA-sekvensen. Det finnes en
AAStringSet() tilsvarende DNAStringSet(). Lag
en visualisering av dette med BrowseSeqs(), og rapporter om
det er noen konserverte områder.
Vi skal neste uke bruke disse dataene videre, så det er lurt å lagre
resultatene herfra med save(). Lagre tabellen med de
translaterte sekvensene (proteinene) samt den multiple
sammenstillingen.