1 Strålingsdata

Vi skal se på et eksempel med data fra en transcriptomics forsøk. Du kan se HTML-dokumentet her for å få en oversikt over data og også bruken av PCA.

La oss bruke dataene til å gjøre et par cluster-analyser. Vi gjør dette i R. Start RStudio, og opprett et R-script. Først laster vi inn data:

load(url("http://arken.nmbu.no/~larssn/teach/bin210/ovinger/data/radiation_small.RData"))

Matrisen som lastes inn heter X. Sjekk dimensjonene

print(dim(X))

og verifiser at det er \(1000\) rader (gener) og \(20\) kolonner (prøver).

2 Hierarkisk clustring av prøvene

La oss først bruke hierarkisk clustring, og se om de \(20\) prøvene grupperer seg på en eller annen måte. Vi må først transponere (snu) matrisen slik at prøvene ligger i radene (og genene i kolonnene).

X.transponert <- t(X)

Funksjonen t() transponerer en matrise.

Det første vi må gjøre ved enhver clustring er å beregne avstander mellom alle objekter (rader). La oss velge euklidsk avstander (se forelesning). Dette beregnes enkelt slik:

d.euclid <- dist(X.transponert)

Funksjonen dist() kan beregne mange typer avstander, men euklidsk avstand er default valget. Siden vi har \(20\) prøver skal det nå bli \(20\cdot 19/2\) parvise avstander, dvs 190. Lag gjerne et histogram over disse, for å se omtrent hvor store avstander vi får.

Vi bruker så hierarkisk clustring til å lage et dendrogram som viser hvordan de ulike prøvene ‘henger sammen’. Vi må da velge hva slags linkage funksjon vi skal bruke (se forelesning). Her brukes "average" linkage:

dendrogram <- hclust(d.euclid, method = "average")
plot(dendrogram, hang = -1, ylab = "Euklidsk avstand", main = "Average linkage")

Kjør dette, og se om dendrogrammet gir et resultat som henger sammen med det vi så ved PCA på de samme prøvene.

Prøv å endre linkage funksjonen til "complete" og deretter til "single", og se hvordan dette endrer bildet. Det er "average" linkage som er det vanlige, men "complete" linkage er også veldig aktuelt å bruke i en del tilfeller. Hva vil vi kunne oppnå med å endre fra "average" til "complete" linkage?

Prøv også Ward’s metode ("ward.D"). Hva er dette?

3 Partisjonering av genene

Vi kan bruke hclust() til å clustre gener også, men et dendrogram med \(1000\) objekter blir ikke veldig pent å se på. La oss derfor prøve noe annet for å gruppere genene.

Med \(k\)-means clustring, eller partisjonering, deles alle objekter i \(k\) grupper (se forelesning). Det finnes en \(k\)-means funksjon i R, men den bruker utelukkende euklidsk avstand. Vi vil også prøve å gruppere gener basert på korrelasjons-avstand, og et alterntiv er å bruke funksjonen pam() (Partitioning Around Medoids) fra cluster-pakken til dette. Installer først cluster-pakken.

Beregn euklidsk avstand mellom alle gener, og lagre i objektet d.euc. Hvor mange parvise avstander blir dette? Hva er en euklidsk avstand egentlig? Hva skal til for at to gener har en euklidsk avstand på (nesten) 0?

Beregn deretter korrelasjons-avstand mellom de samme genene. Dette er ikke innebygget i dist, så vi gjør det slik:

d.corr <- as.dist((1-cor(t(X)))/2)

Hva er en korrelasjons-avstand? Hva skal til for at to gener har en korrelasjons-avstand på (nesten) 0?

3.1 Partisjonering i \(k=2\) grupper

La oss dele genene inn i 2 grupper basert på euklidsk avstand:

library(cluster)
K <- 2
clst <- pam(d.euc, k = K)

Fra objektet clst kan vi nå trekke ut ymse informasjon. Blant annet medoidene i clusterne, dvs. de genene som ligger nærmest sentrum i sitt respektive cluster. La oss plotte disse med verdiene for de 20 koordinatene fortløpende:

library(tidyverse)
t(X[clst$id.med,]) %>% 
  as_tibble(rownames = "Sample") %>% 
  pivot_longer(names_to = "Cluster", values_to = "Expression", -Sample) -> tbl
fig <- ggplot(tbl, aes(x = Sample, y = Expression, color = Cluster)) +
  geom_point() +
  geom_line(aes(group = Cluster)) +
  theme(axis.text.x = element_text(angle = 90))
print(fig)

Lag figuren. De to kurvene angir ekspressjons-verdiene til de ‘typiske’ genene i hvert av de \(2\) clusterne. Basert på denne figuren, hva vil du si at disse to clusterne representerer?

Se litt på koden, og prøv å skjønne hva som gjøres i hvert steg.

3.2 Partisjonering med økende \(k\)

Prøv med økende verdier for \(k\), ihvertfall opp til 6. Basert på medoidene, hva vil du si er det karakteristiske for genene i hver gruppe?

Bytt til korrelasjons-avstander, og prøv verdier for \(k\) fra 2 til 6. Igjen, hva karakteriserer genene i hver gruppe?

I objektet clst ligger mye informasjon. I clst$clustering ligger en vektor som angir hvilken gruppe (cluster) hvert gen tilhører. Hvor mange gener havner i de ulike gruppene? (hint: table())

Vi så tidligere at person \(2\) hadde data som avviker fra resten. Kan du finne en gen-gruppe som viser nettopp dette at person \(2\) avviker fra de andre? Det kunne være interessant å finne hvilke gener som tilhører en slik gruppe, for da har person \(2\) et avvikende gen-uttrykk for disse genene. Kanskje er person \(2\) syk? Dette er noe man muligens burde følge opp videre herfra…

