Kapitel 11 Emner fra eksperimental design
“Amatører sidder og venter på inspiration, resten af os står bare op og går på arbejde.” - Stephen King
Figure 1.1: Image just for fun
11.1 Inledning og læringsmålene
Formålet med dette kapitel er at spørge: Hvordan kan vi anvende de værktøjer, vi har lært i kurset, til at undersøge forskellige emner inden for eksperimentelt design? Dette er ikke en grundig introduktion til eksperimentelt design, men snarere en præsentation af nogle nyttige og interessante emner, der effektivt illustrerer, hvorfor det er vigtigt at lave passende visualiseringer af dine data. Forståelsen af, hvordan batch-effekter påvirker en analyse, er særlig vigtig inden for biologifaget, hvor mange store sekventeringsprojekter involverer data indsamlet eller sekventeret over forskellige batches, sekventeringsmaskiner eller forskellige forberedelsesmetoder.
I skal være i stand til at
- Beskrive randomisering, replikation og blokering
- Beskrive Simpsons paradoks
- Beskrive Anscombes kvartet
- Undersøge for batch-effekter med PCA
- Læs kursusnotaterne om randomisation, replikation, blocking + confounding
- Se videorne til “Simpson’s Paradoks” og “Anscombe’s quartet”
- Læs kursusnotaterne om batch effekt undersøgelse (11.5)
- Quiz - eksperimental (på Absalon)
- Lav problemstillingerne og sige til hvis du har bruge for hjælp :)
11.1.1 Video ressourcer
- Part 1: randomisation, replikation, blocking + confounding
- Ingen video: læs gerne notaterne nedenfor 11.2 “Grundlæggende principper i eksperimentelt design”
- Lav problemstilling Problem 2)
- Part 2: Simpson’s Paradoks
- Se videoen nedenfor samt 11.3 “Case studies: Simpson’s paradoks” i kursusnotaterne
- Lav problemstilling Problem 3)
- Problemstilling Problem 7) giver endnu et eksempel på Simpon’s paradoks
- Part 3: Anscombe’s quartet
- Se videoen nedenfor samt 11.4 “Case studies: Anscombes kvartet” i kursusnotaterne
- Lav Problem 4) i problemstillingerne
- Part 4: Batch effects and principal component analysis.
- Læs 11.5: “Undersøgelse af batch effekts”
- Lav Problem 5) og Problem 6) )
11.2 Grundlæggende principper i eksperimentelt design
11.2.1 Randomisering og replikation
Man laver et eksperiment for at få svar på et bestemt spørgsmål eller en hypotese. Eksperimentet designes ud fra principper, der gør det muligt at fortolke resultaterne fra analysen af datasættet efterfølgende. For at et eksperiment skal være gyldigt, skal det kunne demonstrere hensigtsmæssig replikation og randomisering.
Randomisering
Målet er at udelukke, at konklusionerne simpelthen kan skyldes varians som følge af en faktor, der ikke direkte er interessant i eksperimentet. Derfor bruger man randomisering til at fordele disse faktorer over de forskellige behandlingsgrupper. Et eksempel kan være ‘double-blinding’ i kliniske eksperimenter - både lægen og patienten har ikke kendskab til, hvem der hører til de forskellige grupper, således at forskelsbehandling indenfor grupperne, som kan påvirke de endelige resultater, undgås.
Replikation
Dette sker, når man gentager et eksperiment flere gange - for eksempel ved at have flere patienter i hver behandlingsgruppe. Dette tillader os at beregne variabiliteten i data, som er nødvendig for at kunne konkludere om, der er en forskel mellem grupperne. Man kan altså ikke generalisere resultater, som kun er målt på én person.
Figure 11.1: randomisation og replikation
I ovenstående figur er der 6 replikationer i hver gruppe, der enten er “behandling” eller “kontrol”. I tilfældet “God randomisering” er objekter, som er taget tilfældigt fra populationen, godt matchet mellem de to grupper. I det andet tilfælde, “Dårlig randomisering”, kan man se, at objekternes farver er godt matchet inden for de samme grupper. Dette gør det derfor umuligt at afgøre, om en eventuel forskel mellem “behandling” og “kontrol” i virkeligheden er resultatet af farven frem for de målinger, man gerne vil sammenligne.
Forvekslingsvariabler
Figuren nedenfor illustrerer alder som en forvekslingsvariabel i et eksperiment, hvor man prøver at forstå sammenhængen mellem aktivitetsniveau og vægtøgning. Det kan se ud som om, at et lavt aktivitetsniveau (afhængig variabel) forklarer vægtøgning (uafhængig variabel), men man er nødt til at tage andre variable i betragtning for at sikre, at sammenhængen ikke skyldes noget andet. For eksempel kunne gruppen med det høje aktivitetsniveau bestå af personer, der er yngre end personerne i gruppen med det lave aktivitetsniveau, og deres alder kan påvirke deres vægtøgning (måske på grund af forskelle i stressniveauer, kost osv.).
Blokering
Man kan forsøge at kontrollere for ekstra variable, som vi ikke er interesseret i, gennem “blokering”. Blokering udføres ved først at identificere grupper af individer, der ligner hinanden så meget som muligt. Det kan for eksempel være, at tre forskellige forskere har medvirket til at udføre et stort eksperiment med mange patienter og forskellige behandlingsgrupper. Vi er interesseret i, om der er forskelle mellem behandlingsgrupperne, men ikke i, om der er en forskel i forskernes behandling af patienterne. Derfor vil vi gerne ‘blokere’ efter forsker - altså kontrollere for dem som en “batch” effekt. Man kan også blokere efter fx køn, for at sikre at forskellene i behandlingsgrupperne ikke skyldes forskelle mellem mænd og kvinder. Blokering udføres som en del af en lineær model, efter data er indsamlet, men det er nyttigt at tænke over det fra starten.
11.2.2 Eksempel med datasættet ToothGrowth
Et godt eksempel på et veludført eksperimentelt design er datasættet ToothGrowth, som er baseret på marsvin - de fik forskellige kosttilskud og doser, og derefter blev længden af deres tænder målt.
data(ToothGrowth)
ToothGrowth <- ToothGrowth %>% tibble() %>% mutate(dose = as.factor(dose))
summary(ToothGrowth)#> len supp dose
#> Min. : 4.20 OJ:30 0.5:20
#> 1st Qu.:13.07 VC:30 1 :20
#> Median :19.25 2 :20
#> Mean :18.81
#> 3rd Qu.:25.27
#> Max. :33.90Her kan man se, at for hver gruppe (efter supp og dose) er der 10 marsvin - vi har således replikationer over grupperne, og hver supp (supplement) har hver af de tre mulige værdier for “dose”. Hvis vi for eksempel ikke var interesseret i supp men kun i dosis, kunne vi ‘blokere’ efter supp for at afbøde forskelle i effekten af de to supplementer i supp.
ToothGrowth %>% dplyr::count(supp,dose) %>%
ggplot(aes(x=factor(dose),y=n,fill=factor(dose))) +
geom_bar(stat="identity") +
ylab("Antal marsvin") +
xlab("Dosering") +
facet_grid(~supp) +
theme_bw()
Man skal dog være opmærksom, fordi vi ved ikke, hvordan marsvinene blev tildelt til de forskellige grupper. For eksempel, hvis hanner og hunner ikke er tildelt tilfældigt, kan det forekomme, at supp “OJ” og dosis “0.5” kun har hanmarsvin, og supp “OJ” med dosis “1.0” kun har hunmarsvin. I sådanne tilfælde kunne vi ikke afgøre, om forskellen i dosis “0.5” versus “1.0” er resultatet af dosis eller køn.
11.3 Case studies: Simpsons paradoks
(Se også videoressourcer Del 2).
Simpsons paradoks er et fascinerende statistisk fænomen, der er en vigtig påmindelse om, at korrelation ikke nødvendigvis indebærer kausalitet, og at det er afgørende at forstå de underliggende data og ikke kun de aggregerede statistikker. Simpsons paradoks opstår, når man drager to modsatte konklusioner fra det samme datasæt - på den ene side, når man kigger på dataene samlet, og på den anden side, når man tager visse grupper i betragtning. Vi kan visualisere Simpsons paradoks gennem eksemplet nedenfor - her har vi to variable x og y, som vi kan bruge til at lave et scatterplot, samt nogle forskellige grupper inden for variablen group.
#library(datasauRus)
simpsons_paradox <- read.table("https://www.dropbox.com/s/ysh3qpc7qv0ceut/simpsons_paradox_groups.txt?dl=1",header=T)
simpsons_paradox <- simpsons_paradox %>% tibble()
simpsons_paradoxFALSE # A tibble: 222 × 3
FALSE x y group
FALSE <dbl> <dbl> <chr>
FALSE 1 62.2 70.6 D
FALSE 2 52.3 14.7 B
FALSE 3 56.4 46.4 C
FALSE 4 66.8 66.2 D
FALSE 5 66.5 89.2 E
FALSE 6 62.4 91.5 E
FALSE 7 38.9 6.76 A
FALSE 8 39.4 63.1 C
FALSE 9 60.9 92.6 E
FALSE 10 56.6 45.8 C
FALSE # ℹ 212 more rowsHvis vi ignorerer group og kigger på dataene samlet, kan vi se, at der er en stærk positiv sammenhæng mellem x og y. Men når vi opdeler efter de forskellige grupper ved at skrive colour = group, opnår vi faktisk en negativ sammenhæng indenfor hver af grupperne.
p1 <- simpsons_paradox %>%
ggplot(aes(x,y)) +
geom_point() +
geom_smooth(method="lm",se=FALSE) +
theme_classic()
p2 <- simpsons_paradox %>%
ggplot(aes(x,y,colour=group)) +
geom_point() +
geom_smooth(method="lm",aes(group=group),colour="black",se=FALSE) +
theme_classic()
library(gridExtra)
grid.arrange(p1,p2,ncol=2)
Simpsons paradoks forekommer oftere end man skulle tro, og derfor er det vigtigt at overveje, hvilke andre variable man også er nødt til at tage i betragtning.
11.3.1 Optagelse på Berkeley
Det mest berømte eksempel på Simpsons paradoks drejer sig om optagelsen på Berkeley Universitetet i 1973. Følgende tabel fra Wikipedia (https://en.wikipedia.org/wiki/Simpson%27s_paradox) viser statistikker over antallet af ansøgere samt procentdelen, der blev optaget på universitetet generelt, opdelt efter køn.
Figure 11.2: source: wikipedia
Hvis vi laver et søjlediagram af tallene, kan man se, at der er en højere procentdel af mænd end kvinder, som blev optaget på universitetet (sagen medførte en retssag mod universitetet).
admissions_all <- tibble("sex"=c("all","men","women"),admitted=c(41,44,35))
admissions_all %>% ggplot(aes(x=sex,y=admitted,fill=sex)) +
geom_bar(stat="identity") +
theme_minimal() +
ylab("Procent optaget") +
scale_x_discrete(limits = c("women","men","all")) +
coord_flip()
Da man dog kiggede lidt nærmere på de samme tal, men opdelte efter de forskellige fakulteter på universitetet, fik man et anderledes billede af situationen. I følgende tabel har vi optagelsestallene for mænd og kvinder på hver af de forskellige fakulteter (A til F).
admissions_separate <- tribble(
~department, ~all, ~men, ~women,
#------------|-------|-------|--------#
"A", 64, 62, 82,
"B", 63, 63, 68,
"C", 35, 37, 34,
"D", 34, 33, 35,
"E", 25, 28, 24,
"F", 6, 6, 7
)Man kan se, at for de fleste af fakulteterne er der ikke en markant forskel mellem mænd og kvinder, og i nogle tilfælde havde kvinder faktisk en større sandsynlighed for at blive optaget.
admissions_separate %>%
pivot_longer(-department,names_to="sex",values_to="admitted") %>%
ggplot(aes(x=department,y=admitted,fill=sex)) +
ylab("Procent optaget") +
geom_bar(stat="identity",position = "dodge",colour="black") +
theme_minimal()
Hvad skyldes denne sammenhæng? Det viste sig, at kvinder havde en tendens til at ansøge indenfor de fakulteter, som var sværest at komme ind på. For eksempel kan man se her, at fakultet E har en relativt lav optagelsesprocent. Det samme fakultet var dog et af dem, hvor betydeligt flere kvinder ansøgte end mænd.
applications_E <- tibble("sex"=c("all","men","woman"),applications=c(584,191,393))
applications_E %>% ggplot(aes(x=sex,y=applications,fill=sex)) +
geom_bar(stat="identity") +
theme_minimal() +
ylab("Number of applications to dep. E") +
scale_x_discrete(limits = c("woman","men","all")) +
coord_flip()
Derfor, selvom kvinder ikke havde en lavere sandsynlighed for at blive optaget end mænd i deres fortrukne fag, var antallet af kvinder, der blev optaget i det hele taget på tværs af alle afdelinger, faktisk lavere end antallet af mænd.
11.4 Case studies: Anscombes kvartet
(Se også videoressourcer Del 3).
Anscombes kvartet (se også https://en.wikipedia.org/wiki/Anscombe%27s_quartet) er et meget nyttigt og berømt eksempel fra 1973, der fremhæver vigtigheden af at visualisere datasættet. Vi kan hente dataene fra linket nedenfor - der er x-værdier og y-værdier, som kan anvendes til at lave et scatterplot, og der er også set, der refererer til fire forskellige datasæt (derfor ‘kvartet’).
anscombe <- read.table("https://www.dropbox.com/s/mlt7crdik3eih9a/anscombe_long_format.txt?dl=1",header=T)
anscombe <- anscombe %>% tibble()
anscombe#> # A tibble: 44 × 3
#> set x y
#> <int> <int> <dbl>
#> 1 1 10 8.04
#> 2 1 8 6.95
#> 3 1 13 7.58
#> 4 1 9 8.81
#> 5 1 11 8.33
#> 6 1 14 9.96
#> 7 1 6 7.24
#> 8 1 4 4.26
#> 9 1 12 10.8
#> 10 1 7 4.82
#> # ℹ 34 more rowsFormålet med datasættet er, at vi gerne vil fitte en lineær regressionsmodel for at finde den forventede y-værdi afhængig af x (husk lm(y ~ x)). Da vi har fire datasæt, kan vi opdele datasættet efter set og anvende funktionerne nest og map (se Kapitel 7) til at fitte de fire lineære regressionsmodeller. Vi anvender også tidy og glance for at få resuméstatistikker fra de fire modeller:
my_func <- ~lm(y ~ x, data = .x)
tidy_anscombe_models <- anscombe %>%
group_nest(set) %>%
mutate(fit = map(data, my_func),
tidy = map(fit, tidy),
glance = map(fit, glance))Vi kan anvende unnest på outputtet fra tidy og se på skæringspunktet og hældningen af de fire modeller. Man kan se, at de to parametre er næsten identiske for de fire modeller:
tidy_anscombe_models %>% unnest("tidy") %>%
pivot_wider(id_cols = "set",names_from = "term",values_from="estimate") #> # A tibble: 4 × 3
#> set `(Intercept)` x
#> <int> <dbl> <dbl>
#> 1 1 3.00 0.500
#> 2 2 3.00 0.5
#> 3 3 3.00 0.500
#> 4 4 3.00 0.500Hvad med de andre parametre fra modellen - lad os eksempelvis kigge på r.squared og p.value fra modellerne, som kan findes i outputtet fra glance. Her kan vi igen se, at de er næsten identiske.
#> # A tibble: 4 × 3
#> set r.squared p.value
#> <int> <dbl> <dbl>
#> 1 1 0.667 0.00217
#> 2 2 0.666 0.00218
#> 3 3 0.666 0.00218
#> 4 4 0.667 0.00216Hvad med korrelation? Den er også næsten den samme:
my_func <- ~cor(.x$x,.x$y)
anscombe %>%
group_nest(set) %>%
mutate(cor = map(data, my_func)) %>%
unnest(cor) %>%
select(-data)#> # A tibble: 4 × 2
#> set cor
#> <int> <dbl>
#> 1 1 0.816
#> 2 2 0.816
#> 3 3 0.816
#> 4 4 0.817Kan vi så konkludere, at de fire datasæt, som underbygger de forskellige modeller, er identiske? Lad os lave et scatter plot af de fire datasæt (som vi faktisk burde have gjort i starten af vores analyse).
anscombe %>%
ggplot(aes(x = x, y = y,colour=factor(set))) +
geom_point() +
facet_wrap(~set) +
geom_smooth(method = "lm", se = FALSE) +
theme_minimal()
De fire datasæt er meget forskellige. Vi ved, at de alle har samme bedste tilpasning med rette linjer, men de underliggende data er slet ikke de samme. Det første datasæt ser egnet ud til en lineær regressionsanalyse, men vi kan se i datasæt nummer to, at der ikke engang er en lineær sammenhæng. Og de andre to har outlier værdier, hvilket gør, at den bedst tilpassede rette linje ikke passer særlig godt til punkterne.
11.5 Undersøgelse af “batch-effekter”
En “batch-effekt” er en systematisk teknisk bias, der opstår, når målinger (f.eks. genekspression) udføres i flere omgange eller “batches”. Dette kan potentielt føre til betydelige forvrængninger og fejlagtige konklusioner i dataanalyse, hvis det ikke tages højde for.
Forestil dig, at du udfører et eksperiment, hvor du sammenligner genekspressionen i sundt væv med kræftvæv. Du udfører dine eksperimenter over flere dage, og det viser sig, at alle dine sunde prøver blev behandlet på mandag, og alle dine kræftprøver blev behandlet på tirsdag. Hvis der er en systematisk forskel i, hvordan dine eksperimenter blev udført på de to dage (måske et reagens var lidt anderledes, eller instrumentet blev brugt på en anden måde), vil du se en stor forskel i genekspression mellem dine sunde og kræftprøver. Men denne forskel skyldes faktisk batch-effekten, ikke forskellen mellem sundt og kræftvæv.
Man kan også anvende visualiseringer til at undersøge eventuelle batch-effekter eller forvirrende variabler i datasættet. Dette er især vigtigt i store eksperimenter, hvor forskellige prøver eller dele af datasættet bliver indsamlet på forskellige tidspunkter, steder, eller af forskellige personer. Det er ofte tilfældet i sekvenseringsbaserede datasæt, at man ser batch-effekter, og det kan skyldes mange ting, bl.a.:
- Sekvenseringsdybde
- Grupper af prøver lavet på forskellige tidspunkter af forskellige individer
- Sekvenseringsmaskiner - prøver sekvenseret på forskellige maskiner eller ‘lanes’.
Lad os tage udgangspunkt i nogle genekspressionssekvenseringsdata fra mus (vi så også dette datasæt, da vi lærte om pivot_longer kombineret med left_join).
norm.cts <- read.table("https://www.dropbox.com/s/3vhwnsnhzsy35nd/bottomly_count_table_normalised.txt?dl=1")
coldata <- read.table("https://www.dropbox.com/s/el3sm9ncvzbq6xf/bottomly_phenodata.txt?dl=1")
coldata <- coldata %>% tibble()
norm.cts <- as_tibble(norm.cts,rownames="gene")Jeg begynder med at vælge kun de rækker, der har mindst 50 counts, for at undgå gener med lave ekspressionsniveauer. Det næste jeg gør er at transformere dataene til logaritmisk form (funktion map_if()) for at opnå en bedre fordeling i datasættet. Derefter tager jeg kun de første 1000 rækker med, kun for at gøre det nemmere at køre mine koder på en laptop.
#normalisere og filtrere dataene
norm.cts <- norm.cts %>%
filter(rowSums(norm.cts %>% select(-gene))>50) %>%
map_if(is.numeric,~log(.x+1)) %>% as_tibble()
norm.cts <- norm.cts[1:1000,] #take a subset to save loading times on laptopSå der er 1000 gener i rækkerne, og der er 21 forskellige prøver, der spreder sig over kolonnerne. Vi har også nogle prøveoplysninger - der er to forskellige stammer af mus og også forskellige batches, som vi gerne vil undersøge nærmere.
#> # A tibble: 21 × 5
#> column num.tech.reps strain batch lane.number
#> <chr> <int> <chr> <int> <int>
#> 1 SRX033480 1 C57BL.6J 6 1
#> 2 SRX033488 1 C57BL.6J 7 1
#> 3 SRX033481 1 C57BL.6J 6 2
#> 4 SRX033489 1 C57BL.6J 7 2
#> 5 SRX033482 1 C57BL.6J 6 3
#> 6 SRX033490 1 C57BL.6J 7 3
#> 7 SRX033483 1 C57BL.6J 6 5
#> 8 SRX033476 1 C57BL.6J 4 6
#> 9 SRX033478 1 C57BL.6J 4 7
#> 10 SRX033479 1 C57BL.6J 4 8
#> # ℹ 11 more rowsVi kan se på, hvor mange prøver vi har for hver kombination af stamme og batch i datasættet:
#>
#> 4 6 7
#> C57BL.6J 3 4 3
#> DBA.2J 4 3 4Så man kan se, at både stammen er repræsenteret med tre eller fire prøver i hver af de tre batches. Der er derfor replikation, og da vi har fået repræsenteret hver kombination af stamme og batch, kan vi eventuelt blokere efter batch for at fjerne dens effekt. Her har vi ikke tid til at se på metoder til at fjerne batch-effekter, men det er vigtigt, at vi er i stand til at opdage dem.
11.5.1 Principal component analyse
Man kan undersøge mulige batch-effekter via principal component analyse.
Først skal vi lave en “transpose” af datasættet således at variablerne (samples) bliver i rækkerne og observationerne (genes) er som kolonnerne.
Det er ikke helt ligetil i tidyverse, men her er en måde at gøre det på:
norm.cts.transpose <- as_tibble(t(norm.cts[,-1])) #lave transpose med t() funktion (drop genenavne)
colnames(norm.cts.transpose) <- norm.cts %>% pull(gene) #tilføje gene som kolonnenavne
norm.cts.transpose %>% mutate(sample = colnames(norm.cts)[-1],.before=1) #sample skal være den første kolonne#> # A tibble: 21 × 1,001
#> sample ENSMUSG00000000001 ENSMUSG00000000056 ENSMUSG00000000058
#> <chr> <dbl> <dbl> <dbl>
#> 1 SRX033480 6.35 3.51 3.19
#> 2 SRX033488 6.32 3.56 3.50
#> 3 SRX033481 6.21 3.57 3.08
#> 4 SRX033489 6.29 3.27 3.21
#> 5 SRX033482 6.31 2.99 3.14
#> 6 SRX033490 6.27 3.61 3.09
#> 7 SRX033483 6.30 3.56 3.22
#> 8 SRX033476 6.03 3.69 2.83
#> 9 SRX033478 6.30 3.42 3.73
#> 10 SRX033479 5.88 3.98 3.24
#> # ℹ 11 more rows
#> # ℹ 997 more variables: ENSMUSG00000000078 <dbl>, ENSMUSG00000000088 <dbl>,
#> # ENSMUSG00000000093 <dbl>, ENSMUSG00000000120 <dbl>,
#> # ENSMUSG00000000125 <dbl>, ENSMUSG00000000126 <dbl>,
#> # ENSMUSG00000000127 <dbl>, ENSMUSG00000000131 <dbl>,
#> # ENSMUSG00000000149 <dbl>, ENSMUSG00000000167 <dbl>,
#> # ENSMUSG00000000168 <dbl>, ENSMUSG00000000171 <dbl>, …Heldigvis er der en pakke, som kan gøre processen meget hurtigere, som vi kan benytte os af her:
library(sjmisc)
norm.cts.transpose <- norm.cts %>%
rotate_df(cn=TRUE,rn="sample") %>%
as_tibble()
norm.cts.transpose#> # A tibble: 21 × 1,001
#> sample ENSMUSG00000000001 ENSMUSG00000000056 ENSMUSG00000000058
#> <chr> <dbl> <dbl> <dbl>
#> 1 SRX033480 6.35 3.51 3.19
#> 2 SRX033488 6.32 3.56 3.50
#> 3 SRX033481 6.21 3.57 3.08
#> 4 SRX033489 6.29 3.27 3.21
#> 5 SRX033482 6.31 2.99 3.14
#> 6 SRX033490 6.27 3.61 3.09
#> 7 SRX033483 6.30 3.56 3.22
#> 8 SRX033476 6.03 3.69 2.83
#> 9 SRX033478 6.30 3.42 3.73
#> 10 SRX033479 5.88 3.98 3.24
#> # ℹ 11 more rows
#> # ℹ 997 more variables: ENSMUSG00000000078 <dbl>, ENSMUSG00000000088 <dbl>,
#> # ENSMUSG00000000093 <dbl>, ENSMUSG00000000120 <dbl>,
#> # ENSMUSG00000000125 <dbl>, ENSMUSG00000000126 <dbl>,
#> # ENSMUSG00000000127 <dbl>, ENSMUSG00000000131 <dbl>,
#> # ENSMUSG00000000149 <dbl>, ENSMUSG00000000167 <dbl>,
#> # ENSMUSG00000000168 <dbl>, ENSMUSG00000000171 <dbl>, …Nu kan vi udføre en principal component analyse:
pca_fit <- norm.cts.transpose %>%
select(where(is.numeric)) %>% # behold kun numeriske kolonner
prcomp(scale = TRUE) # udfør PCA på skalerede dataSå kan vi bruge vores standarde måde at visualisere principal components på:
library(ggrepel)
pca_fit_augment <- pca_fit %>%
augment(norm.cts.transpose)
pca_fit_augment %>%
ggplot(aes(x=.fittedPC1,y=.fittedPC2)) +
geom_point() +
geom_text_repel(aes(label=sample)) +
theme_bw()
Anvend left_join til at få oplysningen med:
Når der er så mange kolonnerne i pca_fit_augment bruger jeg tail-funktionen ovenpå colnames så jeg kan se, hvad for nogle oplysningsvariabler vi har fået med:
#> [1] ".fittedPC20" ".fittedPC21" "num.tech.reps" "strain"
#> [5] "batch" "lane.number"Vi er interesseret i at undersøge mulige batch-effektor, så lad os give en forskellige farve efter batch og en forskellige form efter strain:
pca_fit_augment %>%
ggplot(aes(x=.fittedPC1,y=.fittedPC2,colour=as.factor(batch),shape=strain)) +
geom_point(size=3) +
geom_text_repel(aes(label=sample)) +
theme_bw()
Man kan også se på dataene på en anden måde ved at lave boksdiagrammer for de to første principale komponenter, opdelt efter batch. Vi får bekræftet vores observation om, at der er en markant forskel mellem batch 4 og de andre to batches, og det er et problem, som bør korrigeres, før man foretager yderligere analyser af dataene.
p1 <- pca_fit_augment %>%
ggplot(aes(x=factor(batch),y=.fittedPC1,fill=factor(batch))) + geom_boxplot(show.legend = F) + geom_jitter(show.legend = F) + theme_minimal() + ggtitle("Batches PC1")
p2 <- pca_fit_augment %>%
ggplot(aes(x=factor(batch),y=.fittedPC2,fill=factor(batch))) + geom_boxplot(show.legend = F) + geom_jitter(show.legend = F) + theme_minimal() + ggtitle("Batch PC2")
library(gridExtra)
grid.arrange(p1,p2,ncol=2)
11.5.2 EKSTRA: Heatmap
Man kan også lav en heatmap for at undersøge batch effektor. Jeg begynder med at lave en korrelationsmatrix af mine data, som indeholder sammenhængene (pearson’s korrelationskoefficient) mellem de forskellige variabler (husk at drop den første kolonne, da den ikke er numeriske):
Jeg laver en heatmap med heatmap-funktionen i ovenstående:
Samples som har en stækere sammenhænge forekommer tættere på hinanden i rækkerne/kolonnerne. Her kan man se, at samples fra batch 7 er mere ofte sammen oppe på de ydereste rækkere i forhold til batch 4, hvor samples er mere ofte på de nederste rækker.
11.5.3 EKSTRA: Batch effect correction med ComBat()-funktionen
sva-pakken i R har en funktionen der hedder ComBat, der kan hjælpe med at korrigere for batch-effekter, så forskelle mellem prøver skyldes biologisk variation og ikke tekniske forskelle.
Se her for ydereligere om ComBat: https://rdrr.io/bioc/sva/man/ComBat.html
#if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("sva")
library(sva)
corrected <- sva::ComBat(dat=norm.cts %>% select(-gene),batch=coldata %>% pull(batch))
corrected <- as_tibble(corrected) %>% mutate(gene = norm.cts %>% pull(gene),.before=1)Så kan man lave præcis samme procedure som i ovenstående på den korrigede dataframe:
corrected.transpose <- corrected %>%
rotate_df(cn=TRUE,rn="sample")
pca_fit <- corrected.transpose %>%
select(where(is.numeric)) %>% # behold kun numeriske kolonner
prcomp(scale = TRUE) # udfør PCA på skalerede data
pca_fit_augment <- pca_fit %>%
augment(corrected.transpose)
pca_fit_augment <- pca_fit_augment %>%
left_join(coldata %>% rename(sample = column),by="sample")
pca_fit_augment <- pca_fit_augment %>% mutate(batch=as.factor(batch))Og så lav en plot:
PCA1 <- pca_fit_augment %>%
ggplot(aes(x=.fittedPC1,y=.fittedPC2,colour=batch,shape=strain)) +
geom_point(size=3) +
geom_text_repel(aes(label=sample)) +
theme_bw()
PCA2 <- pca_fit_augment %>%
ggplot(aes(x=.fittedPC1,y=.fittedPC2,shape=batch,colour=strain)) +
geom_point(size=3) +
geom_text_repel(aes(label=sample)) +
theme_bw()
grid.arrange(PCA1,PCA2,ncol=2)
–>
11.6 Problemstillinger
Problem 1) Quiz på Absalon - experimental
Problem 2) Eksperimentelt design
Jeg udfører et eksperiment, hvor patienter modtager et af tre forskellige kosttilskud (Gruppe 1, 2 og 3). Der er 5 patienter i hver gruppe, og jeg ønsker at undersøge, om patienternes energiniveau i gennemsnit varierer mellem de tre grupper. Alderen på patienterne i hver af de tre grupper er:
Gruppe 1: 18, 23, 31, 25, 19
Gruppe 2: 24, 29, 35, 21, 30
Gruppe 3: 43, 52, 33, 39, 40
Hvad er problemet med det eksperimentelle design her? Lav boxplots for at illustrere fordelingen af alderne for hver af de tre grupper (du skal starte med at lave en tibble, der indeholder dataene - se tilbage til Emne 5.5 hvis du har glemt, hvordan man laver en tibble).
Hvis man finder en signifikant forskel mellem de tre kosttilskud, kan man så stole på resultaterne?
Hvilke andre variabler end alder kunne være årsag til en eventuel forskel mellem de tre kosttilskud, og som muligvis skulle tages i betragtning?
Hvad kan man gøre for at løse problemet med det ovenstående eksperimentelle design?
Problem 3) Simpsons paradoks Genbesøg af Lung Cap-data
Indlæs LungCapData og tilføj den kategoriske variabel Age.Group:
LungCapData <- read.csv("https://www.dropbox.com/s/ke27fs5d37ks1hm/LungCapData.csv?dl=1")
LungCapData$Age.Group <- cut(LungCapData$Age,breaks=c(1,13,15,17,19),right=FALSE,include.lowest = TRUE)
levels(LungCapData$Age.Group) <- c("<13","13-14","15-16","17+")a) Lav boxplots med smoke på x-aksen og LungCap på y-aksen.
+ Bemærk hvilken gruppe der har den højeste lungekapacitet.
b) Lav samme plot, men adskilt efter Age.Group, og beskriv, hvordan det er et eksempel på Simpsons paradoks.
c) Lav et boxplot med Age på y-aksen og Smoke på x-aksen for at understøtte forklaringen på, hvorfor man ser Simpsons paradoks i dette datasæt.
Problem 4) Anscombes analyse Gentag Anscombes analyse med dinosaurdatasættet:
FALSE # A tibble: 1,846 × 3
FALSE dataset x y
FALSE <chr> <dbl> <dbl>
FALSE 1 dino 55.4 97.2
FALSE 2 dino 51.5 96.0
FALSE 3 dino 46.2 94.5
FALSE 4 dino 42.8 91.4
FALSE 5 dino 40.8 88.3
FALSE 6 dino 38.7 84.9
FALSE 7 dino 35.6 79.9
FALSE 8 dino 33.1 77.6
FALSE 9 dino 29.0 74.5
FALSE 10 dino 26.2 71.4
FALSE # ℹ 1,836 more rowsa) Fit en lineær regressionsmodel for hvert af datasættene (brug group_by, nest og map i kombination med en custom funktion), hvor y er den afhængige variabel og x er den uafhængige variabel. Anvend også tidy og glance på alle modellerne.
b) Brug resultaterne fra tidy til at undersøge hældning og skæring med y-aksen for de forskellige modeller - ligner de hinanden?
c) Brug også resultaterne fra glance til at undersøge r.squared og p.value.
d) Er de alle det samme datasæt? Lav et scatter plot opdelt efter de forskellige datasæt. Hvordan ser de bedste rette linjer ud på plottene?
Problem 5) Vi vil gerne undersøge eventuelle batcheffekter i det følgende datasæt. Det er simulerede “single cell”-sekvenserings count data (dataframen cse50) samt dataframen batches, som angiver hvilken batch hver af de 500 celler tilhører.
cse50 <- read.table("https://www.dropbox.com/s/o0wzojpcsekeg6z/cell_mix_50_counts.txt?dl=1")
batches <- read.table("https://www.dropbox.com/s/4t382bfgro46ka5/cell_mix_50_batches.txt?dl=1")
batches <- tibble("batch"=batches %>% pull(batch),cell=paste0("cell_",c(1:500)))
cse50 <- as_tibble(cse50,rownames="gene") %>% mutate(gene = paste0("gene_",gene))a) Transformere værdierne i cse50 til log-skala (tilføj 1 først og tag logaritmen bagefter).
b) For at undersøge mulige batch-effektor bruger vi samme process som i kursusnotaterne ovenstående. Lav en transpose af cse50 således at cellerne er i rækkerne og generne er i kolonnerne.
c) Udfør en prinicipal component analyse og tilføj din transposed-datasæt til dit resultat
d) Tilføj batch oplysning til din nye dataframe fra c (som har dine principal components).
e) Lav et tilpassende plot af dine første to principal components, hvor du farver punkterne efter de forskellige batches.
f) Som alternativ visualisering lav også boxplots for de første to prinpical components fordelt efter batch. Kommenter kort på eventuelle batch effekts i datasættet ud fra e) og f).
Problem 6) Som tilføjelse til Problem 5), prøve
a) Lav en heatmap af cellerne i datasættet med heatmap-funktionen()
b) Anvend ComBat-function på datasættet (husk at installere sva-pakken - når den ikke er på CRAN se kode for gøre det i notaterne ovenpå) og gentage samme undersøgelsen som du lavede i Problem 6 for at se, om batch-effekterne findes stadig i datasættet.
Problem 7) Yderligere Simpson’s paradoks
Kør følgende kode for at indlæse og bearbejde det følgende datasæt airlines.
airlines <- read.table("http://www.utsc.utoronto.ca/~butler/d29/airlines.txt",header=T)
airlines <- airlines %>%
pivot_longer(-airport) %>%
separate(name,sep="_",into = c("airline","status")) %>%
mutate(airline = recode(airline, aa = "Alaska", aw = "American")) %>%
pivot_wider(names_from=status,values_from=value) %>%
mutate("ontime" = ontime + delayed) %>%
rename(flights = ontime)a) For hvert flyselskab opsummer antallet af flights og antallet af delayed.
b) Beregn også andelen (proportionen) af flyvninger, der er forsinkede, i hvert flyselskab. Lav et søjlediagram for at vise proportionerne.
c) Denne gang, lave samme opsummerings beregninger som i a og b men til hver lufthavn (i stedet for hvert flyselskab). Lav et plot.
d) Denne gang, beregn andelen af flyvninger, der er forsinkede, for hver kombination af både lufthavn og flyselskab. Omsæt igen dette til et plot.
e) Kan du forklare dette? Hint: kig eksempelvis på rådataene og især lufthavnen “Phoenix”.
11.7 Yderligere læsning
Simpson’s paradox og airlines: http://ritsokiguess.site/docs/2018/04/07/simpson-s-paradox-log-linear-modelling-and-the-tidyverse/
Batch effekt correction: https://en.wikipedia.org/wiki/Batch_effect#Correction