real_data.Rmd

---
title: "macrophage"
author: "zhu"
date: "20/12/2018"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F,message = F,warning = F,results = F)
```

```{r data transformation}
library(tidyverse)
library(readxl)
library(knitr)

colnames <- colnames(read_xlsx('macrophage%20polarization.xlsx',range = 'PeakArea!I1:BL1'))
df <- read_xlsx('macrophage%20polarization.xlsx',range = 'PeakArea!I4:BL43',col_names = colnames)
#df.raw <- read_xlsx('macrophage%20polarization.xlsx',range = 'PeakArea!I4:BL43',col_names = colnames)
#rm(colnames)

i=2
j=5
r=4
time=c(0,2,6,10,24)

df$togroup <- rep(1:10,each=4)
df[20,] <- df[17:19,] %>% 
  summarise_all(.,funs(mean))
df.clean <- df %>% 
  group_by(togroup) %>% 
  mutate_all(~ifelse(is.na(.x), mean(.x, na.rm = TRUE), .x))#fill NA with mean of replicates. 

df.clean[1:4,"Fructose-1,6-bisphosphate"] <- df.clean$`Fructose-1,6-bisphosphate`[5]#F16B is missing values at entire time 1 control. filled with values at time 1 treatment. the assumption is that value at time 0 for both control and treatment is the same

plot.mean_sd <- function(df,...){
  mean_sd <- df %>% 
    group_by(togroup) %>% 
    summarise_all(.,funs(mean,sd)) %>% 
    select(.,-togroup) %>% 
    flatten()
  
  plot(mean_sd[1:560],mean_sd[561:1120],type = 'p',xlab = 'mean',ylab = 'sd',...)
  
}
df.log <- log(df.clean)#log transform

plot.mean_sd(df.log, main = 'log transformed', sub = 'log transformed shows homoscedastic')
plot.mean_sd(df.clean, main = 'original data', sub = 'heteroscedastic')#check the homoscedasticity

mean_sd.df <- df.log %>% 
    group_by(togroup) %>% 
    summarise_all(.,funs(mean,sd))

df.clean$togroup <- NULL
df.clean$Leucine <- NULL#leucine has lots of NA. so removed
df.log$Leucine <- NULL
df.log$togroup <- NULL

df.final.masigpro <- t(df.log)#for masigpro
df.final <- t(scale(df.log))#for asca
colnames <- rownames(df.final)

```




1. original mean vs sd shows linear dependency. log transform corrected it therefore it is needed.
2. MaSigpro is a univariate method, differences in sd between features has no effects. mean centering is also not needed. Therefore, to perform Masigpro with only log transformed data.
3. ASCA-gene however is multivariate therefore need both mean centering and scaling of the features.

```{r ssq variances}
calculate.ss <- function(tag,df.clean){
  mean.grand <- df.clean %>% 
    ungroup() %>% 
    mutate_all(funs(mean))
  df.clean$tag <- tag
  mean.sub <- df.clean %>%
    group_by(tag) %>%
    mutate_all(funs(mean)) %>% 
    ungroup() %>% 
    select(.,-tag)
  (mean.sub-mean.grand)^2
  # mean.grand
  # mean.sub
  
}
df.clean <- df.log
mean.grand <- df.clean %>% 
  ungroup() %>% 
  mutate_all(funs(mean))
ss.treatment <- calculate.ss(rep(1:2,each=20),df.clean)
ss.time <- calculate.ss(rep(rep(1:5,each=4),2),df.clean)
ss.replicate <- df.clean %>% 
  mutate_all(.,funs((.-mean(.))^2))
total.sum1 <- calculate.ss(1:40,df.clean)
total.sum <- sum(ss.treatment,ss.time,ss.replicate)
total.sum1 <- sum(total.sum1)
percent_var.treatment <- sum(ss.treatment)/total.sum
percent_var.time <- sum(ss.time)/total.sum
percent_var.rep <- sum(ss.replicate)/total.sum

