1 Fylogenetiske trær på ymse vis

Vi vil nå re-konstruere fylogentiske trær på litt ulike måter. Poenget er å bli litt kjent med hvordan vi i praksis kan lage slike trær.

1.1 Data

Mitokondriene har eget DNA som brukes som markør for å se på slektskap langt tilbake i tid, ihvertfall i et humant perspektiv. Grunnen er at mitokondrielt DNA (mtDNA) er mer beskyttet og derfor ofte bedre bevart i gamle prøver (beinrester og lignende). mtDNA arves kun fra mor, og kan derfor brukes til å følge en ”mor-linje” bakover i tid.

Last ned filen primates_mtDNA.fasta, som inneholder sekvenser fra mtDNA for en rekke primater.

Lag en multippel sammenstilling av sekvensene, og bruk Muscle hos EBI slik vi gjorde i uke 7. Last ned ei FASTA-fil med de sammenstilte sekvensene.

Vi har sett at ulike programvarer gir litt ulike sammenstillinger. Slik sett burde vi kjøre de samme dataene gjennom flere, slik vi gjorde i uke 7. Men, siden fokuset nå er på hvordan vi lager trær, så bruker vi ikke tid på det.

1.2 Distansebasert tre i R

Dette blir repetisjon fra forrige uke. Fra sammenstillingen, lag et fylogenetisk tre med Neighbor Joining i R. Vi bruker pakken ape som sist uke.

Les inn sammenstillingen med read.dna() funksjonen, og lagre i et objekt som du kaller msa.ape. Beregn avstander fra dette med dist.dna() og velg Jukes-Cantor modellen model = "JC69". Lagre i et objekt du kaller d.jc.

Bruk så objektet d.jc til å lage et tre med bionj() funksjonen, og lagre i objektet tre.nj.

Til slutt kan du plotte dette med plot(). Vi skal nedenfor se andre måter å plotte trær på i R, men den vanlige plot() funksjonen duger nå. Den kan også varieres en del, prøv f.eks. plot(tre.nj, type = "unrooted", edge.width = 2).

OTU’ene kan, utfra sine navn, deles inn i fire grupper. Hvordan harmonerer kladene i treet med dette? Hva er våre (European human) nærmeste slektninger blant disse?

1.3 Maximum parsimony tre i R

De diskrete metodene for tre-konstruksjon bruker sammenstillingen direkte, uten noen avstands-tabell. De har også det til felles at de scorer en tre-topologi på en eller annen måte, og så er ideen å søke gjennom mange trær for å finne det som gir best score.

Et maximum parsimony tre betyr at vi scorer et tre med en parsimony score som angir antall mutasjoner som må til for å forklare sammenstillingen, gitt et tre. Maximum parsimony betyr ‘størst enkelhet’, og vi vil derfor ha en så liten score som mulig (færrest antall mutasjoner).

La oss se hvordan MP er implementert i R. Installer R-pakken phangorn. Fra oppgaven foran har vi allerede et tre (tre.nj), og la oss score dette:

library(phangorn)
msa.ph <- read.phyDat("C:/BIN210/data/primates_mtDNA_muscle.fasta", format = "fasta")
parsimony.score <- parsimony(tre.nj, msa.ph)

Legg merke til at vi leser sammenstillingen inn på nytt (endre filnavn til ditt eget), nå med en funksjon kalt read.phyDat(). Denne gir oss sammenstillingen på enda et format, forskjellig fra det vi fikk med read.dna() over. Deretter bruker vi funksjonen parsimony() til å beregne scoren for tre.nj gitt dataene i msa.ph. Kjør koden, hvor liten parsimont-score får du for dette treet og med denne sammenstillingen?

Med MP-metoden vil vi ha så lav parsimony score som mulig. Er det mulig å finne et alternativt tre som gir en lavere score? En funksjon for å optimere tre-topologien er optim.parsimony():

