1 Newick-formatet

Forsikre deg om at du skjønner Newick-formatet (se forelesning). Tegn for hånd kladogrammet som er angitt her: (((OTU3,OTU4),OTU5),((OTU2,OTU6),OTU1))

Motsatt, skriv følgende kladogram i Newick-format:

Installer R-pakken ape, som inneholder mange funksjoner vi kan bruke for å gjøre fylogenetiske analyser i R. Her er et lite R-script som du kan kopiere, og som viser hvordan vi kan lese trær i newick-formatet inn i R, og plotte de:

library(ape)
tekst.med.newick.tre <- "tre1(((OTU3,OTU4),OTU5),((OTU2,OTU6),OTU1));"
tre <- read.tree(text = tekst.med.newick.tre)

Hvis newick-teksten over ligger på en fil med navn mitt_tre.txt, så leser vi det inn med

tre <- read.tree(file = "mitt_tre.txt")

Sleng så på en plot(tre) til slutt i scriptet, og du får fasiten på første oppgave over…

2 Evolusjonære avstander

Vi skal her se hvordan vi fra et utvalg sekvenser beregner evolusjonære avstander mellom de, og deretter grupperer sekvensene basert på disse avstandene. Vi er denne gangen opptatt av antallet grupper vi får.

Når man skal utforske et mikrobielt miljø med moderne sekvenserings-teknologi, så starter man ofte med det som kalles en ‘16S profiling’. Det betyr at vi forsøker å sekvensere alt vi finner av genet 16S rRNA. Dette genet finnes hos alle prokaryoter, og er den mest brukte markøren for å gjenkjenne bakterier og archaer. En ‘profiling’ betyr at vi forsøker å gruppere det vi finner i grupper som tilsvarer omtrent en art. Antallet slike grupper sier da noe om den biologiske diversiteten i det miljøet vi studerer, jo flere grupper jo høyere diversitet.

En vanlig tommelfinger-regel er å gruppere med en \(p\)-distanse på opp til \(0.03\) når vi studerer 16S-genet. Se ukens forelesning om \(p\)-distanser. Alle sekvenser som er mindre enn denne \(p\)-distansen fra hverandre skal tilhøre samme gruppe. Et 16S gen er omtrent 1500 baser langt. Omtrent hvor mange baser kan da være ulike mellom to slike gener fra samme gruppe?

Vi ser på et bittelite utvalg på 100 reads (sekvenser) her, fra magesekken hos laks, og altså finne omtrent hvor mange arter er det i laksens mage. Hver read er sekvensen til en kort bit av hele 16S genet, men alle reads kommer fra samme område av genet. Grovt sett kan vi skissere prosedyren slik:

Lag multippel sammenstilling av sekvensene
Beregn avstander mellom alle par av sekvenser
Lag et dendrogram med hierarkisk clustring og complete linkage (mer om dette nedenfor)
Kutt dendrogram-treet i klader ved avstand \(0.03\), slik at alle grupper får dette som største interne avstand

Vi gjør dette i R, og vil bruke funksjoner fra pakkene microseq og ape. Begge disse må være installert.

2.1 Last ned data

Vi har laget en multippel sammenstilling ved å bruke TCoffee hos EMBL-EBI, slik vi gjorde tidligere i BIN210. Vi kjørte også de samme sekvensene gjennom Muscle og Clustal på samme webside. Resultatene foreligger på disse FASTA-filene:

Last ned disse, og husk hvor du lagrer de på din PC.

2.2 Les inn data i R

De filene du lastet ned er FASTA-filer, og vi kan som vanlig lese slike filer inn i med readFasta() fra microseq pakken. Dette gir oss da en tabell med kolonnene Header og Sequence:

library(tidyverse)
library(microseq)
library(ape)
msa.tbl <- readFasta("C:/BIN210/data/tcoffee_OTU_msa.fasta")  # endre sti hos deg

Vi vil nå bruke funksjoner i R-pakken ape til å beregne evolusjonære avstander. Disse vil ikke ha dataene som en tabell! Vi må derfor gjøre om dette til et format som kalles DNAbin:

msa.tbl %>%               # fra tabell...
  msa2mat() %>%           # ...til matrise...
  as.DNAbin() -> msa.bin  # ...til DNAbin

Det skal sies at dersom du leser inn fila med funksjonen read.dna() fra ape-pakken, så hadde du fått dataene på dette formatet direkte, uten å gå veien om tabell.

2.3 Beregne evolusjonære avstander

Vi beregner nå enkelt de evolusjonære avstandene basert på innholdet i msa.bin med funksjonen dist.dna() som også finnes i pakken ape:

avstander <- dist.dna(msa.bin, model = "raw")