df.var.realdata <- data.frame(sourc=c('treatment','time','rep'),ss.percentage=c(percent_var.treatment*100,percent_var.time*100,percent_var.rep*100))
knitr::kable(df.var.realdata)
```


```{r function plot metabolites}
# plot_metabolites <- function(df,range,...){
#   df.final <- df
#   df.toplot <- data.frame(t(df.final[range,]))
#   df.toplot$time <- rep(rep(time,each=r),i)
#   df.toplot$treatment <- rep(1:2,each=r*j)
#   a <- df.toplot %>% 
#     gather(key = metabolites,value = value,1:length(range))
#   
#   output <- ggplot(a,aes(x=time,y=value,color=factor(treatment)))+
#     geom_point(size=0.5)+
#     stat_summary(fun.y = mean,geom = 'line')+
#     theme(legend.title = element_text(size = 5))+
#     facet_wrap(metabolites~., scales = 'free')+
#     geom_smooth(method = 'lm', formula = y~poly(x,2),se = F,linetype = '3313')+
#     labs(...)
#   output
# }
# plot_a_metabolite <- function(df,FUN,which){
#   df.final <- df
#   trend.toplot <- data.frame(replicate=rep(1:(i*j),each=r),time=rep(time,each=r),treatment=factor(rep(1:i,each=j*r)),metabolite=df.final[which,])
#   ggplot(trend.toplot,aes(x=time,y=metabolite,color=treatment))+
#     geom_point()+
#     stat_summary(fun.y = FUN,geom = 'line')+
#     geom_smooth(method = 'lm', formula = y~poly(x,2),se = F,linetype = '3313')
# }
```

```{r plot 1:27}
source('utils.R')
ascaonly <- plot_metabolites(df.final.masigpro,c('CMP','AMP','ATP'),formula = y~x,tag = 'A')
masigproonly <- plot_metabolites(df.final.masigpro,c("Glycerol-phosphate","Isocitric acid","Lactic Acid"),tag = 'B')
both <- plot_metabolites(df.final.masigpro,c("CTP","GSH","IMP"),tag = 'B')
# plot_metabolites(df.final.masigpro,37:51)
plot(ascaonly)
both
```

```{r}
# plot_a_metabolite(df.final.masigpro,mean,"ADP")
# plot_a_metabolite(df.final.masigpro,median,"Fructose-1,6-bisphosphate")
```

```{r,results=T}

masigpro.design <- design_matrix(i,j,r,time)
degree <- 2
design <- make.design.matrix(edesign = masigpro.design,degree = degree)
dis <- design$dis
edesign <- design$edesign
colnames(df.final.masigpro) <- rownames(edesign)

masigpro.fit <- maSigPro(df.final.masigpro,masigpro.design,degree=degree, rsq = 0.4)
kable(masigpro.fit$summary)
toplot <- rownames(masigpro.fit$sig.genes$treatment2vstreatment1$sig.profiles)
```

```{r,include=F}
# AA <- df.final.masigpro["GSSG",]
# PlotGroups(AA,edesign = edesign, show.fit = T, dis = dis, groups.vector = design$groups.vector,main = 'GSSG')
```

```{r}
mx <- asca.design.matrix(i,j,r,time)

Fac=c(1,1,2,2)
type=1
source('ASCA-genes.1.2.1/sourceASCA.R')
asca.fit <- ASCA.2f_leverage(X = t(df.final),Designa = mx$j, Designb = mx$i,type = type, Fac = Fac)
```

```{r}
source('RMSECV.R')
rmsecv <- RMSECV(asca.fit$Model.bab$data,8,4)
rmsecv
```


```{r,results=T}
permut_wrapped_0.05 <- wrap.permutation.info(df.final = df.final,asca.fit = asca.fit,groups = colnames,R = 100,which_leverage = ASCA.2f_leverage)
permut_wrapped_0.01 <- wrap.permutation.info(df.final = df.final,asca.fit = asca.fit,groups = colnames,R = 100,which_leverage = ASCA.2f_leverage,alpha=0.01)

Nulldist <- permut_wrapped_0.05$Null_distribution_plot+
  geom_hline(yintercept = permut_wrapped_0.01$lev_limit,color='red')
#ggsave('Nulldist.png',path = '../../Google Drive/BDA_internship/thesis/Mythesis/data/')
ggsave('Nulldist.png',path = 'D:/google drive/BDA_internship/thesis/Mythesis/data/')
# Nulldist_0.01 <- permut_wrapped_0.01$Null_distribution_plot
# #print(sort(masigpro.fit$summary[,2]))
# # venn_list <- list(maSigpro = sort(masigpro.fit$summary[,2]), ASCA.gene = sort(colnames[asca.fit$Model.bab$leverage>permut_wrapped$lev_limit]))
# # venn::venn(venn_list)
# Nulldist_0.05
```

in asca-gene but not masigpro

```{r,include=F}
AA <- df.final.masigpro["ATP",]
PlotGroups(AA,edesign = edesign, show.fit = T, dis = dis, groups.vector = design$groups.vector,main = 'Glycerol-phosphate')
```

```{r,include=F}
AA <- df.final.masigpro["Isocitric acid",]
PlotGroups(AA,edesign = edesign, show.fit = T, dis = dis, groups.vector = design$groups.vector,main = 'Isocitric acid')
```

```{r,include=F}
AA <- df.final.masigpro["Lactic Acid",]
PlotGroups(AA,edesign = edesign, show.fit = T, dis = dis, groups.vector = design$groups.vector,main = 'Lactic Acid')
```

```{r,include=F}
AA <- df.final.masigpro["Cytidine",]
PlotGroups(AA,edesign = edesign, show.fit = T, dis = dis, groups.vector = design$groups.vector,main = 'Cytidine')
```

in both

```{r,include=F}
AA <- df.final.masigpro["Phosphoglyceric acid",]
PlotGroups(AA,edesign = edesign, show.fit = T, dis = dis, groups.vector = design$groups.vector,main = 'Tryptophan')
```

```{r,include=F}
AA <- df.final.masigpro["Ascorbic Acid",]
PlotGroups(AA,edesign = edesign, show.fit = T, dis = dis, groups.vector = design$groups.vector,main = 'Ascorbic Acid')
```