tre.mp <- optim.parsimony(tre.nj, msa.ph)
parsimony.score.opt <- parsimony(tre.mp, msa.ph)

Ble scoren mindre nå? Lag et plot av treet (slik du plottet tre.nj over). Blir treet annerledes enn NJ-treet?

En alternativ søke-metode er bab() som står for branch-and-bound. Denne er mer omstendelig, og tregere, men kan være verdt å prøve:

bab.obj <- bab(subset(msa.ph, 1:12))
parsimony.score.bab <- parsimony(bab.obj[[1]], msa.ph)

Blir scoren bedre? Lag et plot at treet, som ligger i bab.obj[[1]]. Er treet annerledes?

1.4 Maximum likelihood tre

Istedenfor a parsimony score bruker vi nå en likelihood score, som regel en log-likelihood. En likelihood er sannsynligheten for data (=sammenstillingen) gitt et tre. Mens vi med MP ville ha minst mulig score, vil vi nå ha størst mulig score (=log-likelihood).

Igjen kan vi bruke funksjoner fra phangorn til å score et gitt tre:

pml.obj.nj <- pml(tre.nj, data = msa.ph)
cat("Log-likelihood score for NJ-tree:", pml.obj.nj$logLik, "\n")

Score-verdien finner du i pml.obj.nj$logLik. Hva er verdien?

En log-likelihood er alltid negativ, og vi vil altså ha den så stor som mulig (minst mulig negativ). Igjen kan vi søke etter det optimale treet, denne gangen med optim.pml():

pml.obj <- optim.pml(pml.obj.nj, model = "JC", optNni = T)

Kjør koden. Hva blir scoren nå? Lag et plot av treet, som ligger i pml.obj$tree.


I eksemplene over ble resultatet av de diskrete metodene ikke så annerledes enn NJ-treet som vi starte ut med. Enten var NJ-treet vi starte ut med nær optimalt, eller så finnes en helt annen løsning som er den optimale, og at tre-søkingen ikke fant denne. Vi har hele tiden brukt den svært enkle Jukes-Cantor modellen, kan hende de diskrete metodene ville ha gitt større forbedringer med andre modeller.

1.5 ML-tre fra IQ-TREE

Siden ML-trær krever mye beregninger er det en rekke spesial-programvare for dette. La oss prøve IQ-TREE. Her er en kort oppskrift:

  • Under Input Data, les inn filen med MUSCLE-sammenstillingen som du brukte over.
  • Under Substitution Model Options, velg JC under DNA

La resten være som det er, og Submit Job.

Du tas så til en side for resultater. Den ser omtrent slik ut:

Hvis du klikker på linken øverst i det høyre vinduet så oppdaterer siden seg. I det venstre vinduet listes dine job’er, og om beregninger pågår (waiting, running) eller er ferdige (success) som i bildet over. Under det venstre vinduet skal det være en Download-knapp. Klikk på den, og du skal få lastet ned et zip-arkiv. Pakk det ut. Nede i mappen guest ligger ei ny mappe med et kryptisk navn. Nede i den finner du blant annet ei .pdf-fil med treet som en figur. Ble dette treet annerledes enn de andre du laget i R?

I samme arkiv er det også ei .treefile-fil med treet på Newick format. Denne kan du lese inn i R hvis du vil plotte treet der.

1.6 MEGA

MEGA er en gratis programvare for å lage trær som du kan installere på din egen PC. De som får problemer med installasjonen kan hoppe over denne oppgaven, men det kan hende du får glede av MEGA senere, f.eks, under oblig 2, så prøv. Dette er en typisk ‘trykk-på-knapp’ programvare, og det skal være rimelig enkelt å finne fram i dette på egenhånd. Vi lager derfor kun en lettvekts-forklaring på prosedyrene. Last ned fra http://www.megasoftware.net/.

1.6.1 Sammenstillinger

