---
title: 'Genoomgrootte'
output:
  html_document:
    code_download: yes
    highlight: tango
    number_sections: yes
    theme: cosmo
    toc: yes
    toc_float: yes
  word_document:
    toc: yes
---

# Genoomgrootte

Sinds de jaren '70 worden de genomen van verscheidene organismen bestudeerd. Deze zijn niet enkel de sleutel tot de evolutietheorie, maar geven ook een beeld van ziekten veroorzaakt door genetische afwijkingen, laten het ontdekken van nieuwe functies van bacteriën toe of kunnen de taxonomie van organismen drastisch wijzigen.
In de vroege jaren 2000 was er een stevige opmars van "next-generation sequencing" methoden die het mogelijk maakten om complete genomen te sequencen en te bestuderen. Dit maakte het ook mogelijk om een predictie te doen over de genoomgrootte. Hierbij kwam men er al snel achter dat de complexiteit van een organisme niet meteen gelinkt is met de grootte van zijn genoom. Sommige ééncelligen hebben bijvoorbeeld meer DNA dan mensen.

Bij sommige groepen van eukaryoten, zoals de zoogdieren (mammalia), wordt de grootte van het genoom niet enkel beïnvloed door genetische factoren, maar ook door omgevingsfactoren. Zo zou men de vraag kunnen stellen of het voedingspatroon van dieren gelinkt is met de grootte van het genoom.

Genoomgrootte wordt typisch uitgedrukt in picograms (massa), ook wel de C-waarde genoemd. Een database die alle, momenteel gekende (2017), waarden voor genoomgrootten bevat is de "Animal genome size library". Van deze database hebben wij een subset gemaakt die enkel de zoogdieren bevat (totaal van 810 observaties). De dataset bestaat uit de volgende variabelen:

- `Common name`: de naam van het organisme
- `C value`: maat voor de grootte van het genoom
- `Order`: Indelingsmaat voor de classificatie (niveau: Orde)

# Libraries inladen
```{r}
library(ggplot2)
library(dplyr)
#install.packages("tidyr")
library(tidyr)
#install.packages("multcomp")
library(multcomp)
library(readr)
```

# Dataset genome_sizes inlezen

```{r}
# Importeer de data
genomes <- read_delim("https://raw.githubusercontent.com/statOmics/statistiekBasisCursusData/master/practicum5/genome_sizes.txt", delim = "\t", col_types = list(col_character(), col_double(), col_factor()))
#Vervang de spaties in de namen met een punt.
colnames(genomes) <-  make.names(colnames(genomes))

```

# Onderzoeksvraag: Is er een effect van voedingspatroon?

We zijn geinteresseerd of er een link is tussen het voedingspatroon van zoogdieren en de grootte van hun genoom. Hierbij wensen we volgende groepen (Orden) binnen de zoogdieren te bestuderen:

- Primaten (omnivoren)
- Carnivora (carnivoren)
- Artiodactyla (=evenhoevingen, herbivoren).

Om deze vraag te beantwoorden, zullen we eerst een subset opstellen van deze soorten. Deze subset noemen we `interest1`. Vanaf nu zullen we met deze subset verder werken.

```{r}
interest1 <- genomes %>% filter(Order %in% c("Artiodactyla","Carnivora","Primates"))

#De dataset filteren past niet de levels aan van de factors. Dit passen we hier manueel aan bij Order
interest1<-interest1 %>% mutate(Order=factor(Order,levels= c("Primates", "Carnivora", "Artiodactyla")))
```

# Vraag 1: Inspecteer de variabelen Common.Name, C.value en Order aan de hand van een frequentietabel (leesopdracht)

Hieronder wordt een samenvattende tabel weergegeven van de variabelen in de dataset `interest1`. Voor de variabelen van interesse, namelijk Orde en genoomgrootte (C.value). Voor een categorische variabele kiezen we voor `table()`, bij een continue variabele houdt zo'n tabel geen steek en opteren we voor `range()`.

Op basis van de frequentietabel zien we bijvoorbeeld dat er meer dan dubbel zoveel waarnemingen zijn van de primaten in vergelijking met de carnivoren en de evenhoevigen.