4 R-lab

4.1 Tabeller

I forbindelse med transkriptom-data er det typisk store tabeller med ekspressjons-verdier vi sitter med. La oss repetere en del ting vi har sett i forbindelse med slike tabeller.

Typisk har vi en tabell med transkriptom-data på en tekst-fil, der hver kolonne er adskilt med tabulator eller et annet tegn. Vi har tidligere lest en tabell inn slike filer i R med funksjonen read_delim(). Jeg har lagt en fil med transkriptom-data ut på nettet, og vi kan lese den direkte inn i R slik:

library(tidyverse)
adr <- url("http://arken.nmbu.no/~larssn/teach/bin210/Rlab/data/radiation_small.txt")
transkr.tbl <- read_delim(adr, delim = "\t")

Legg merke til følgende

  • Vi oppretter forbindelsen til en web-addresse med funksjonen url(). Dette blir det vi kaller en connection. Tenk på det som en forbindelse eller ‘rør’ ut til addressen. Vi kan deretter lese innholdet på websiden gjennom dette ‘røret’.
  • Når vi bruker read_delim() må vi vite hva slags symbol som skiller kolonnene i fila vi skal lese inn. Her er det tabulator, og da angir vi dette ved delim = "\t". Hvis det hadde vært semikolon eller komma (.csv-filer) hadde vi skrevet delim = ";" eller delim = ",".

4.2 Skal vi konvertere til matrise?

Når vi leser inn med read_delim() så blir resultatet (her kalt transkr.tbl) alltid en tabell, der hver kolonne er en vektor. De ulike kolonnene kan inneholde ulike typer data (tekst eller tall).

Tabellen vi leste inn over inneholder først en kolonne med tekst (gen-navnene) og deretter bare tall. Da er det i mange sammenhenger lurt å lagre kolonnene med tall i en matrise istedenfor i en tabell, og legge gen-navnene inn som rownames() til matrisen. I en matrise er alle kolonner og alle rader vektorer. Det betyr at alle celler i en matrise må inneholde data av samme type. Her er jo alt tall (bortsett fra første kolonne), så det er mulig å endre dette til en matrise.

Hva er poenget? Vel, hvis du henter ut rad nummer 1 fra en tabell, så blir dette en ny tabell, ikke en vektor. Hvis du henter ut rad 1 fra en matrise så er det en vektor. Det er en god del operasjoner vi kan gjøre på vektorer, men ikke på tabeller. Du kan heller ikke kjøre PCA på en tabell, men på en matrise. Ved å gjøre om fra tabell til matrise så har vi flere muligheter senere. I det hele tatt så bør vi tenke at en tabell er fin for å lagre data, mens en matrise egner seg bedre når vi skal utføre operasjoner på dataene.

Så lenge vår tabell består bare av tall, er det super-enkelt å endre dette til en matrise:

transkr.tbl %>% 
  select(-1) %>%               # velger bort kolonne 1
  as.matrix() -> transkr.mat   # gjør om til matrise
rownames(transkr.mat) <- transkr.tbl$Gene  # bruker tekstene fra kolonne 1 som rad-navn

Husk å fjerne kolonnen med tekster før du konverterer til matrise (slik vi gjorde med select(-1) over), ellers blir det en matrise med bare tekster!

Det er også fint å kunne konvertere en matrise tilbake til en tabell. Dette gjør vi da typisk slik:

transkr.mat %>% 
  as_tibble(rownames = "Gener") -> ny.tbl

Argumentet rownames = "Gener" betyr at matrisas rad-navn skal bli en ny kolonne i tabellen, og denne skal hete Gener.

4.3 Noen oppgaver

Prøv litt på følgende (enkle?) oppgaver.

  • Hva er gjennomsnittlige ekspressjons-verdi for alle gener og alle prøver? Finn også median, maksimum og minimum verdi.

  • Finn gjennomsnittlig ekspressjons-verdi for Radiated og Untreated for hvert gen, og beregn differansen mellom de. Lagre dette i en vektor (hint: løkke).

  • Sorter ekspressjons-matrisa slik at de genene med størst differanse fra oppgaven over kommer først, og deretter i gradvis minkende rekkefølge.

5 Noen potensielle eksamensoppgaver

Vel, kanskje ikke siden det er hjemme-eksamen dette året. Men, test deg selv, fasiten får du ved å kjøre koden i R.s

5.1 Oppgave 1

txt <- "Dette er en felle"
x <- length(str_split(txt, pattern = " "))

Hva slags innhold får x?

  1. "Dette" "er" "en" "felle"
  2. 1
  3. 4
  4. 17

5.2 Oppgave 2

data <- readFasta("sekvensdata.fasta")
x <- length(data$Header) == length(data$Sequence)

Hva slags innhold får x?

  1. TRUE
  2. 2
  3. FALSE
  4. 10

5.3 Oppgave 3

data <- readFasta("sekvensdata.fasta")
x <- length(data$Header == data$Sequence)

Hva slags innhold får x?

  1. Ett enkelt tall som angir antall sekvenser
  2. En vektor av tall, der elementene angir lengden til hver av sekvensene
  3. En vektor av logiske (TRUE eller FALSE) verdier
  4. Bare TRUEeller FALSE avhengig av sekvensene som ble lest inn

5.4 Oppgave 4

v <- c("Dette","er","ikke","noen","felle")
x <- str_length(v)

Hva slags innhold får x?

  1. 5
  2. 1
  3. 5 2 4 4 5
  4. "Dette er ikke noen felle"