Bruk Data knappen øverst i MEGA-vinduet til å åpne fasta-filen vi brukte over (primates_mtDNA.fasta). Velg Align (ikke Analyze). Sekvensene åpnes i et under-vindu. På dette under-vinduet er det en Alignment-meny, velg Muscle, og bruk default parametre til å lage en multippel sammenstilling. Bruk så Data-menyen til å lagre (Export Alignment) i FASTA-format.

1.6.2 Avstander

Bruk Distance knappen øverst i MEGA-vinduet, og velg Compute Pairwise Distances, og åpne fasta-filen fra over. Beregn Jukes-Cantor distanser (Model/Method under Substitution Model). Sekvensene er ikke protein-kodende.

1.6.3 Trær i MEGA

Lag et tre med Neighbor Joining metoden basert på distansene. Bruk meny-knappen Phylogeny på MEGA-vinduet. Bruk 500 bootstrap-samples (Test of Phylogeny).

Lag så (minst) ett tre basert på Maximum Likelihood. Bruk samme data og modeller som før. Prøv med 500 bootstrap-samples her også, men reduser til 100 hvis det går fryktelig tregt.

Blir trærne laget med NJ og ML like de som du fikk i oppgavene foran?

2 Eksempel på R-spørsmål fra tidligere eksamener

Riktignok blir det fire svar-alternativer på eksamen, men her kan dere tenke litt på egenhånd. Fasit finner dere selvsagt ved bare å kopiere dette inn i R og kjøre det…

Hva slags innhold har a etter at følgende kode er kjørt:

b <- 1
for(i in 1:3){
  a <- b + i
}

Hva slags innhold har a etter at følgende kode er kjørt:

a <- 1
for(i in 1:3){
  a <- a + i
}

Hva slags innhold har a etter at følgende kode er kjørt:

b <- "ACTAGGCTA"
a <- str_replace(b, "A", "a")

Hva slags innhold har a etter at følgende kode er kjørt:

b <- "ACTG"
a <- str_c("A", b, b, "A")



3 R-lab

3.1 Plotte trær med ggtree

R pakken ggtree kan brukes for å visualisere trær i R. Du må altså først ha gjort alle beregninger, enten i R eller med en annen programvare, slik at du har et tre av den typen vi får når vi bruker ape pakken. Trær fra andre programvarer (f.eks. MEGA eller IQ-TREE) kan lagres på ei fil på Newick-formatet, og leses inn i R.

Pakken ggtree må innstalleres på en litt annerledes måte enn vi har sett før. Den ligger ikke på CRAN, men på det som kalles Bioconductor. Du finner den her. Alle Bioconductor-pakker installeres på samme måte, og du finner kode-snuttene på web-siden i linken. Kopier og kjør direkte i Console-vinduet.

For å plotte må vi ha et tre, og vi bruker bare noe fra oppgaven over:

library(ape)
msa <- read.dna("C:/BIN210/data/primates_mtDNA_muscle.fasta", format = "fasta")
d.jc <- dist.dna(msa, model = "JC69")
tre.nj <- bionj(d.jc)

Vi kan lage et plot på enkleste vis slik:

library(ggtree)
mitt.gg.tre <- ggtree(tre.nj, layout = "rect") +
  geom_tiplab(size = 4) +
  geom_treescale(width = 0.10) +
  ggplot2::xlim(0, 0.3)
print(mitt.gg.tre)

La oss gå gjennom koden. Vi oppretter figur-objektet med ggtree() på samme måte som vi lager en figur med ggplot(). Input er treet, og i tillegg angir vi layout som rektangulær. Vi kan lage trær med mange ulike utseender. Deretter legger vi på ymse geomer, akkurat som når vi bruker ggplot(). Her legger vi tekster på treets blad-noder (geom_tiplab()). Deretter vil geom_treescale() legge inn en strek som viser en avstand på 0.10 i treet, slik at vi har en peiling på hvor lange grenene er. Til slutt bruker vi funksjonen xlim() fra ggplot2pakken til å skalere x-aksen. Når figur-objektet er laget blir det plottet ved å print()e det.