Her må vi angi hva slags evolusjonær modell vi vil bruke. Ved å angi model = "raw" så får vi \(p\)-distanser. Ved å angi model = "JC69" så får vi Jukes-Cantor avstander. Det finnes mange flere å velge blant, sjekk Help-fila for denne funksjonen. Vi klarer oss fint med \(p\)-distanser her.

Objektet avstander inneholder nå alle parvise avstander mellom de 100 sekvensene. Det blir \(100\cdot 99/2\) som er 4950 verdier. Du kan lage et enkelt histogram over de for å se hvor mange som er store og små:

hist(avstander, col = "slategray3")

2.4 Gruppere sekvenser

Vi vil gruppere sekvensene med en metode som heter hierarkisk clustring. Dette er en variant av UPGMA-metoden (se boka og forelesning). Grovt sett så betyr dette å slå sammen sekvenser som til enhver tid er de ‘mest like’ og på denne måten bygge opp en slags tre-struktur som kalles et dendrogram. En linkage funksjon angir hva vi mener med ‘mest like’. Vi vil bruke det som kalles complete linkage denne gangen. Vi skal snakke mer om dette mot slutten av BIN210. Her er koden for å lage en slik tre-struktur:

dendrogram <- hclust(avstander, method = "complete")

Funksjonen hclust() er en av basis-funksjonene i R.

Ofte vil man plotte et slikt tre, men det hjelper oss ikke så mye denne gangen. Prøv plot(dendrogram). Du kan fjerne sekvens-navnene på denne måten: plot(dendrogram, labels = FALSE).

Vi ønsker heller å kutte treet opp i klader slik at avstanden mellom alle blader i en og samme klade er mindre enn \(0.03\). Hver slik klade tilsvarer da en gruppe ettersom alle sekvenser i samme klade ligger med avstand mindre enn \(0.03\) fra hverandre.

klader <- cutree(dendrogram, h = 0.03)
cat("Fant", length(unique(klader)), "grupper\n")

Sett dette sammen til ett script, og se hvor mange grupper du finner i disse dataene når vi bruker \(p\)-distanse grensen på \(0.03\) (h = 0.03). Prøv å endre grensen til \(0.02\) og \(0.04\), og se om det blir en (stor) endring i antall grupper.

Prøv til slutt å gjenta prosedyren, men med de andre multiple sammenstillingene (muscle og clustal) som input. Får vi andre resultater basert på de andre sammenstillingene? Det er de samme sekvensene vi starter ut med, så alle forskjellene du eventuelt får må derfor skyldes måten de ble sammenstilt på av de tre metodene.

3 Eksamens-spørsmål

En smakebit fra tidligere eksamens-spørsmål. Årets eksamen blir av samme type, dvs flervalgs-oppgaver.

3.1 Spørsmål 1

Hva er ikke riktig om evolusjonær avstand og observert avstand (p-distanse)?

Evolusjonær avstand måles i mega-år.
Forskjellen mellom p-distanse og evolusjonær avstand øker med økende p-distanse.
Evolusjonær avstand er lik eller større enn p-distanse.
p-distansen er alltid et tall mellom 0 og 1.

3.2 Spørsmål 2

Seks sekvenser, A, B, C, D, E og F, skal sammenstilles. Guide-treet er angitt som ((D,(C,((B,F),E))),A).

Hvilke to sekvenser er de første til å sammenstilles?

D og A
B og F
C og E
F og E

3.3 Spørsmål 3

Hvilken påstand er ikke riktig?

Ingen paraloger er orthologer.
Alle orthologer er homologer.
Alle paraloger er homologer.
Alle homologer er orthologer.

4 R-lab

4.1 Mer om løkker

Denne uka vil vi se litt mer på løkker. Dette er noe som faller litt vanskelig for mange, og samtidig er dette noe av det aller viktigste å skjønne og beherske. Det ligger en enorm kraft i dette at vi, på noen få linjer kode, kan kverne gjennom enorme jobber som det ville tatt evigheter å gjøre manuelt. Vi så bruk av løkker allerede i uke 5, når vi simulerte mange sekvenser (U’er). Da brukte vi en løkke for å gjenta kode mange ganger.

4.2 Samle data fra mange filer

Ganske ofte får vi bruk for å lese inn innholdet av mange filer, og samle dette opp i en datastruktur i R, for eksempel i en tabell. La oss se på et slikt eksempel.

I en mappen ligger det et antall FASTA-filer. Hver FASTA-fil inneholder noen gen-sekvenser (16S-gener) fra noen genomer. Vi ønsker nå å

Lese alle filene inn i R, en etter en
Samle opp alle sekvensene i en stor tabell

For å kjøre eksemplene under laster ned disse filene (høyre-klikk, lagre som) til en mappe på din PC:

