inst/doc/NPARC.R

In NPARC: Non-parametric analysis of response curves for thermal proteome profiling experiments

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----dependencies, message=FALSE----------------------------------------------
library(dplyr)
library(magrittr)
library(ggplot2)
library(broom)
library(knitr)
library(NPARC)

## ----load_data----------------------------------------------------------------
data("stauro_TPP_data_tidy")

## -----------------------------------------------------------------------------
df <- stauro_TPP_data_tidy

## ----data_head----------------------------------------------------------------
df %>% 
  mutate(compoundConcentration = factor(compoundConcentration), 
         replicate = factor(replicate), 
         dataset = factor(dataset)) %>% 
  summary()


## ----filter_qupm_or_NAs, results='asis'---------------------------------------
df %<>% filter(uniquePeptideMatches >= 1)
df %<>% filter(!is.na(relAbundance))

## ----count_curves_per_protein-------------------------------------------------
# Count full curves per protein
df %<>%
  group_by(dataset, uniqueID) %>%
  mutate(n = n()) %>%
  group_by(dataset) %>%
  mutate(max_n = max(n)) %>% 
  ungroup

table(distinct(df, uniqueID, n)$n)

## -----------------------------------------------------------------------------
# Filter for full curves per protein:
df %<>% 
  filter(n == max_n) %>%
  dplyr::select(-n, -max_n)

## ----select_stk4--------------------------------------------------------------
stk4 <- filter(df, uniqueID == "STK4_IPI00011488")

## -----------------------------------------------------------------------------
stk4 %>% filter(compoundConcentration == 20, replicate == 1) %>% 
  dplyr::select(-dataset) %>% kable(digits = 2)

## ----plot_stk4----------------------------------------------------------------
stk4_plot_orig <- ggplot(stk4, aes(x = temperature, y = relAbundance)) +
  geom_point(aes(shape = factor(replicate), color = factor(compoundConcentration)), size = 2) +
  theme_bw() +
  ggtitle("STK4") +
  scale_color_manual("staurosporine (mu M)", values = c("#808080", "#da7f2d")) +
  scale_shape_manual("replicate", values = c(19, 17))

print(stk4_plot_orig)

## ----fit_null_stk4------------------------------------------------------------
nullFit <- NPARC:::fitSingleSigmoid(x = stk4$temperature, y = stk4$relAbundance)

## ----summarize_null_stk4------------------------------------------------------
summary(nullFit)

## ----augment_null_stk4--------------------------------------------------------
nullPredictions <- broom::augment(nullFit)

## ----head_augment_result------------------------------------------------------
nullPredictions %>% filter(x %in% c(46, 49)) %>% kable()

## ----add_null_resids_stk4-----------------------------------------------------
stk4$nullPrediction <- nullPredictions$.fitted
stk4$nullResiduals <- nullPredictions$.resid

stk4_plot <- stk4_plot_orig + geom_line(data = stk4, aes(y = nullPrediction))

print(stk4_plot)

## ----fit_alternative_stk4-----------------------------------------------------
alternativePredictions <- stk4 %>%
# Fit separate curves per treatment group:
  group_by(compoundConcentration) %>%
  do({
    fit = NPARC:::fitSingleSigmoid(x = .$temperature, y = .$relAbundance, start=c(Pl = 0, a = 550, b = 10))
    broom::augment(fit)
  }) %>%
  ungroup %>%
  # Rename columns for merge to data frame:
  dplyr::rename(alternativePrediction = .fitted,
                alternativeResiduals = .resid,
                temperature = x,
                relAbundance = y)

## ----add_alternative_resids_stk4----------------------------------------------
stk4 <- stk4 %>%
  left_join(alternativePredictions, 
            by = c("relAbundance", "temperature", 
                   "compoundConcentration")) %>%
  distinct()

## ----plot_Fig_2_A_B-----------------------------------------------------------
stk4_plot <- stk4_plot +
  geom_line(data = distinct(stk4, temperature, compoundConcentration, alternativePrediction), 
            aes(y = alternativePrediction, color = factor(compoundConcentration)))

print(stk4_plot)

## ----compute_rss_stk4---------------------------------------------------------
rssPerModel <- stk4 %>%
  summarise(rssNull = sum(nullResiduals^2),
            rssAlternative = sum(alternativeResiduals^2))