Prøv å fjerne noen av geomene eller den siste skaleringen for å se hva som da skjer.

Det er svært mange muligheter for å lage fine trær med ggtree, men la oss nå se på en egenskap som åpner mange muligheter.

Ikke sjelden har vi en tabell med data omkring de sekvensene som treet bygger på. La oss lage en slik tabell for dette eksempelet:

library(microseq)
readFasta("C:/BIN210/data/primates_mtDNA_muscle.fasta") %>% 
  mutate(Type = "Human") %>% 
  mutate(Type = if_else(str_detect(Header, "Gorilla"), "Gorilla", Type)) %>% 
  mutate(Type = if_else(str_detect(Header, "Chimp"), "Chimpanzee", Type)) %>% 
  mutate(Type = if_else(str_detect(Header, "Orangutan"), "Orangutan", Type)) -> tbl

Vi fikk også her repetert dette med bruk av if_else() som vi har sett før! Vi laget her en ny kolonne Type som sier noe om hva slags dyr de ulike organismene er. Det er lett å tenke seg at vi ofte har en tabell med slike ekstra opplysninger omkring hver sekvens vi lager trær av.

Vi vil nå bruke kolonnen Type fra denne tabellen til å fargelegge bladene i treet. Da gjør vi slik:

mitt.gg.tree <- ggtree(tre.nj) %<+% tbl +
  geom_tiplab(aes(color = Type), hjust = -0.05) +
  geom_tippoint(aes(color = Type), size = 5) +
  geom_treescale(width = 0.10) +
  theme(legend.position = "right") +
  ggplot2::xlim(0, 0.35)
print(mitt.gg.tree)

Legg først merke til hvordan vi bruker operatoren %<+% for å ‘lime’ tabellen tbl til ggtree-objektet. For at dette skal fungere må første kolonnen i tabellen inneholde de samme tekstene som er labels i treet. Det kommer nærmest av seg selv så lenge vi laget tabellen ut fra samme FASTA-fil som treet. Når vi leser fila med readFasta() så blir jo Header første kolonne, og dette er de tekstene treet har som blad-labels.

Merk også at vi nå i geom_tiplab() angir en mapping med aes() eksakt slik vi gjør i ggplot(). Vi mapper farge til kolonnen Typeslik at hver unike tekst i denne kolonnen gir en farge på tekstene i treet. Vi legger også på et nytt geom, geom_tippoint() som lager en dot på hvert blad, og legger samme farging på disse. En legend som angir fargene dukker opp, og den plasseres til høyre med theme(), og så til slutt en skalering av x-aksen igjen. Dette siste er noe man typisk må prøve seg fram med.

Vi kan mappe flere egenskaper til ulike kolonner vi måtte ha i tabellen vi limte inn. Dette gjør at vi elegant kan få fram informasjon om sekvensene fra tabellen i treet. Prøv f. eks. å sette label = Type inne i aes() til geom_tiplab().

3.2 Mer om logicals

Ved siden av løkker er dette med logicals var noe som kan være litt trøbblete å fordøye. Se kort tilbake på Rlab uke 3 der vi så dette for første gang.

Vi har sett hvordan vi kan opprette et objekt, og gi den et innhold, ved å bruke en tilordning, som angis enten ved <- eller enkelt likhets-tegn =. Stort sett er det tall eller tekster vi håndterer i R. Et unntak er logicals. Et objekt som inneholder en logisk verdi kaller vi en logical. Det er kun to logiske verdier: TRUE eller FALSE. Vi kan opprette et slik objekt på samme vis som de andre:

g <- TRUE

…men det vanlige er at en logical opptrer som resultat av en sammenligning. Se på denne:

a <- 4
b <- (a < 5)