I virkeligheten har vi ofte mange flere filer, men dette holder for illustrere problemet.

4.2.1 Få tak i alle filnavn

Hvis vi har mange filer kan vi ikke manuelt skrive inn alle filnavn. Her er måten vi avleser navnet på alle filer og under-mapper som ligger i en spesiell mappe:

filnavn <- list.files(path = "C:/BIN210/data",
                      pattern = "_16S.fasta",
                      full.names = T) # endre til din sti

Funksjonen list.files() gir oss navnet på alt innholdet i den mappen vi angir navnet til i path. Argumentet pattern kan vi bruke til å avgrense med, her avleses nå kun filer som inneholder teksten "_16S.fasta". Hvis du vil avlese alle filer i mappen så angir du ikke noe pattern. Til slutt setter vi full.names = T, som angir at vi vil ha med hele stien i alle filnavn.

Sjekk at objektet filnavn blir en vektor med 3 tekster, og inneholder de riktige filnavnene, innkludert sti.

4.2.2 Lese filene og lagre innholdet

Vi utvider nå scriptet med ei løkke der vi leser inn filene etter tur:

library(microseq)   # trenger readFasta nedenfor
for(i in 1:length(filnavn)){
  cat("Leser fila:", filnavn[i], "\n")
  sekvens.tbl <- readFasta(filnavn[i])
}

## Leser fila: C:/BIN210/data/noatunensis_2_16S.fasta 
## Leser fila: C:/BIN210/data/persica_9_16S.fasta 
## Leser fila: C:/BIN210/data/tularensis_18_16S.fasta

Når vi leser inn en FASTA-fil med readFasta(), så får vi en tabell 2 kolonner (Headerog Sequence). Dette lagres her i objektet sekvens.tbl. Problemet er bare at hver gang løkka går, så leser vi inn sekvensene fra en ny fil, men disse overskriver det forrige innholdet i sekvens.tbl. Objektet fylles med nytt innhold, og vi sletter dermed det som lå der fra før. Slik sett vil resultatet av dette scriptet nå bli at sekvens.tbl bare inneholder resultatene fra den siste fila vi leste inn. Hvordan får vi tatt vare på alle de tidligere innleste sekvensene?

4.2.3 Legge til rader i en tabell

Dersom vi har to tabeller, og disse har de samme kolonnene, så kan vi skjøte de sammen med bind_rows(). Koden bind_rows(tbl1, tbl2) vil da legge tbl1 først, og deretter tbl2 under denne, og gi oss resultatet som en stor tabell. Tabellene stables oppå hverandre. Tabellene må ha like kolonne-navn for at dette skal bli vellykket.

I vår kode så leser vi inn en og en tabell med data. Da vil vi skjøte den sist innleste på de foregående. Men, når vi leser den første fila så har vi ingen ‘foregående’.

Vi oppretter derfor en tom tabell før løkka starter, og det er denne vi skjøter på i hver runde i løkka. Hele scriptet blir da noe slikt:

library(microseq)

filnavn <- list.files(path = "C:/BIN210/data", pattern = "_16S.fasta", full.names = T)
sekvens.tbl <- NULL
for(i in 1:length(filnavn)){
  sekvens.tbl <- bind_rows(sekvens.tbl, readFasta(filnavn[i]))
}

Når vi oppretter objektet sekvens.tbl og fyller det med NULL så betyr det bare at navnet registreres, men objektet har ennå ikke noe innhold. I første runde i løkka leser vi inn en tabell, og skjøter dette på dette NULL-objektet. Resultatet av dette er da kun den innleste tabellen. I alle påfølgende runder i løkka så legges siste innleste tabell under de foregående. På denne måten får vi nå tatt vare på alle innleste sekvenser. Tabellen sekvens.tbl vokser nå for hver runde i løkka.

4.3 Vi leker med løkker

La oss se nærmere på for-løkker slik at vi kan forstå litt bedre. En ting er å kopiere og endre litt på halv-ferdig kode, slik som over, men for å kunne bruke dette som et verktøy så må man nesten ha det litt mer under huden.

4.3.1 Definisjoner

Ei løkke skrives i R alltid på formen:

for(loop.objekt in vektor){
  
}

De fargede ordene (for og in) er nøkkelord i R-språket, og har bare denne ene spesielle betydningen. Du kan ikke bruke disse ordene som objekt-navn. Både parentesene og krøll-parentesene inngår også alltid i denne syntaksen, og du må skrive det eksakt slik hver gang. Det som kan byttes ut er loop.objekt og vektor.

Inni parentesen angis det alltid først et objekt, og over kalt loop.objekt. Vi har stort sett kalt denne i, men du velger selv hva du kaller denne. La oss benevne denne som løkke-objektet heretter.