```{r}
summary(interest1)
interest1 %>% count(Order) # categorische variabele
range(interest1$C.value) # continue variabele
```

# Vraag 2: Data exploratie (boxplot)

Gebruik nu een ggplot visualisatie om een boxplot te maken van de genoomgrootte per orde. Toon de individuele observaties op de boxplot met een jitter en geef de boxplots een kleur per orde Sla het resultaat op in het object `boxplot`.

```{r}
boxplot <- ggplot(data=...,aes(...)) +
  ... + # boxplot maken. Gebruik het argument outlier.shape = NA
  ... + 
  theme_bw() +
  ggtitle("Genoomgrootte per orde")
boxplot
```


# Vraag 3: De hypotheses (Meerkeuzevraag)
Stel dat we de primaten gelijk stellen aan groep 1, carnivoren aan groep 2 en Artiodactyla aan groep 3. De gemiddelde genoomgroottes van deze groepen stellen we respectievelijk gelijk aan $\mu_1$, $\mu_2$ en $\mu_3$. De nulhypothese kan dan geformuleerd worden als:

1. $H_0$: De gemiddelde genoomgrootte van de verschillende groepen is gelijk aan nul, ofwel $\mu_1=\mu_2=\mu_3=0$
2. $H_0$: De gemiddelde genoomgrootte is gelijk tussen de verschillende groepen, ofwel $\mu_1=\mu_2=\mu_3$
3. $H_0$: Er is geen verband tussen de verschillende behandelingsgroepen
4. $H_0$: Het intercept is gelijk aan nul.

Geef in Dodona het nummer van het correcte antwoord.

# Vraag 4: De hypotheses (Meerkeuzevraag)

De alternatieve hypothese kan dan geformuleerd worden als

1. $H_A$: De gemiddelde genoomgrootte van de primaten is verschillend van de twee andere ordes.
2. $H_A$: De gemiddelde genoomgrootte van minstens 2 ordes is verschillend van elkaar.
3. $H_A$: De gemiddelde genoomgrootte van alle ordes is verschillend van elkaar.
4. $H_A$: De gemiddelde genoomgrootte van alle ordes is gelijk aan elkaar.

Geef in Dodona het nummer van het correcte antwoord.

# Vraag 5: Fit het model
Fit een lineair model met de variabelen `Order`en `C.value`. Bepaal welke de (on)afhankelijke variabele is.
```{r}
model = lm(...)
summary(model)
```

# Vraag 6: Opmerking lineair model (leesopdracht)
Het lineair model uit de vorige vraag zag er zo uit:

```{r}
summary(model)
```

Volgens de output lijken er verschillen te zijn tussen de genoomgroottes van de drie verschillende groepen. Let wel op: hier is nog geen correctie gebeurd voor meervoudig testen!

# Vraag 7: Voorwaarden ANOVA (leesopdracht)

We zullen de assumpties van de ANOVA test nagaan. We weten dat de volgende vier voorwaarden voldaan moeten zijn alvorens de test uit te voeren:

- normale verdeling binnen elke groep
- elke groep heeft een gelijke variantie
- De gegevens zijn onafhankelijk 

De data uit de database is een verzameling van data van verschillende studies, die onafhankelijk van elkaar verzameld zijn. Echter, sommige soorten zijn meermaals geobserveerd (bv. binnen de groep van de carnivoren, zien we bijvoorbeeld dat er 10 waarnemingen zijn van honden, 5 van katten, 2 van de Bengaalse tijger en dat van de andere soorten slechts 1 waarneming aanwezig is). Het is sterk te verwachten dat de genoomgroottes van twee honden meer gelijkaardig zullen zijn dan de genoomgroottes van een hond versus een kat. Er zit dus een additionele correlatiestructuur in de data waardoor de gegevens binnen elke groep niet onafhankelijk zijn.

Om didactische redenen zullen we de test echter toch uitvoeren.

# Vraag 8: Homoscedasticiteit nagaan

Ga de voorwaarde van de homoscedasticiteit na door middel van een boxplot. Het is niet nodig om de boxplots nogmaals aan te maken. We kunnen de boxplot van vraag 2 gebruiken.

