In protViz/SRMService: Report Quantitative Mass Spectrometry Data

knitr::opts_chunk$set(echo = FALSE, fig.width = 10, fig.height = 10, message = FALSE)

library(readr)
library(SRMService)
library(quantable)
library(reshape2)
library(plyr)
library(limma)
library(dplyr)
library(corrplot)

for(i in params$grp2$getConditions()$condition){
  quantable::mypairsSmooth(params$grp2$getNormalizedConditionData(i))
}

dd <- as.matrix(params$grp2$getNormalized()$data)
dd <- dd[complete.cases(dd),]
simpleheatmap(dd, margins = c(10,5))

M <- cor(params$grp2$getNormalized()$data, use="pairwise.complete.obs")
corrplot.mixed(M^2,upper="square", lower="pie",  diag="u", tl.pos="d")

Fit model

dir.create(params$OUTPUT_DIR)
intmat <- as.matrix(params$grp2$getNormalized()$data)
x<-quantable::matrix_to_tibble(intmat)
write.table(x , file = file.path( params$OUTPUT_DIR,"NormalizedIntensities.txt") )

fit <- lmFit(intmat , designMatrix)

Compute contrasts

cont <- limma::makeContrasts(contrasts = contrasts, levels =  levels)
lmfit.cont <- contrasts.fit(fit, cont)
lmfit.cont.ebayes <- eBayes(lmfit.cont)

res <- list()
for(i in 1:length(colnames(cont)))
{
  name <- colnames(cont)[i]
  print(name)
  res[[name]] <- data.frame(Condition = name, topTable(lmfit.cont.ebayes, coef=name, number=Inf))
}


res <- lapply(res,tibble::rownames_to_column,var="ProteinID")
res <- rbind.fill(res)
res$uprotID <- split2table(res$ProteinID,split="\\|")[,3]

write_tsv(res, path=file.path(params$OUTPUT_DIR,"limmaPValues.txt"))

cc <- data.frame(
  sl=c(0,0),
  p = c(0.01,0.05), 
  Area = c('p=0.01','p=0.05')
  )


res$Name <- split2table(res$ProteinID,split="\\|")[,3]
res$Name <- split2table(res$Name,split="\\_")[,1]

res$colour <- ifelse( grepl("REV_",res$ProteinID), "REV", "Forward")
quantable::multigroupVolcano(res, effect = "logFC",
                       type="adj.P.Val" ,
                       condition = "Condition", colour="colour",label="Name",
                       xintercept = c(-1.5,1.5))

Meta analysis of p-values and foldchanges

Foldchanges

logFCMap <- acast(res, ProteinID  ~ Condition, value.var = "logFC")
quantable::mypairsSmooth(logFCMap)

corrplot(cor(logFCMap, use="pairwise.complete.obs"),method="square", type="upper",  diag=TRUE, tl.pos="d")

PVals

adjPValMap <- acast(res, ProteinID  ~ Condition, value.var = "adj.P.Val", sep=".")
qvalue <- 0.05
nrSigSample <- adply(adjPValMap,2, function(x){ c( signifQ= sum(!is.na(x) & x <qvalue))},.id="Contrast")

p <- ggplot(nrSigSample, aes(x = Contrast, y= signifQ )) + geom_bar(stat="identity") 
p <- p + theme_classic()
p <- p + theme(axis.text.x = element_text(angle = 90, hjust = 1))
p

respVals <- merge(adjPValMap, logFCMap, by="row.names", suffixes = c(".adjPVal", ".log2FC"))
write_tsv(respVals, path=file.path(params$OUTPUT_DIR,"_pValFC_WideFormat.txt"))

nrSigProtein <- adply(adjPValMap,1, function(x){c("nrSigComparisons" = sum( !is.na(x) & x  < qvalue))}, .id="proteinID")
x <- table(nrSigProtein$nrSigComparisons[nrSigProtein$nrSigComparisons > 0])

x <- res %>% dplyr::group_by_at( "ProteinID" ) %>% dplyr::summarise(n=n(),nrSigComparisons = sum(!is.na(adj.P.Val) & adj.P.Val <qvalue)) 

x <- x %>% dplyr::filter( nrSigComparisons > 0) %>% group_by(nrSigComparisons) %>% dplyr::summarize(nr = n())
class(x$nrSigComparisons) <- "character"

p <- ggplot(x, aes(x = nrSigComparisons, y= nr )) + geom_bar(stat="identity") 
p <- p + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))
p

\newpage

Table of Significant Fold changes

colnames(res)
res2 <- res
res2$P.Value <- NULL
res2$ProteinID <-NULL
res2$AveExpr <- NULL
res2$t <- NULL
res2$B <- NULL
knitr::kable(subset(res2,adj.P.Val <0.05))