---
title: "Exercice 4"
output: html_document
date: "2023-06-23"
---


<script>
function myFunction(ID) {
    var x = document.getElementById(ID);
    if (x.style.display === "none") {
        x.style.display = "block";
    } else {
        x.style.display = "none";
    }
}
</script>


```{r,echo=F,warning=F,message=F}
if (!require("scales")) install.packages("scales")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("gplots")) install.packages("gplots")
```

# Load data

- First, load the results from Exerice 2:

```{r,warning=F}
load("data/Ex2.Results.Rdata")
```


# Difference in expression

- For this part, you only care about Control condition. **SD condition is filtered out !**.
  - Filter accordingly metadata and expression data

- Using boxplot (basic or ggplot2), show *Spcs1* gene expression: boxplots must be separated by Strain and Time

- Assuming this gene expression follows a normal distribution, you want to write a statistical. 
  - What would happend if you use Time as numerical or a factor ? (here Time is based on 24h cycle, with light-dark transition at Time 12)
  - Add a new variable in your metadata representing Time as factor (named Time_as_Factor).

- you write to linear models to fit *Spcs1* expression:
 -  ~ Strain * Time_as_Factor
 -  ~ Strain + Strain:Time_as_Factor
 - What is the difference between these two models ? What the coefficients represent ?

<button onclick="myFunction(ID='myDIV2')" style="background-color: red"><span style='color:white'>Show answer !</span></button>

<div id="myDIV2" style="display:none">

<br>

First, we filter our data to only keep control condition

```{r}
idx<-metadata$Condition == "Ctr"
metadata_ctr<-metadata[idx,]
expression_ctr<-expression[,idx]

# We refactor Strain, to remove "d2" and "NA" levels
metadata_ctr$Strain<-factor(metadata_ctr$Strain)
```

We use boxplot to plot gene expression for each Strain and Time

```{r}
boxplot(expression_ctr["Spcs1",]~metadata_ctr$Time + factor(metadata_ctr$Strain),xlab="Time.Strain",ylab="Spcs1 expression")
```

Using ggplot2 you need your x value to be factor, and fill by Strain. 

```{r}
MyData<-cbind.data.frame(metadata_ctr,Expression=expression_ctr["Spcs1",])
bp<-ggplot(aes(x=factor(Time),y=Expression,fill=Strain),data=MyData) + geom_boxplot() + 
  scale_x_discrete(breaks=c(0,6,12,18),labels=c(0,6,12,18))
bp
```

Or you can using *facet_grid()* function to separate by Strain

```{r}
MyData<-cbind.data.frame(metadata_ctr,Expression=expression_ctr["Spcs1",])
bp<-ggplot(aes(Time,Expression,group=Time,fill=Strain),data=MyData) + geom_boxplot()
bp + facet_grid(. ~ Strain) + scale_x_continuous(breaks=c(0,6,12,18),labels=c(0,6,12,18))
```

If we use Time as a numerical value, we consider that gene expression will increase or decrease linearly with time (linear regression). This is often not the case on 24h cycles, where expression will be cyclic. Transforming Time into factors allows to have any pattern of expression.

We add a new variable as Time_as_Factor

```{r}
metadata_ctr$Time_as_Factor<-factor(metadata_ctr$Time)
```

We write two test, model1 and model2

```{r}
model1<-lm(expression_ctr["Spcs1",]~ Time_as_Factor * Strain,data=metadata_ctr)
model2<-lm(expression_ctr["Spcs1",]~ Time_as_Factor + Strain:Time_as_Factor,data=metadata_ctr)
```

See the summary of these two models:

```{r}
summary(model1)
```


```{r}
summary(model2)
```

The fitted value of these two models are exactly identical.

```{r}
fitted(model1) == fitted(model2)
```

And we see that the first 5 coefficients are identical, only the last 3 coefficients change.

To model the expression for Strain 'D2' at time 6-12-18, the first model use the following coefficients (excluding 0)

model1 (D2-Time 6): $expression = \beta_0 + \beta_1 + \beta4 + \beta5$

```{r}
sum(coef(model1)[c(1,2,5,6)])
mean(expression_ctr["Spcs1", metadata_ctr$Strain == "D2" & metadata_ctr$Time_as_Factor == "6"])
```

while model 2 uses: $expression = \beta_0 + \beta_1 + \beta5$

```{r}
sum(coef(model2)[c(1,2,6)])
```

In model1, the last 3 coefficients represent the difference of expression between D2 and B6 at time 6,12,18 compared to the difference at time 0
In model2, the last 3 coefficients represent the difference of expression between D2 and B6 at time 6,12,18, independently of the difference at time 0


In model1, the $\beta_4$ coefficient (`r names(coef(model1))[5]`) is included in all time-points of D2 expression

See. model matrix

```{r}
model.matrix(model1)
```


```{r}
model.matrix(model2)
```

</div>


# Heatmap

- You want to show a heatmap of the expression data.
- Using gplots show expression
- Add red-green color, scale data (Z-score transformation) and show Mouse Strain on the top of the heatmap 

<button onclick="myFunction(ID='myDIV3')" style="background-color: red"><span style='color:white'>Show answer !</span></button>

<div id="myDIV3" style="display:none">

<br>


Before using heatmap.2, you want to scale the data:

```{r}
expression_scaled<-t(apply(expression,1,scale))
```

You can then plot your data:

```{r,warning=F}
library(gplots)
heatmap.2(expression_scaled,col=redgreen(75),trace="none",Colv = T,
          ColSideColors = c("red","blue")[factor(metadata$Strain)] )
```

If you use scale="row" instead, heatmap.2 will use hclust on non-scaled data, then represent data as scaled

```{r}
heatmap.2(expression,col=redgreen(75),trace="none",Colv = T,scale="row",
          ColSideColors = c("red","blue")[factor(metadata$Strain)] )
```


</div>