```{r}
boxplot

```

# Vraag 9: Conclusie homoscedasticiteit (leesopdracht)

Op basis van de boxplots kunnen we besluiten dat er geen reden is om aan te nemen dat de populatievarianties verschillend zullen zijn (de groottes van de boxen zijn vrij gelijkaardig). We zullen de homoscedasticiteitsassumptie dus aanvaarden.

# Vraag 10: Normaliteit nagaan
Ga na of ook de assumptie van normaliteit binnen de verschillende groepen hier voldaan is.

```{r}
plot_qq <- interest1 %>% ggplot(aes(sample = ...)) +
                        ... + # qq-punten
                        ... + # qq-lijn
                        theme_bw() +
                        facet_wrap(~Order) #qq-plot per Orde
plot_qq
```

# Vraag 11: Conclusie normaliteit binnen elke groep (leesopdracht)

De gegevens voor de Artiodactyla lijken normaal verdeeld. Bij de andere groepen zijn er afwijkingen van een normale verdeling. Voor de primaten is de data wel licht scheef naar rechts verdeeld (de datapunten op de QQ-plot liggen langs de rechterkant boven de eerste bissectrice; de geobserveerde waarden liggen dus hoger dan men zou verwachten onder de normale). Omdat het aantal observaties hier > 100 kunnen we beroep doen op de centrale limietstelling. Deze zegt dat bij een voldoende groot aantal observaties het gemiddelde van de observaties bij benadering een normale verdeling volgt.

**Merk op:** we kunnen ons ook de vraag stellen of deze database een goede voorstelling is van de werkelijke populatie. Als we binnen de groep van de carnivoren kijken, zien we bijvoorbeeld dat er 10 waarnemingen zijn van honden, 5 van katten, 2 van de Bengaalse tijger en dat van de andere soorten slechts 1 waarneming aanwezig is. Het is hoogst te betwijfelen dat in de werkelijke populatie van "Carnivoren" de frequentieverdeling van de groepen gelijk is als deze dataset. Bij het verzamelen van genetische data wordt namelijk eerder gekeken naar soorten van interesse (of soorten waar men gemakkelijk toegang toe krijgt). In het optimale geval zou de steekproef (database) de werkelijke frequentie van de populatie moeten benaderen. De steekproef is dus niet aselect gebeurd waardoor het resultaat vertekend zal zijn door de overrepresentatie van bepaalde dieren en de onderrepresentatie van andere dieren.

# Vraag 12: Nulhypothese ANOVA test in termen van de hellingsparameters (meerkeuzevraag)

Anova wordt uitgevoerd aan de hand van het lineair regressiemodel.
De nulhypothese kunnen we bijgevolg als volgt stellen:

1. \[ H_0: \beta_1 = \beta_2 = 0 \]
2. \[ H_0: \beta_1 = \beta_2  \]
3. \[ H_0: \beta_1 = \beta_2 \neq 0 \]
4. \[ H_0: \beta_1 \neq \beta_2 \]

# Vraag 13: Nulhypothese ANOVA test in termen van de hellingsparameters van het model

Met de ANOVA test testen we simultaan voor beide hellingsparameters. In woorden kan de nulhypothese op de twee hellingsparameters als volgt worden formuleerd: Het verschil in gemiddelde genoomgrootte tussen Carnivora en Primates alsook het verschil in gemiddelde genoomgrootte tussen Artiodactyla en Primates zijn beide gelijk aan nul.

Merk op dat de nulhypothese in termen van de hellingsparameters van het model uiteraard ook impliceert dat het verschil in gemiddelde genoomgrootte tussen Carnivora en Artiodactyla gelijk is aan nul.

# Vraag 14: Alternatieve hypothese ANOVA test
Anova wordt uitgevoerd aan de hand van het lineair regressiemodel.

1. $$H_A:$$ minstens één hellingsparameter is gelijk aan nul.
2. $$H_A:$$ alle hellingsparameters verschillen van nul.
3. $$H_A:$$ minstens één hellingsparameter verschilt van nul.
4. $$H_A:$$ alle hellingsparameters zijn gelijk aan nul.