Hva skjer i siste tilordning? På høyre siden gjør vi en sammenligning, vi kommer på en måte med utsagnet: “a er mindre enn 5”. Dette utsagnet er enten sant (TRUE) eller usant (FALSE). Resultatet av dette tilordnes som verdien til objektet b:

print(b)
## [1] TRUE

Vi har flere operatorer som kan brukes til slike sammenligninger

  • a < b betyr “a er mindre enn b”. Tilsvarende betyr > “større enn”
  • a <= b betyr “a er mindre enn eller lik b”. Tilsvarende har vi også >=.
  • a == b betyr “a er lik b”. Husk dobble likhets-tegn, ellers blir det en tilordning!
  • a != b betyr “a er ikke lik b”.

Det er typisk tall vi sammenligner slik, men == og != brukes også ofte for å sammenligne tekster. Et eksempel:

ord <- c("Nei","Nej","No","Njet","Nein")
er.engelsk <- (ord == "No")

Her lages først en vektor med fem ulike varianter av et ord. Så sammenlignes vektoren med teksten "No". Da vil hvert element i vektoren sammenlignes med "No", og vi får en logisk vektor, som så lagres i objektet er.engelsk. Den får disse verdiene:

print(er.engelsk)
## [1] FALSE FALSE  TRUE FALSE FALSE

Siden sammenligningen innvolverte en vektor på fem elementer så blir resultatet av sammenligningen også en vektor på fem elementer. Men, de fem elementene er logiske: enten var sammenligningen sann eller usann for hvert element.

Vi kan bruke funksjonen which() til å finne hvilke(t) element i en logisk vektor som har verdien TRUE:

idx <- which(er.engelsk)
print(idx)
## [1] 3

La oss se hvordan vi kan bruke logicals i R til å raskt finne fram til enkelte elementer i en datastruktur.

3.3 Eksempel - bruk av logicals

Jeg har en fil ved navn samle_fil_16S_msa.fasta som inneholder noen 16S-sekvenser som har blitt sammenstilt med Muscle. Fila er på FASTA-format. Vi så forrige uke hvordan vi bruker funksjoner fra ape pakken for å lese inn og beregne p-distanser basert på en slik fil:

library(ape)
data <- read.dna("C:/BIN210/data/samle_fil_16S_msa.fasta", format = "fasta")
p.distanser <- dist.dna(data, model = "raw")

Objektet p.distanser er en avstands-tabell på et spesielt format. Her er et histogram over alle avstandene:

hist(p.distanser, col = "tan3", xlab = "p-distanse verdi",
     ylab = "Antall p-distanser", main = "")

3.3.1 Alle p-distanser fra samme art

Vi vil nå ikke ha alle avstandene, men kun de som er mellom sekvenser innenfor en og samme art. La oss finne avstandene innenfor arten tularensis.

For å få til dette er det først lurt å gjøre om avstands-tabellen til en matrise. Da er det nemlig enkelt å avlese rad- og kolonne-navn, og plukke ut bestemte rader/kolonner fra matrisa.

Hvis jeg har \(N\) sekvenser så skal jeg få en \(N\times N\) matrise som viser meg p-distansene mellom alle par av sekvenser. Dette får jeg til slik:

p.matrise <- as.matrix(p.distanser)

En slik konvertering fra en tabell til en matrise er stort sett alltid mulig. Både tabeller og matriser er jo ‘firkanter’ bestående av rader og kolonner. La oss ta en titt på en liten bit av denne matrisa:

print(p.matrise[16:17,16:17])    # skriver ut rad 16 til 17 og kolonne 16 til 17
##                      persica_9_kopi1 hispaniensis_1_kopi1
## persica_9_kopi1           0.00000000           0.01264045
## hispaniensis_1_kopi1      0.01264045           0.00000000