Etter nøkkel-ordet in kommer en vektor. Vi har som regel skrevet 1:5 eller liknende. Dette er også en vektor, men den opprettes der og da. Vi kunne godt ha brukt en vektor som allerede finnes, vi skal se på det nedenfor. La oss benevne denne løkke-vektoren heretter.

Det som kommer mellom { og } er den koden som utføres hver gang løkka går. Her står det nå intet, og da er jo hele greia bortkastet. Vi lager en løkke utelukkende fordi den koden som står inne i løkka skal utføre gjentatte ganger. La oss benevne denne løkke-koden heretter.

4.3.2 Antall runder - løkke-vektoren

Følgende sannhet gjelder alltid:

Ei løkke går like mange runder som det er elementer i løkke-vektoren.

Kopier og kjør kode-snuttene under, og prøv å gjennomskue hva som skjer før du kjører:

for(v in 1:5){
  cat("Hei\n")
}

vek <- rep(2.5, times = 7)
for(OlaConny in vek){
  cat("Hei\n")
}

w <- c("The","power","of","looping")
for(i in w){
  cat("Hei\n")
}

Legg merke til at vi kan enten opprette en vektor der og da, inne i parenetesen (f.eks. 1:5) eller vi kan benytte en vektor som allerede finnes fra før. Uansett, antall elementer i løkke-vektoren bestemmer antall runder løkka gjør.

Sjekk innholdet i vektorene vek og w (bare skriv objektene i Console-vinduet, og return). Sjekk også hva slags verdier v, OlaConny og i har etter at løkkene er ferdige.

4.3.3 Løkke-objektet

I kode-snuttene over kalte vi løkke-objektet litt ymse (v, OlaConny, i). Vi velger selv dette navnet.

Vi brukte aldri denne løkke-objektet til noen verdens ting inne i løkke-koden i snuttene over. Det er lov. Her bestod løkke-koden bare av den samme utskriften hver gang, og den ble gjentatt for hver runde løkka gjorde. Men, som oftest bruker vi løkke-objektet til noe inne i løkke-koden. Derfor er det litt viktig å skjønne hva slags verdier denne får. Følgende sannhet gjelder alltid:

Løkke-objektet får verdiene til løkke-vektoren etter tur

La oss endre litt på kode-snuttene over. Kopier disse og kjør, igjen prøver du å tenke etter hva som vil skje før du kjører:

for(v in 1:5){
  cat("Løkke-objektet inneholder nå:", v, "\n")
}

vek <- rep(2.5, times = 7)
for(OlaConny in vek){
  cat("Løkke-objektet inneholder nå:", OlaConny, "\n")
}

w <- c("The","power","of","looping")
for(i in w){
  cat("Løkke-objektet inneholder nå:", i, "\n")
}

Løkke-vektoren fyller altså to funksjoner:

Dens lengde bestemmer antall runder i løkka
Dens innhold bestemmer innholdet til løkke-objektet i hver runde

4.3.4 Bruk av løkker

Vi så over at vi brukte en løkke til å gjennomløpe en datastruktur (en vektor med filnavn) og så utføre noe kode knyttet til hvert element i datastrukturen. I slike tilfeller er det naturlig at løkke-vektoren inneholder hel-tallene fra 1 opp til antall elementer i datastrukturen. På den måten kan vi bruke løkke-objektet som indeks for å angi hvilket element i datastrukturen vi vil gjøre noe med inne i løkke-koden.

Når vi skal gjennomløpe en datastruktur, så vet vi alltid på forhånd eksakt hvor mange ganger løkka skal gå. Dersom vi, i en eller sammenheng, har bruk for å gjenta kode, men ikke vet hvor mange ganger, så bruker vi ikke for-løkker. Det finnes andre løsninger for slike tilfeller, men det er ikke tema denne gangen.

4.3.5 Oppgave

I oppgave 2 over så leste vi inn en multippel sammenstilling, og fant antall grupper. Dette gjorde vi for 3 ulike sammenstillings-metoder, der sammenstillingene lå på 3 ulike filer. Prøv å utvide scriptet slik at du bruker ei løkke for å gjøre jobben for hver av sammenstillingene.

Ekstra, for de som vil strekke seg litt: For hver av sammenstillingene, beregn evolusjonære avstander med modellene "raw", "JC69", "K80", "F81"og "TN93", se Help-filen for dist.dna(). Grupper deretter som før (kutt ved h=0.03). NB! Bruk pairwise.deletion = TRUE i dist.dna(), ellers vil du få problemer i noen av tilfellene. Hva er effekten av denne opsjonen? Ta vare på antall grupper du får med hver av modellene, for hver av sammenstillingene. Hva påvirker antall grupper mest, valg av sammenstillings-metode eller valg av evolusjonær modell?

BIN210 - øvinger uke 8