2_proteomics.Rmd

---
title: "Proteomics: data analysis"
date: "`r Sys.Date()`"
author:
  - name: Arun Seetharam
    affiliation: Tuteja Lab
    affiliation_url: https://www.tutejalab.org
output:
  #pdf_document: default
  rmdformats::readthedown:
    self_contained: true
    thumbnails: false
    lightbox: true
    gallery: true
    highlight: tango
---

```{r setup, include=FALSE}
options(max.print = "125")
knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>",
  fig.path = "assets/",
  fig.width = 8,
  prompt = FALSE,
  tidy = FALSE,
  message = FALSE,
  warning = TRUE
)
knitr::opts_knit$set(width = 75)
```

# Prerequisites 

Load the packages needed for the analyses and read-in the relavant data

```{r libload, warnings=TRUE, message=FALSE}
library(scales)
library(data.table)
library(vctrs)
library(tidyverse)
library(plotly)
library(TissueEnrich)
library(ggrepel)
library(RColorBrewer)
library(pheatmap)
setwd("/work/LAS/geetu-lab/arnstrm/mouse.trophoblast.smallRNAseq")
dat <-
  read.csv("assets/data_processed.csv",
           row.names = NULL,
           stringsAsFactors = TRUE)
orig_data <-
  read.csv("assets/orig_data.csv",
           row.names = NULL,
           stringsAsFactors = TRUE)
source("assets/theme_clean.R")
```

quick check:

```{r qcdata, warnings=TRUE, message=FALSE}
dat_tidy <- dat %>% drop_na()
dat %>% summarise(Number_of_proteins = n())
```
# Volcano plot

```{r pvalPlot, fig.cap="Figure 1: Volcano plot (proteomics)", fig.width=8, fig.height=6}
df <- orig_data
exp <- log(1.5, 2)
df$log2fc <- log(df$Ratio, base = 2)
df$negLog10p = -log10(df$Qvalue)
df$diffexpressed[df$log2fc <= -exp &
                   df$Qvalue <= 0.05] <- "up in pTGC"
df$diffexpressed[df$log2fc >= exp &
                   df$Qvalue <= 0.05] <- "up in TSC"
df$diffexpressed[df$log2fc >= -exp &
                   df$log2fc <= exp] <- "other genes"
df$diffexpressed[df$Qvalue > 0.05] <- "other genes"
df$newGenes <-
  str_replace_all(string = df$Genes,
                  pattern = "\\;.*$",
                  replacement = "")
data.table::setDT(df)[diffexpressed == "other genes", newGenes := NA]

g <-
  ggplot(data = df,
         aes(
           x = log2fc,
           y = negLog10p,
           col = diffexpressed,
           label = newGenes
         )) +
  geom_point(alpha = 0.5) +
  theme_classic() +
  scale_color_manual(name = "Expression",
                     values = c("grey", "#c6007b", "#a0b600")) +
  ggtitle(paste("pTGC vs. TSC (proteomics)")) +
  #  geom_text_repel(show.legend = F) +
  xlab("Log2 fold change") +
  ylab("-log10 pvalue") +
  theme(legend.text.align = 0)
ggplotly(g)
#ggsave("prot_volc.png", dpi=900, width = 8, height = 6)
```

# PCA plot

This was run on another R version due to incompatibility with existing packages.
The command is shown below, but it is not executed when knitting the document.

```{r pcaAnalysis, eval = FALSE}
library("MetaboAnalystR")
library(ggplot2)
library(ggrepel)

setwd("~/metabo")
mSet <- InitDataObjects("conc", "stat", FALSE)
mSet <- Read.TextData(mSet, "data_processed.csv", "colu", "disc")

mSet <- SanityCheckData(mSet)
mSet <- RemoveMissingPercent(mSet, percent = 0.5)
mSet <- ImputeMissingVar(mSet, method = "min")
mSet <- PreparePrenormData(mSet)
mSet <-
  Normalization(mSet,
                "NULL",
                "LogNorm",
                "NULL",
                ratio = FALSE,
                ratioNum = 20)
mSet <-
  PlotSampleNormSummary(mSet, "snorm_0_", "png", 300, width = NA)
mSet <- PCA.Anal(mSet)
saveRDS(mSet, "mSet.rds")
```

