Strålingsdata

Vi vil se på data fra et forsøk der man ville se om bestråling av hud gir endringer i genenes uttrykk i hudcellene. I forsøket målte man gen-uttrykket hos alle menneskelige gener, men vi ser bare på et utvalg på \(1000\) gener her nå.

I forsøket målte man gen-uttrykket i hudceller hos 10 personer, og for hver person både i ubehandlet vev ("Untreated") og i vev som er bestrålt ("Radiated"). Fra person \(1\) har vi en måling med navn "Radiated 1" og en med navn "Untreated 1", og tilsvarende for person \(2,3,...,10\). Alle målingene er normalisert mot en referanse, og log-transformert. Dette betyr at en positiv verdi her angir at genet er høyere uttrykt i forhold til referansen, og en negativ verdi at genet er lavere uttrykt i forhold til referansen. Alle tall er således sammenlignbare.

Tallene er samlet i en matrise X med 1000 rader (en for hvert gen) og 20 kolonner (en for hver prøve). La oss laste inn dataene direkte fra serveren

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

Vi kan ta en kort kikk på verdiene til de første \(4\) genene over de første \(5\) prøvene:

print(X[1:4,1:5])
##                 Radiated 1 Radiated 2 Radiated 3 Radiated 4 Radiated 5
## Gene_25k_809938       0.14      -0.04       0.35      -0.21       0.03
## Gene_206632          -2.60      -3.99      -1.46      -3.10      -3.60
## Gene_295923           0.57       0.51       0.12       0.52       0.58
## Gene_358736          -0.72      -0.06      -0.59      -0.38      -0.68

Oversikt over prøvene

Noe av det aller første man gjør når man sitter med en slik data-matrise er å lage noen figurer der vi visuelt kan danne oss et bilde av hvordan helheten ser ut, mer detaljerte studier kan komme senere. Det kan jo hende alt er bare rør, og at hele forsøket må forkastes…

I et forsøk som dette vil vi vel typisk se etter

Når vi vil inspisere prøvene (ikke genene) så må vi snu tabellen slik at radene er prøver og gener er kolonner. Dette kalles å transponere matrisen, og gjøres i R med funksjonen t(). Et PCA-plot vil alltid gi oss et bilde av radene, og vi må selv velge hva vil vil studere. Vi transponener matrisen, og tar en kikk på de første 4 kolonnene:

X.transp <- t(X)
print(X.transp[,1:4])
##              Gene_25k_809938 Gene_206632 Gene_295923 Gene_358736
## Radiated 1              0.14       -2.60        0.57       -0.72
## Radiated 2             -0.04       -3.99        0.51       -0.06
## Radiated 3              0.35       -1.46        0.12       -0.59
## Radiated 4             -0.21       -3.10        0.52       -0.38
## Radiated 5              0.03       -3.60        0.58       -0.68
## Radiated 6             -0.55       -3.69        0.45       -0.32
## Radiated 7             -0.09       -3.96        0.77       -0.42
## Radiated 8              0.23       -4.78        0.54       -0.58
## Radiated 9             -0.58       -2.69        0.54       -1.02
## Radiated 10             0.01       -3.35        0.39       -0.46
## Untreated 1             0.35       -2.81        0.47        0.17
## Untreated 2             0.38       -3.34        0.75        0.12
## Untreated 3             0.15       -3.00        0.69       -0.04
## Untreated 4             0.21        0.09        0.41        0.09
## Untreated 5            -0.07       -3.09        0.38       -0.31
## Untreated 6             0.13       -3.11        0.48        0.15
## Untreated 7             0.30       -3.99        0.47        0.13
## Untreated 8             0.97       -2.18        0.54       -0.08
## Untreated 9            -0.66       -2.30        0.34       -0.29
## Untreated 10            0.15       -4.82        0.46       -0.16

Hvis vi bare ser på tallene for de første \(4\) genene over, så virker det ikke som noen av prøvene er veldig spesielle. Første koordinat (genet "Gene_25k_809938") har verdier rundt \(0\) for alle prøver, andre koordinat er tydelig negativ for stort sett alle prøver (unntak for prøve "Untreated 4"), tredje er svakt positiv for alle prøver, etc. Men, vi må ta i betraktning alle \(1000\) genene, og kan ikke inspisere tall for tall på denne måten.

Vi plotter prøvene som punkter i et plan