kable(rssPerModel, digits = 4)

## -----------------------------------------------------------------------------
BPPARAM <- BiocParallel::SerialParam(progressbar = FALSE)

## ----fit_all_proteins---------------------------------------------------------
fits <- NPARCfit(x = df$temperature, 
                 y = df$relAbundance, 
                 id = df$uniqueID, 
                 groupsNull = NULL, 
                 groupsAlt = df$compoundConcentration, 
                 BPPARAM = BPPARAM,
                 returnModels = FALSE)

str(fits, 1)

## ----show_fit_results---------------------------------------------------------
fits$metrics %>% 
  mutate(modelType = factor(modelType), nCoeffs = factor(nCoeffs), nFitted = factor(nFitted), group = factor((group))) %>% 
  summary

fits$predictions %>% 
  mutate(modelType = factor(modelType), group = factor((group))) %>% 
  summary

## ----show_stk4_results--------------------------------------------------------
stk4Metrics <- filter(fits$metrics, id == "STK4_IPI00011488")

rssNull <- filter(stk4Metrics, modelType == "null")$rss
rssAlt <- sum(filter(stk4Metrics, modelType == "alternative")$rss) # Summarize over both experimental groups

rssNull
rssAlt

## ----plot_stk4_results--------------------------------------------------------
stk4Predictions <- filter(fits$predictions, modelType == "alternative", id == "STK4_IPI00011488")

stk4_plot_orig +
  geom_line(data = filter(stk4Predictions, modelType == "alternative"), 
            aes(x = x, y = .fitted, color = factor(group))) +
    geom_line(data = filter(stk4Predictions, modelType == "null"), 
            aes(x = x, y = .fitted))

## ----testDOFiid---------------------------------------------------------------
modelMetrics <- fits$metrics 
fStats <- NPARCtest(modelMetrics, dfType = "theoretical")

## ----plotF_DOFiid-------------------------------------------------------------
fStats %>% 
  filter(!is.na(pAdj)) %>%
  distinct(nFittedNull, nFittedAlt, nCoeffsNull, nCoeffsAlt, df1, df2) %>%
  kable()

## ----plotFiid-----------------------------------------------------------------
ggplot(filter(fStats, !is.na(pAdj))) +
  geom_density(aes(x = fStat), fill = "steelblue", alpha = 0.5) +
  geom_line(aes(x = fStat, y = df(fStat, df1 = df1, df2 = df2)), color = "darkred", size = 1.5) +
  theme_bw() +
  # Zoom in to small values to increase resolution for the proteins under H0:
  xlim(c(0, 10))

## ----plotPiid-----------------------------------------------------------------
ggplot(filter(fStats, !is.na(pAdj))) +
  geom_histogram(aes(x = pVal, y = ..density..), fill = "steelblue", alpha = 0.5, boundary = 0, bins = 30) +
  geom_line(aes(x = pVal, y = dunif(pVal)), color = "darkred", size = 1.5) +
  theme_bw()

## ----testDOFemp---------------------------------------------------------------
modelMetrics <- fits$metrics 
fStats <- NPARCtest(modelMetrics, dfType = "empirical")

## ----plotFnew-----------------------------------------------------------------
ggplot(filter(fStats, !is.na(pAdj))) +
  geom_density(aes(x = fStat), fill = "steelblue", alpha = 0.5) +
  geom_line(aes(x = fStat, y = df(fStat, df1 = df1, df2 = df2)), color = "darkred", size = 1.5) +
  theme_bw() +
  # Zoom in to small values to increase resolution for the proteins under H0:
  xlim(c(0, 10))

## ----plotPnew-----------------------------------------------------------------
ggplot(filter(fStats, !is.na(pAdj))) +
  geom_histogram(aes(x = pVal, y = ..density..), fill = "steelblue", alpha = 0.5, boundary = 0, bins = 30) +
  geom_line(aes(x = pVal, y = dunif(pVal)), color = "darkred", size = 1.5) +
  theme_bw()

## ----selectHits---------------------------------------------------------------
topHits <- fStats %>% 
  filter(pAdj <= 0.01) %>%
  dplyr::select(id, fStat, pVal, pAdj) %>%
  arrange(-fStat)

## ----showHits, results='asis'-------------------------------------------------
knitr::kable(topHits)

## ----session------------------------------------------------------------------
devtools::session_info()