# Vraag 15: Alternatieve hypothese ANOVA test in termen van de hellingsparameters

In woorden kan de alternatieve hypothese ook als volgt worden formuleerd: Minstens één van beide verschillen (het verschil in gemiddelde genoomgrootte tussen Carnivora en Primates en/of het verschil in gemiddelde genoomgrootte tussen Artiodactyla en Primates) is niet gelijk aan nul.

Merk op dat de alternatieve hypothese in termen van de hellingsparameters ook impliceert dat er een verschil kan zijn in gemiddelde genoomgrootte tussen Carnivora en Artiodactyla.

# Vraag 16: Voer de anovatest uit.

Sla het resultaat van de anovatest op in `anova_test`.

```{r}
anova_test <- anova(...)
```

# Vraag 17: Conclusie ANOVA (leesopdracht)

We besluiten dat de nulhypothese verworpen kan worden ($p = 4.86e-06$) en dat op het 5% significantieniveau de gemiddelde grootte van het genoom extreem significant verschilt tussen minstens twee van de drie bestudeerde groepen van mammalia: Primaten, Carnivoren en Artiodactyla.

# Vraag 18: Post-hoc analyse

Ga nu na door middel van een post-hoc analyse tussen welke groepen de genoomgrootte verschilt. Maak daarbij ook een betrouwbaarheidsinterval voor iedere paarsgewijze vergelijking.

```{r}
#Tukey test voor paarsgewijze vergelijking
library(multcomp)
mcp <- glht(...,linfct=mcp(... ="Tukey"))
summary(mcp)

#95% betrouwbaarheidsintervallen
set.seed(0)
conf <- confint(...)
par(mfrow=c(1,1), mar=c(3,10,2,1))  #pas margins aan binnen de figuur
plot(mcp)
```

# Vraag 19: Interpretatie post-hoc analyse (leesopdracht)
Interpretatie van de effectgroottes en de aangepaste p-waarden:

- De gemiddelde genoomgrootte (C-waarde) van Carnivora is 0,40 kleiner dan de gemiddelde genoomgrootte van Primates. Dit verschil is extreem significant op het globaal 5% significantieniveau (Tukey test, aangepaste p-waarde < 1e-04, 95% BI: [0.18; 0.61]).

- De gemiddelde genoomgrootte (C-waarde) van Artiodactyla is 0,32 kleiner dan de gemiddelde genoomgrootte van Primates. Dit verschil is extreem significant op het globaal 5% significantieniveau (Tukey test, aangepaste p-waarde = 0.0004, 95% BI: [0.13; 0.53]).

- De gemiddelde genoomgrootte (C-waarde) van Artiodactyla is 0,07 groter dan de gemiddelde genoomgrootte van Carnivora. Dit verschil is niet significant op het globaal 5% significantieniveau (Tukey test, aangepaste p-waarde = 0.81, 95% BI: [-0.18; 0.32]).

# Vraag 20: Conclusie (leesopdracht)

De gemiddelde grootte van het genoom verschilt extreem significant  tussen minstens twee van de drie bestudeerde groepen van mammalia ($p = 4.86e-06$).

Op een globaal 5% significantieniveau vinden we significante verschillen in gemiddelde genoomgrootte tussen de primaten enerzijds en de  carnivoren en evenhoevigen anderzijds.
De gemiddelde genoomgrootte van Primaten is respectievelijk 0.40 pg en 0.33 pg groter dan dat Carnivora en Artiodactyla Carnivora (Carnivora: Tukey test, aangepaste p-waarde < 1e-04, 95% BI: [0.18; 0.61] en Artiodactyla: Tukey test, aangepaste p-waarde = 0.0004, 95% BI: [0.13,0.53]).

De gemiddelde genoomgrootte is niet-significant verschillend tussen carnivoren en eenhoevigen op het 5% significantie-niveau (aangepaste p-waarde = 0.066).

Merk op: aangezien onze data niet onafhankelijk zijn, was onze analyse mogelijk te liberaal. Onze p-waarden zullen dus mogelijk te klein geschat zijn. De genomen zijn ook geen representatieve steekproef uit de populatie van zoogdieren waardoor de resultaten mogelijks vertekend zijn! 