Dersom vi hadde hatt data fra kun 2 gener så kunne vi ha plottet prøvene som punkter i et koordinat-system med 2 akser, en for hvert gen. La oss gjøre dette for de 2 første genene i matrisen.

Før plotting skyfler jeg data fra matrisen over i en tabell, og splitter rad-navnene i to ulike vektorer; en som angir behandling og en som angir person:

X.transp[,1:2] %>% 
  as_tibble(rownames = "Treatment") %>% 
  separate(Treatment, into = c("Treatment", "Person"), sep = " ") -> tbl

Da er vi klare for å plotte:

ggplot(tbl, aes(x = Gene_25k_809938, y = Gene_206632)) +
  geom_point(aes(color = Treatment), size = 5) +
  geom_text(aes(label = Person), nudge_x = 0.06)

Det var ikke så lett å se noe forskjell på bestrålte og ubestrålte prøver her, og dessuten er dette bare for 2 av genene. Vi har jo \(1000\) gener, og hvis vi skal plotte alle par mot hverandre slik må vi lage \((1000\cdot 999)/2\) ulike figurer! Det er en håpløs strategi…

Enda en håpløs figur

Vi kan forsøke å plotte hver prøve som en slags kurve, der vi på den horisontale aksen viser verdiene for gen 1,2,…,1000. Da får vi ihvertfall med alle data i samme figur. Dette blir seende slik ut:

matplot(t(X.transp), pch = 16, cex = 0.1, xlab = "Gen nummer", ylab = "Ekspresjons-verdi")
matlines(t(X.transp))

Dette er en håpløs figur, der de 20 kurvene legger seg mer eller mindre oppå hverandre, og alle hopper opp og ned hele tiden. Det er umulig å se om noen prøver (kurver) skiller seg fra de andre her.

PCA i R

La oss kjøre en PCA i R, med funksjonen prcomp():

pca <- prcomp(X.transp)

En PCA betyr at vi “ser” prøvene i et rotert koordinat-system, som er valgt slik at så mye variasjon som mulig er forklart langs de første retningene. Det første man derfor ser på er vanligvis hvor mye av den totale variasjonen i dataene vi har langs de ulike retningene:

relativ.forklart.varians <- pca$sdev^2/sum(pca$sdev^2)
barplot(relativ.forklart.varians[1:20], names.arg = 1:20, col = "tan4",
        xlab = "Prinsipal komponent", ylab = "Andel varians forklart")

Vi ser at første retning forklarer over 30% av all variansen, og at andre retning er nær 15%. Det betyr at disse to tilsammen viser oss nesten halvparten av all variansen i dataene. La oss derfor se hvordan prøvene legger seg i forhold til hverandre hvis vi ser på de gjennom disse to ‘vinduene’. MERK: Dette viser altså ikke den hele historien som data kan fortelle, men omfatter såpass mye av variasjonen at det trolig har verdi.

Først kopierer vi score-verdiene fra vår PCA inn i tabellen vi har fra før:

tbl <- tbl %>%
  mutate(PC1 = pca$x[,1]) %>% 
  mutate(PC2 = pca$x[,2])

Merk at i objektet pca finnes en matrise kalt x (elendig navn!) der kolonnene inneholder score-verdiene for prinsipal komponentene. Første kolonne er PC1, andre kolonne er PC2, osv.

Vi lager så samme plot som over, men bruker PC1 og PC2 som nye koordinater:

ggplot(tbl, aes(x = PC1, y = PC2)) +
  geom_point(aes(color = Treatment), size = 5) +
  geom_text(aes(label = Person), nudge_x = 0.7) +
  coord_fixed()

Legg merke til at vi bruker coord_fixed() her. Denne bør du alltid legge inn når du lager PCA-plot. Det sikrer at begge aksene har samme ‘bredde’, dvs at enhetene er like store. Dette gjør at avstandene vi ser mellom punktene blir korrekte. Dersom den ene aksen trekkes voldsomt ut kan vi bli lurt av det.

Vi ser med en gang to ganske klare effekter:

Dette betyr at det er en tydelig effekt av bestråling på det store bildet av gen-ekspresjonen hos hudceller, og vi kan nå lete etter disse effektene i større detalj, og finne hvilke gener dette angår. Vi ser også at Person \(2\) er et avvik fra resten, og vi bør vurdere om data fra denne personen skal tas ut allerede nå ettersom dette vil kunne tilsløre effektene vi ser hos resten.