I rad 16 og kolonne 17 ligger nå p-distansen mellom sekvens 16 og 17. Generelt, i celle fra rad \(i\) og kolonne \(j\) ligger p-distansen mellom sekvens \(i\) og \(j\). Navnet til sekvensene er angitt i radene/kolonnene i matrisa. På alle diagonal-elementer skal det være p-distansen 0, fordi alle sekvenser er eksakt identisk like seg selv. I celle p.matrise[i,j] skal vi ha eksakt samme verdi som i p.matrise[j,i].

La oss ta en titt på alle rad-navnene til matrisen:

print(rownames(p.matrise))
##  [1] "noatunensis_2_kopi1"   "noatunensis_3_kopi1"   "noatunensis_4_kopi1"  
##  [4] "noatunensis_5_kopi1"   "noatunensis_6_kopi1"   "noatunensis_7_kopi1"  
##  [7] "noatunensis_8_kopi1"   "philomiragia_13_kopi1" "philomiragia_11_kopi1"
## [10] "philomiragia_12_kopi1" "philomiragia_17_kopi1" "philomiragia_15_kopi1"
## [13] "philomiragia_16_kopi1" "philomiragia_14_kopi1" "persica_10_kopi1"     
## [16] "persica_9_kopi1"       "hispaniensis_1_kopi1"  "tularensis_24_kopi1"  
## [19] "tularensis_25_kopi1"   "tularensis_18_kopi1"   "tularensis_20_kopi1"  
## [22] "tularensis_21_kopi1"   "tularensis_22_kopi1"   "tularensis_23_kopi1"  
## [25] "tularensis_26_kopi1"   "tularensis_27_kopi1"   "tularensis_19_kopi1"

Det er 27 rader (og 27 kolonner) i denne matrisen, og derfor 27 rad-navn (og kolonne-navn). Vi merker oss at første ‘ord’ i alle tekster synes å være arts-navnet, og at disse ‘ordene’ er adskilt med en "_".

Vi vil nå lage en ny avstands-matrise som kun inneholder avstandene mellom sekvenser fra arten tularensis. Vi må da ha tak i hvilke rader og kolonner som inneholder disse. Her er en måte å gå fram på, der vi bruker funksjonen word() som vi har sett før:

rad.navn <- rownames(p.matrise)
art <- word(rad.navn, 1, 1, sep = "_")

Vektoren art inneholder nå første ordet (før "_") i hvert rad-navn. Dette ordet skal være arts-navnet.

Legg merke til at vi her forteller word() at ord skal skilles med "_", altså ikke en blank, som er default. Merk også at bruk av word() bare fungerer fordi alle arts-navn alltid ligger på samme sted (først) i alle tekstene.

3.3.2 Vi bruker logiske vektorer

Nå gjør vi en sammenligning, får en logisk vektor, og bruker which() på denne:

er.tularensis <- (art == "tularensis")
idx <- which(er.tularensis)
print(idx)
##  [1] 18 19 20 21 22 23 24 25 26 27

Sammenligningen i første linje resulterer i en vektor av logicals, er.tularensis. Ta en titt på den, og verifiser at den er TRUE for de elementene der art er "tularensis", og FALSE ellers.

Funksjonen which() gir oss hvilke elementer som er TRUE, og objektet idx forteller oss derfor hvilke rader og kolonner som inneholder distanser mellom tularensis sekvensene. Vi kan nå bruke denne til å lage matrisen

p.matrise.tularensis <- p.matrise[idx,idx]

Til slutt kan vi konvertere denne matrisa tilbake til en avstands-tabell, og for eksempel lage et histogram over bare disse p-distansene:

p.distanser.tularensis <- as.dist(p.matrise.tularensis)
hist(p.distanser.tularensis, col = "darkgreen", main = "tularensis-avstander",
     xlab = "p-distanse verdi", ylab = "Antall p-distanser")

Vi kunne ha gjort dette med mindre kode, ved å bruke funksjonen str_detect() som vi har sett før. Hvordan ville du ha gjort det?