Read in the saved RDS file from the previous step and plot the PCA for samples.

```{r pcaPlot, fig.cap="Figure 2: PCA plot (proteomics)", fig.width=10, fig.height=6}
mSet <- readRDS("assets/mset.rds")
group <- sapply(strsplit(rownames(mSet$analSet$pca$x), "_"), "[", 1)
intgroup.df <- as.data.frame(group)
d <- data.frame(
  PC1 = mSet$analSet$pca$x[, 1],
  PC2 = mSet$analSet$pca$x[, 2],
  intgroup.df,
  name = rownames(mSet$analSet$pca$x)
)
# rename
d$sample <- "NA"
d$rep <- gsub("$", ")", gsub("^rep", "(rep", str_extract(d$name,"rep.")))
d$sample[which(str_detect(rownames(d), "Diff_rep"))] <- "pTGC"
d$sample[which(str_detect(rownames(d), "Undiff_rep"))] <- "TSC"
d$group <- d$sample
d$name <- paste0(d$sample, d$rep, sep = " ")
d <- d %>% dplyr::select(PC1:name)
g <- ggplot(d, aes(PC1, PC2, color = group)) +
  scale_shape_manual(values = 1:10) +
  scale_color_manual(values = c('pTGC'      = '#c6007b',
                                'TSC'  = '#a0b600')) +
  theme_classic() +
  theme(legend.title = element_blank()) +
  geom_point(size = 2, stroke = 2) +
  geom_text_repel(aes(label = name)) +
  xlab(paste("PC1", round(mSet$analSet$pca$variance[1] * 100, 2), "% variance")) +
  ylab(paste("PC2", round(mSet$analSet$pca$variance[2] * 100, 2), "% variance"))
g
```
# Tissue Enrichment 

```{r tissueEnrich, warning=TRUE, message=TRUE}
plotTE <- function(inputGenes = gene.list,
                   myColor = "color") {
  gs <-
    GeneSet(geneIds = inputGenes,
            organism = "Mus Musculus",
            geneIdType = SymbolIdentifier())
  output <- teEnrichment(inputGenes = gs, rnaSeqDataset = 3)
  en.output <-
    setNames(data.frame(assay(output[[1]]),
                        row.names = rowData(output[[1]])[, 1]),
             colData(output[[1]])[, 1])
  en.output$Tissue <- rownames(en.output)
  logp <- -log10(0.05)
  en.output <-
    mutate(en.output,
           significance = ifelse(Log10PValue > logp,
                                 "colored", "nocolor"))
  en.output$Sig <- "NA"
  ggplot(en.output, aes(reorder(Tissue, Log10PValue),
                        Log10PValue,
                        fill = significance)) +
    geom_bar(stat = 'identity') +
    theme_clean() + ylab("- log10 adj. p-value") + xlab("") +
    scale_fill_manual(values = c("colored" = myColor, "nocolor" = "gray")) +
    scale_y_continuous(expand = expansion(mult = c(0, .1)),
                       breaks = scales::pretty_breaks()) +
    coord_flip()
}
exp <- log(1.2, 2)
diff <- df$newGenes[df$log2fc <= -exp & df$Qvalue <= 0.05]
undiff <- df$newGenes[df$log2fc >= exp & df$Qvalue <= 0.05]
df.diff <- dat_tidy %>%
  rowwise() %>%
  mutate(Diff = mean(c(
    Diff_rep1,   Diff_rep2,  Diff_rep3,  Diff_rep4,  Diff_rep5
  ))) %>%
  dplyr::select(proteinID, Diff) %>%
  dplyr::filter(Diff > 0)  %>%
  ungroup() %>%
  mutate(quart = ntile(Diff, 4)) %>%
  mutate(decile = ntile(Diff, 10))
df.undiff <- dat_tidy %>%
  rowwise() %>%
  mutate(Undiff = mean(
    c(
      Undiff_rep1,
      Undiff_rep2,
      Undiff_rep3,
      Undiff_rep4,
      Undiff_rep5
    )
  )) %>%
  dplyr::select(proteinID, Undiff) %>%
  dplyr::filter(Undiff > 0)  %>%
  ungroup() %>%
  mutate(quart = ntile(Undiff, 4)) %>%
  mutate(decile = ntile(Undiff, 10))

filterCuts <- function(dataIn = df.undiff,
                       cutOff = 4,
                       type = decile) {
  type <- enquo(type)
  dataIn %>%
    dplyr::filter(!!type == cutOff) %>%
    dplyr::select(proteinID)
}
undiff.75pc <- filterCuts(df.undiff, 4, type = quart)
diff.75pc <- filterCuts(df.diff, 4, type = quart)
undiff.90pc <- filterCuts(df.undiff, 10, type = decile)
diff.90pc <- filterCuts(df.diff, 10, type = decile)
heat_colors <- brewer.pal(9, "YlOrRd")
plotTE.heatmap.full <-
  function(inputGenes = gene.list,
           inputTissue = "E14.5-Brain",
           GeneNames = FALSE) {
    gs <-
      GeneSet(geneIds = inputGenes,
              organism = "Mus Musculus",
              geneIdType = SymbolIdentifier())
    output <- teEnrichment(inputGenes = gs, rnaSeqDataset = 3)
    en.output <-
      setNames(data.frame(assay(output[[1]]),
                          row.names = rowData(output[[1]])[, 1]),
               colData(output[[1]])[, 1])
    en.output$Tissue <- rownames(en.output)
    seExp <- output[[2]][[inputTissue]]
    exp <-
      setNames(data.frame(assay(seExp), row.names = rowData(seExp)[, 1]),
               colData(seExp)[, 1])
    g <- pheatmap(
      exp,
      color = heat_colors,
      cluster_rows = F,
      cluster_cols  = T,
      show_rownames = GeneNames,
      border_color = NA,
      fontsize = 10,
      scale = "row",
      fontsize_row = 8
    )
    g
  }

```

