---
title: "FAXtraction analysis"
output: html_document
---

## R-packages
```{r, eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
library(plyr)
library(tidyverse)
library(broom)
library(scales)

predictvals <- function(model, xvar, yvar, xrange=NULL, samples=100, ...) {
  if (is.null(xrange)) {
    if (any(class(model) %in% c("lm", "glm")))
      xrange <- range(model$model[[xvar]])
    else if (any(class(model) %in% "loess"))
      xrange <- range(model$x)
  }
  newdata <- data.frame(x = seq(xrange[1], xrange[2], length.out = samples))
  names(newdata) <- xvar
  newdata[[yvar]] <- predict(model, newdata = newdata, ...)
  newdata
}
```

## Import data
```{r, echo=FALSE, message=FALSE, warning=FALSE}
d <- read_delim(file.choose(), "\t")
```

## Format data
```{r, eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
df <- d %>%
  select(id, `Fasta headers`, starts_with("LFQ")) %>%
  gather("h", "value", 3:ncol(.)) %>%
  separate(h, c("fac", "rep"), "_") %>%
  mutate(fac = as.integer(gsub("LFQ intensity ", "", fac)), rep = gsub("FA", "", rep)) %>%
  group_by(id, rep) %>%
  mutate(norm = value/max(value)) %>%
  filter(!is.na(norm))

meta <- df %>%
  group_by(id) %>%
  select(1, fasta = 2, 5) %>%
  mutate(abundance = max(value)) %>%
  select(-value) %>%
  ungroup() %>%
  distinct(id, .keep_all = T)
```

## Modeling
```{r, eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
make_model <- function(data) {glm(norm ~ fac, data, family=binomial)}

predvals <- dlply(df, .(id, rep), make_model) %>%
  ldply(.fun=predictvals, xvar="fac", yvar="norm", type="response")

m <- df %>%
  group_by(id, rep) %>%
  do(exp = make_model(.))

coef <- tidy(m, exp) %>%
  select(1:4) %>%
  spread("term", "estimate") %>%
  select(1, 2, a = 3, b = 4) %>%
  mutate(f50 = -a/b, slope = (b * exp(a + b * f50 )/(exp(a + b * f50)+1)^2), test = f50 > 60 & f50<100 & slope > 0.025)

hit <- coef %>%
  group_by(id) %>%
  summarise(hit = sum(test)) %>%
  left_join(meta, by = "id") %>%
  filter(hit > 0, abundance > 0) %>%
  select(1, 2, 4, 3)
```

## Subset data based on hits
```{r, eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
dfhit <- merge(x=hit, y=df, by = "id")
predvalshit <- merge(x=hit, y=predvals, by="id")
```

## Plot data

```{r, eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
p <- ggplot(data = df, mapping = aes(x = fac, y = norm, color=as.factor(rep))) +
  geom_line(data=predvals, alpha=1, size=1) +
  geom_point(size=2) +
  scale_x_continuous(breaks=pretty_breaks(n=6)) +
  theme_light() +
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5)) +
  ylab("Normalized abundance (AU)") + xlab("Formic acid (%)") +labs(color="Replicate") +
  facet_wrap(~id, ncol=25)

phit <- ggplot(data = dfhit, mapping = aes(x = fac, y = norm, color=as.factor(rep))) +
  geom_line(data=predvalshit, alpha=1, size=1) +
  geom_point(size=2) +
  scale_x_continuous(breaks=pretty_breaks(n=6)) +
  theme_light() +
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5)) +
  ylab("Normalized abundance (AU)") + xlab("Formic acid (%)") +labs(color="Replicate") +
  facet_wrap(~id, ncol=25)
```

## Export data
```{r, eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
write_delim(hit, "hits.txt", delim = "\t")
ggsave(filename="plot_all.pdf", plot=p, width=25, height=nrow(meta)/25, useDingbats=FALSE, limitsize=FALSE)
ggsave(filename="plot_hits.pdf", plot=phit, width=ifelse(nrow(hit)<25, nrow(hit), 25), height=nrow(hit)/25+1.5, useDingbats=FALSE, limitsize=FALSE)
```