# Enrichment plots (`TissueEnrich`) {.tabset}

## pTGC

```{r te1, fig.cap="Figure 3A: Differentially expressed proteins (up-regulated in pTGC) enrichment in Mouse ENCODE data", fig.width=8, fig.height=6}
plotTE(unique(diff), myColor = "#c6007b")
```

## TSC

```{r te2, fig.cap="Figure 3B: Differentially expressed proteins (up-regulated in TSC) enrichment in Mouse ENCODE data", fig.width=8, fig.height=6}
plotTE(unique(undiff), myColor = "#a0b600")
```

## pTGC upper quartile

```{r te3, fig.cap="Figure 3C: Upper quartile proteins (pTGC) enrichment in Mouse ENCODE data", fig.width=8, fig.height=6}
plotTE(as.character(diff.75pc$proteinID), "#c6007b")
```

## TSC upper quartile

```{r te5, fig.cap="Figure 3E: Upper quartile proteins (TSC) enrichment in Mouse ENCODE data", fig.width=8, fig.height=6}
plotTE(as.character(undiff.75pc$proteinID), "#a0b600")
```


# Enrichment Heatmaps for E14.5-Placenta {.tabset}

## pTGC

```{r hmte1, fig.cap="Figure 4A: Heatmaps for E14.5-Brain tissue from the enrichment test using up-regulated proteins of pTGC", fig.width=8, fig.height=4}
plotTE.heatmap.full(
  inputGenes = unique(diff),
  inputTissue = "E14.5-Placenta",
  GeneNames = TRUE
)
```

## TSC

```{r hmte2, fig.cap="Figure 4B: Heatmaps for E14.5-Brain tissue from the enrichment test using up-regulated proteins of TSC", fig.width=8, fig.height=2}
plotTE.heatmap.full(
  inputGenes = unique(undiff),
  inputTissue = "E14.5-Placenta",
  GeneNames = TRUE
)
```

# Session Information

```{r sessioninfo}
sessionInfo()
```