DuguetAnalysis_2UpepMax.Rmd

---
title: "Proteomic comparison of regulatory and conventional mouse T-cells (Fanny Duguet)"
author: "Marie Locard-Paulet"
date: '`r Sys.Date()`'
output: html_document
---

```{r, include=FALSE}
require(ggplot2)
require(gplots)
library(corrplot)
require(reshape2)
require(gridExtra)
require(knitr)

colfunc <- colorRampPalette(c("darkblue", "blue", "deepskyblue", "white", "yellow", "orange", "red", "darkred"))
```

```{r multiplot, echo = F}
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
```

```{r pairwise plots, echo=F}
panel.cor <- function(x, y, digits = 2, cex.cor, ...)
{ # modified from https://gist.github.com/arsalvacion/1ba2373bbe89b2d3c023#file-pairscor-r
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  # correlation coefficient
  r <- cor(x, y, use = "pairwise.complete.obs")
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste("r= ", txt, sep = "")
  text(0.5, 0.6, txt)

  # p-value calculation
  p <- cor.test(x, y)$p.value
  txt2 <- format(c(p, 0.123456789), digits = digits)[1]
  txt2 <- paste("p= ", txt2, sep = "")
  if(p<0.01) txt2 <- paste("p= ", "<0.01", sep = "")
  text(0.5, 0.4, txt2)
}
``` 

# The data

This document presents the analysis of Fanny Duguet's data set. It is constituted of 7 biological replicates, with mouse Foxp3+ Tregs (that include both CD25+ and CD25- Tregs) or of Foxp3- naïve T cells, refered to as regulatory T-cells (Treg) and conventional T-cells (Tconv), respectively. All the files have been analysed in Maxquant using matching between runs (groupped per experiment) and normalised across all samples for label free quantification. The peptides with a peptide score under 20 have been removed from the analysis. I work with the output table called `proteinGroups-Score15.txt` located in the folder `/Raw/Duguet`.

Here are the different steps of the analysis:

1. I remove all the contaminants by getting rid of the `Majority.protein.IDs` containing `REV\_` or `CON\_`. I also remove the proteins identified with less than 2 unique peptides.
2. The quantification is performed with the LFQ normalised intensity values (column names starting  with `LFQ.intensity`).
3. 0 are replaced by `NAs`. 
4. Proteins with `NAs` in all conditions are removed.
5. I log2-transform the LFQ intensities.
7. All the NA values are replaced by the noise (1% quantile of all the LFQ values of the file).
8. I keep the proteins where there are at least 3 values from independent experiments in the Tconv or the Treg.
9. The means per biological replicate are calculated for each cell population.
10. A paired 2-sided t-test is performed to identify proteins differencially regulated in the 2 cell populations.
11. Proteins with a p-value under or equal to 0.05 and with a fold change (log2(Treg/Tconv)) above 1 or under -1 (i.e. 2-fold regulation) are selected, their p-values are adjusted (BH) and the proteins with an adjusted-pvalue under or equal to 0.05 are considered significantly regulated.
 
The biological replicates are independent experiments. These are labelled:

* M5 (2 injection replicates) - Replicate 1
* M6 (3 injection replicates) - Replicate 2
* M11_R1 (3 injection replicates) - Replicate 3
* M11_R2 (3 injection replicates) - Replicate 4
* M12_1 (1 injection replicate) - Replicate 5
* M12_2 (1 injection replicate) - Replicate 6
* M12_3 (1 injection replicate) - Replicate 7

A file `FilesList.txt` in the folder `MappingTables/` allows the mapping of all raw files, gradient times and replicate names for the analysis.

```{r input, include = F}
mapping <- read.delim2(file = "MappingTables/FilesList.txt", header = T, sep = "\t")
```

`r kable(mapping)`

```{r input2, include = F}
tab <- read.delim("./Raw/Duguet/proteinGroups-Score15.txt", stringsAsFactors = F, sep = "\t")
input <- length(unique(tab$Majority.protein.IDs))
tab <- tab[!grepl("CON_", tab$Majority.protein.IDs) & !grepl("REV_", tab$Majority.protein.IDs),]
input2 <- length(unique(tab$Majority.protein.IDs))
vec <- sapply(tab$Peptide.counts..unique., function(x){strsplit(x, ";", fixed = T)[[1]][1]})
tab <- tab[as.numeric(vec)>=2,]
input3 <- length(unique(tab$Majority.protein.IDs))
mat <- tab[,grepl("LFQ.", names(tab))]
mat[mat==0] <- NaN
row.names(mat) <- tab$Majority.protein.IDs
mat <- log2(mat)

ms <- tab[,grepl("Identification.type", names(tab))]
row.names(ms) <- tab$Majority.protein.IDs

log <- sapply(1:nrow(mat), function(x) sum(is.na(mat[x,])) != length(mat[x,]))
mat <- mat[log,]
ms <- ms[log,]

dimnames(mat)[[2]] <- gsub("M5", "1", dimnames(mat)[[2]])
dimnames(mat)[[2]] <- gsub("M6", "2", dimnames(mat)[[2]])
dimnames(mat)[[2]] <- gsub("M11_R1", "3", dimnames(mat)[[2]])
dimnames(mat)[[2]] <- gsub("M11_R2", "4", dimnames(mat)[[2]])
dimnames(mat)[[2]] <- gsub("M12_1", "5", dimnames(mat)[[2]])
dimnames(mat)[[2]] <- gsub("M12_2", "6", dimnames(mat)[[2]])
dimnames(mat)[[2]] <- gsub("M12_3", "7", dimnames(mat)[[2]])
```

From the MaxQuant export table, there are `r input` unique "Majority.protein.IDs". After removing the contaminants, `r input2` remain, and `r input3` remain when we keep the protein with a minimum of 2 unique peptide. After removing the rows with only NAs, there are `r nrow(mat)` unique "Majority.protein.IDs" remaining.

Pairwise plots of the raw values per biological experiment with technical replicates:

```{r paiwise, echo=FALSE}
#Exp <- c("M5", "M6", "M11_R1", "M11_R2")
Exp <- c(".1", ".2", ".3", ".4", ".5", ".6", ".7")
pdf("Figures/2PepOrMore/QCTechnicalByReplicates.pdf", useDingbats=FALSE, 11.69, 8.27) # width and height in inches.
for (el in Exp) {
  sub <- mat[,grepl(el, dimnames(mat)[[2]])]
  pairs(sub, pch = ".", main = paste0("Biological replicate ", el), upper.panel = panel.cor)
}
dev.off()

for (el in Exp) {
  sub <- mat[,grepl(el, dimnames(mat)[[2]])]
  pairs(sub, pch = ".", main = paste0("Biological replicate ", el), upper.panel = panel.cor)
}

M <- cor(mat, use = "complete.obs")
col1 <- colorRampPalette(c("darkred", "red", "indianred1", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "cornflowerblue", "blue", "darkblue", "black")) 
corrplot.mixed(M, lower = "square", upper = "number", col=col1(600)) 

pdf("Figures/2PepOrMore/QCTechnicalReplicates.pdf", useDingbats=FALSE, 11.69, 8.27) # width and height in inches.
corrplot.mixed(M, lower = "square", upper = "number", col=col1(600)) 
dev.off()
dev.off()
```

```{r ReplacementNoise, echo=FALSE}
par(mar = c(15, 3, 4, 1))
boxplot(mat, las = 2, main = "log2-transformed LFQ intensities\nbefore replacing missing values")

par(mar = c(5, 3, 1, 1))
noise <- sapply(1:ncol(mat), function(x) quantile(mat[,x], probs=0.01, na.rm = T))
```

1% quantile per file used to replace missing values:

```{r, echo=FALSE}
kable(data.frame("file" = dimnames(mat)[[2]], "Replacement value" = noise))
#for (i in 1:length(noise)) {
#  hist(mat[,i], breaks = 30, main = paste0("Noise value (1% quantile): ", dimnames(mat)[[2]][i]))
#  abline(v = noise[i], col = "red")
#}
mat <- mat[,order(names(mat))]

#Exp2 <- c("M5", "M6", "M11_R1", "M11_R2", "M12_1", "M12_2", "M12_3")
Exp2 <- Exp

vec <- rep(FALSE, nrow(mat))
for (i in 1:nrow(mat)) {
  lconvVal <- mat[i, grepl("Tconv", dimnames(ms)[[2]])]
  lregVal <- mat[i, grepl("Treg", dimnames(ms)[[2]])]
  v <- sapply(Exp2, function(x) lconvVal[grepl(x, names(lconvVal))])
  conv <- sapply(v, function(x) sum(!is.na(x)))
  v <- sapply(Exp2, function(x) lregVal[grepl(x, names(lregVal))])
  reg <- sapply(v, function(x) sum(!is.na(x)))
  if (length(conv[conv!=0]) >= 3 | length(reg[reg]) >= 3) {vec[i] <- TRUE}
}

mat <- mat[vec,]
filter1 <- nrow(mat)

mat2 <- mat
for (i in 1:ncol(mat)) {
  v <- mat[,i]
  v[is.na(v)] <- noise[i]
  mat2[,i] <- v
}

mat <- mat2
```

After having kept only the proteins having a minimum of 3 values from independent biological condition (Tconv or Treg), `r filter1` proteins remain.

```{r PCA, echo=FALSE}
par(mar = c(15, 3, 4, 1))
boxplot(mat, las=2, main = "log2(LFQ) intensities after\nreplacement of missing values")
```

# Mean per biological replicate

The table with the mean values per condition and experiment is saved as `Duguet_Means.txt` in the folder `/OutputTables`. I make the pairwise plots out of curiosity. 

```{r Means, echo=FALSE}
Means <- matrix(nrow = nrow(mat), ncol = 14)
conditions <- c("Tconv.M11_R1", "Tconv.M11_R2", "Tconv.M12_1", "Tconv.M12_2", "Tconv.M12_3", "Tconv.M5", "Tconv.M6", "Treg.M11_R1", "Treg.M11_R2", "Treg.M12_1", "Treg.M12_2", "Treg.M12_3", "Treg.M5", "Treg.M6")
conditions <- gsub("M5", "1", conditions)
conditions <- gsub("M6", "2", conditions)
conditions <- gsub("M11_R1", "3", conditions)
conditions <- gsub("M11_R2", "4", conditions)
conditions <- gsub("M12_1", "5", conditions)
conditions <- gsub("M12_2", "6", conditions)
conditions <- gsub("M12_3", "7", conditions)

dimnames(Means)[[2]] <- conditions
dimnames(Means)[[1]] <- row.names(mat)
for (i in 1:ncol(Means)) {
  cond <- conditions[i]
  sub <- mat[,grepl(cond, dimnames(mat)[[2]], fixed = T)]
  if (is.data.frame(sub)) {
    m <- rowMeans(sub, na.rm=T)
  } else {
    m <- sub
  }
  Means[,i] <- m
}

par(mar = c(9, 3, 3, 1))
boxplot(Means, las = 2, main = "Mean LFQ values per experiment and condition")

#write.table(cbind(Means, "Protein" = row.names(Means)), "OutputTables/2PepOrMore/Duguet_Means.txt", sep = "\t", row.names = F)

Exp <- c("Tconv", "Treg")
for (el in Exp) {
  sub <- Means[,grepl(el, dimnames(Means)[[2]])]
  pairs(sub, pch = ".", main = el, upper.panel = panel.cor)
}

pdf("Figures/2PepOrMore/QCMeanPerCondition.pdf", useDingbats=FALSE, 8.27, 8.27) # width and height in inches.
Exp <- c("Tconv", "Treg")
for (el in Exp) {
  sub <- Means[,grepl(el, dimnames(Means)[[2]])]
  pairs(sub, pch = ".", main = el, upper.panel = panel.cor)
}
dev.off()
```

CV values for QC plots:

```{r, echo = F}
CV <- matrix(nrow = nrow(mat), ncol = 14)
conditions <- c("Tconv.M11_R1", "Tconv.M11_R2", "Tconv.M12_1", "Tconv.M12_2", "Tconv.M12_3", "Tconv.M5", "Tconv.M6", "Treg.M11_R1", "Treg.M11_R2", "Treg.M12_1", "Treg.M12_2", "Treg.M12_3", "Treg.M5", "Treg.M6")
conditions <- gsub("M5", "1", conditions)
conditions <- gsub("M6", "2", conditions)
conditions <- gsub("M11_R1", "3", conditions)
conditions <- gsub("M11_R2", "4", conditions)
conditions <- gsub("M12_1", "5", conditions)
conditions <- gsub("M12_2", "6", conditions)
conditions <- gsub("M12_3", "7", conditions)

dimnames(CV)[[2]] <- conditions
dimnames(CV)[[1]] <- row.names(mat)
for (i in 1:ncol(CV)) {
  cond <- conditions[i]
  sub <- mat[,grepl(cond, dimnames(mat)[[2]], fixed = T)]
  if (is.data.frame(sub)) {
    m <- sapply(1:nrow(sub), function(x) sd(sub[x,], na.rm=T))
  } else {
    m <- rep(NA, nrow(CV))
  }
  CV[,i] <- m
}
CV <- CV/Means
hist(CV)
medCV <- median(CV, na.rm = T)
q099 <- quantile(CV, 0.99, na.rm = T)
q090 <-quantile(CV, 0.9, na.rm = T)


ggplot(melt(CV), aes(x = value)) + geom_histogram() + theme_bw() + ggtitle("histogram of CVs obtained from technical replicates")
pdf("Figures/2PepOrMore/QCCVs.pdf", useDingbats=FALSE, 8.27, 4) # width and height in inches.
ggplot(melt(CV), aes(x = value)) + geom_histogram() + theme_bw() + ggtitle("histogram of CVs obtained from technical replicates")
dev.off()
```

The median of technical CVs is `r medCV*100`, 99% of the CVs are under `r q099*100` and 90% are under `r q090*100`.

Correlation between biological replicates and conditions:

```{r, echo = F}
#col1 <- colorRampPalette(c("red", "white", "blue")) 
M <- cor(Means, use = "complete.obs")
corrplot.mixed(M, lower = "circle", upper = "number") # , col=col1(100)
pairs(Means, pch = ".", upper.panel = panel.cor)

pdf("Figures/2PepOrMore/QCMeanReplicates.pdf", useDingbats=FALSE, 8.27, 8.27) # width and height in inches.
pairs(Means, pch = ".", upper.panel = panel.cor)
dev.off()
```

CV between biological replicates and conditions:

```{r, echo = F}
CV <- matrix(nrow = nrow(mat), ncol = 2)
conditions <- c("Tconv", "Treg")

dimnames(CV)[[2]] <- conditions
dimnames(CV)[[1]] <- row.names(mat)
for (i in 1:ncol(CV)) {
  cond <- conditions[i]
  sub <- mat[,grepl(cond, dimnames(mat)[[2]], fixed = T)]
  if (is.data.frame(sub)) {
    m <- sapply(1:nrow(sub), function(x) sd(sub[x,], na.rm=T))
    m2 <- sapply(1:nrow(sub), function(x) mean(as.numeric(sub[x,]), na.rm=T))
    m <- m/m2
  } else {
    m <- rep(NA, nrow(CV))
  }
  CV[,i] <- m
}
hist(CV)

medCV <- median(CV, na.rm = T)
q099 <- quantile(CV, 0.99, na.rm = T)
q090 <-quantile(CV, 0.9, na.rm = T)



ggplot(melt(CV), aes(x = value)) + geom_histogram() + theme_bw() + ggtitle("histogram of CVs obtained between means of biological replicates")
pdf("Figures/2PepOrMore/QCCVMeans.pdf", useDingbats=FALSE, 8.27, 4) # width and height in inches.
ggplot(melt(CV), aes(x = value)) + geom_histogram() + theme_bw() + ggtitle("histogram of CVs obtained between means of biological replicates")
dev.off()
```

The median of biological CVs is `r medCV*100`, 99% of the CVs are under `r q099*100` and 90% are under `r q090*100`.

# Statistical analysis

In the table with mean values, there are `r nrow(Means)` protein quantified.


```{r stat1, echo = F}
MeansFiltered <- Means
```

We then calculate the pvalue of two-sided paired t-test. The p-values from proteins that have a minimum of 2-fold change between Treg and Tconv (log2(mean(Treg))-log2(mean(Tconv)) above or equal to 1 or log2(mean(Treg))-log2(mean(Tconv)) under or equat to -1) are adjusted (BH). The fold change (FC) coresponds to log2(mean(LFQ in Treg) / Mean(LFQ in Tconv)). Everything is saved with the raw table in the table `OutputTables/2PepOrMore/TotalAnalysisTconvTreg.txt`.

```{r stat2, echo = F}
stat <- matrix(nrow = nrow(MeansFiltered), ncol = 3)
dimnames(stat)[[2]] <- c("FC", "pval", "adjpval")
dimnames(stat)[[1]] <- dimnames(MeansFiltered)[[1]]
for (i in 1:nrow(MeansFiltered)) {
  reg <- MeansFiltered[i, grepl("Treg", dimnames(MeansFiltered)[[2]])]
  conv <- MeansFiltered[i, grepl("Tconv", dimnames(MeansFiltered)[[2]])]
  if(length(reg[!is.na(reg)])>=2 & length(conv[!is.na(conv)])>=2) {
    vec <- is.na(reg)+is.na(conv)
    if(length(vec[vec==0])>=2) {
      stat[i,1] <- mean(reg, na.rm = T)-mean(conv, na.rm = T)
      stat[i,2] <- t.test(x = reg, y = conv, alternative = "two.sided", paired = T)$p.value
    }
  }
}
ap <- stat[,2][abs(stat[,1]) >= 1]
ap2 <- p.adjust(ap, method = "BH")
stat[,3] <- ap2[match(rownames(stat), names(ap2))]

df <- data.frame(stat, "Majority.protein.IDs" = row.names(stat))
df2 <- data.frame(Means, "Majority.protein.IDs" = row.names(Means))
export <- merge(df, df2, by = "Majority.protein.IDs")
export <- merge(export, tab, by = "Majority.protein.IDs")
#write.table(export, "OutputTables/2PepOrMore/TotalAnalysisTconvTreg.txt", row.names=F, sep = "\t")
```

```{r stat3, echo=FALSE}
sign <- stat[((stat[,1]>=1) | (stat[,1]<=-1)) & (stat[,3] <= 0.05),]
sign <- sign[complete.cases(sign),]

sign <- data.frame(sign, "Majority.protein.IDs" = dimnames(sign)[[1]])
signM <- merge(sign, tab, by = "Majority.protein.IDs")
#write.table(signM, "OutputTables/2PepOrMore/SignificativeHits.txt", row.names = F, sep = "\t")
```

## Volcano plot:

For this volcano plot, I use marker names from the file `MappingTables/MappingMarkers.tab` for colours and shapes. The idea is to localise proteins known to be differentially regulated in the Terg Vs Tconv in the bulk of our proteins.

The output figure is saved as a .pdf: `Figures/2PepOrMore/Volcano.pdf`.

```{r, echo = F}
IDMapper <- matrix(ncol = 2, nrow = nrow(tab))
dimnames(IDMapper)[[2]] <- c("Majority.Protein.IDs", "Gene.names")
IDMapper[,1] <- tab$Majority.protein.IDs
IDMapper[,2] <- tab$Gene.names
#IF WE WANT TO USE THE PROTEIN NAMES INSTEAD OF GENE NAMES:
#pn <- sapply(1:nrow(tab), function(x) strsplit(tab$Fasta.headers[x], "_MOUSE", fixed = T)[[1]][1])
#pn <- sapply(1:nrow(tab), function(x) strsplit(pn[x], "|", fixed = T)[[1]][3])
#IDMapper[,2] <- pn

map <- read.delim("./MappingTables/MappingMarkers.tab", header = T, row.names = NULL)
candidate <- as.character(unique(map$Entry[map$value=="candidate"]))
down <- as.character(unique(map$Entry[map$value=="down"]))
up <- as.character(unique(map$Entry[map$value=="up"]))
map <- cbind(map, "pn" = IDMapper[match(map$Entry, IDMapper[,1]),2])

df <- as.data.frame(stat)

par(mfrow=c(3,2))
title <- "Treg Vs Tconv"
col <- ifelse(df$adjpval<=0.05 & df$FC >=1, "up", "no regulation")
col[df$adjpval<=0.05 & df$FC <= -1] <- "down"
numbers <- table(col)
shape <- rep("20", nrow(df))
col[col=="down"] <- paste0("down-regulated", " (", numbers[1], ")")
col[col=="no regulation"] <- paste0("no regulation", " (", numbers[2], ")")
col[col=="up"] <- paste0("up-regulated", " (", numbers[3], ")")
col[row.names(df) %in% up] <- "Markers up"
col[row.names(df) %in% down] <- "Markers down"
col[row.names(df) %in% candidate] <- "THEMIS"
shape[row.names(df) %in% up] <- "21"
shape[row.names(df) %in% down] <- "21"
shape[row.names(df) %in% candidate] <- "21"
text <- rep(NA, nrow(df))
text[row.names(df) %in% candidate] <- "THEMIS"
text[row.names(df) %in% down] <- as.character(map$pn[match(row.names(df)[row.names(df) %in% down], map$Entry)])
text[row.names(df) %in% up] <- as.character(map$pn[match(row.names(df)[row.names(df) %in% up], map$Entry)])


df <- data.frame(df, "colour" = col, "shape"=shape, "label"=text)

g <- ggplot(df, aes(x=FC, y=-log(pval,10), col = colour, shape = factor(shape), fill = colour, label = label)) + theme_minimal() + geom_vline(xintercept = 0) + geom_point(alpha = 0.65, size = 4) + geom_hline(yintercept = 0) + geom_point(alpha = 0.65, size = 4) + scale_colour_manual(values = c("deepskyblue", "darkblue", "darkred", "black", "black", "darkorange")) + scale_fill_manual(values = c("deepskyblue", "blue", "darkred", "black", "yellow", "darkorange")) + xlab("log2(Fold change)") + xlim(c(-6.4, 6.4)) + ylim(c(0,8.6)) + ggtitle(title) + scale_shape_manual(values = c(20, 21)) + geom_text(hjust=0.5, vjust=-0.8)
print(g)

pdf("Figures/2PepOrMore/Volcano.pdf", useDingbats=FALSE, 11.69, 8.27) # width and height in inches.
g
dev.off()
```

## Bar plot figures with most regulated proteins:

These two plots present the 40 most regulated proteins in the Treg Vs Tconv (up in orange and down in blue). The proteins known to present a different quantity in the two cell populations are highlighted.

The figures are saved as .pdf in the folder `/Figures/2PepOrMore`: `BarPlotUp.pdf` and `BarPlotDown.pdf`.


```{r, echo=F}
# Make bar plot:
df <- signM[order(signM$FC),]
df <- data.frame(df, "GeneName" = IDMapper[match(df$Majority.protein.IDs, IDMapper[,1]),2])
df$GeneName <- factor(as.character(df$GeneName), levels = as.character(df$GeneName))

col <- rep("NA", nrow(df))
col[df$Majority.protein.IDs %in% up] <- "Markers up"
col[df$Majority.protein.IDs %in% down] <- "Markers down"
col[df$Majority.protein.IDs %in% candidate] <- "THEMIS"

df <- data.frame(df, "colour"=col)

p1 <- df[(df$FC >= 1) & (df$adjpval <= 0.05),]
p2 <- df[(df$FC <= -1) & (df$adjpval <= 0.05),]
p1 <- p1[((nrow(p1)-39):nrow(p1)),]
p2 <- p2[1:40,]
print(ggplot(p1, aes(x = GeneName, y = FC, fill = colour)) + theme_bw() + geom_bar(stat = "identity", width = 0.6) + theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 5)) + scale_fill_manual(values = c("darkred", "darkorange"))  + scale_y_continuous(expand=c(0, 0.1)) + ggtitle("Upregulated proteins"))
print(ggplot(p2, aes(x = GeneName, y = FC, fill = colour)) + theme_bw() + geom_bar(stat = "identity", width = 0.6) + theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 5)) + scale_fill_manual(values = c("darkblue", "dodgerblue2", "gold"))  + scale_y_continuous(expand=c(0, 0.1)) + ggtitle("Downregulated proteins"))

pdf("Figures\\2PepOrMore\\BarPlotUp.pdf", useDingbats=FALSE, 11.69, 8.27) # width and height in inches.
ggplot(p1, aes(x = GeneName, y = FC, fill = colour)) + theme_bw() + geom_bar(stat = "identity", width = 0.6) + theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) + scale_fill_manual(values = c("darkred", "darkorange"))  + scale_y_continuous(expand=c(0, 0.1)) + ggtitle("Upregulated proteins")
dev.off()
pdf("Figures\\2PepOrMore\\BarPlotDown.pdf", useDingbats=FALSE, 11.69, 8.27) # width and height in inches.
ggplot(p2, aes(x = GeneName, y = FC, fill = colour)) + theme_bw() + geom_bar(stat = "identity", width = 0.6) + theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) + scale_fill_manual(values = c("darkblue", "dodgerblue2", "gold"))  + scale_y_continuous(expand=c(0, 0.1)) + ggtitle("Downregulated proteins")
dev.off()

```

I generate boxplots with the 80 proteins in the barplots (highest fold changes).

```{r, echo = F, message=F}
mat1 <- export[,c(24,5:18)]
mat1 <- mat1[mat1$Gene.names %in% p1$Gene.names,]

j <- 3
n <- 1
for (k in round(seq(from = 1, to = (nrow(mat1)-(nrow(mat1)/j)+1), by = nrow(mat1)/j), 0))  {
  plots <- list()
  mat2 <- mat1[mat1$Gene.names %in% mat1$Gene.names[k:(k+13)],]
  for (i in 1:nrow(mat2)) {
    title <- as.character(mat2$Gene.names[i])
    mat <- mat2[i,2:(ncol(mat2))]
    gtab <- melt(mat)
    vec <- ifelse(grepl("Tconv", gtab$variable), "Tconv", "Treg")
    gtab <- cbind(gtab, "CellType" = vec)
    g <- ggplot(gtab, aes(x = CellType, y = log2(value), col = CellType)) + geom_boxplot(width = 0.8) + geom_point(size = 1) + theme_bw() + xlab("") + ylab("LFQ value per replicate") + ggtitle(title) + scale_colour_manual(values = c("black", "red3")) + theme(legend.position="none", axis.text=element_text(size=8), axis.title=element_text(size=11,face="bold"))
    plots[[length(plots)+1]] <- g
  }
  na <- paste0("Figures\\2PepOrMore\\MultipleBoxPlotsUp", n, ".pdf")
  pdf(na, useDingbats=FALSE, 8.27, 11.69) # width and height in inches.
  multiplot(plotlist = plots, cols = 4)
  dev.off()
  n <- n+1
}

mat1 <- export[,c(24,5:18)]
mat1 <- mat1[mat1$Gene.names %in% p2$Gene.names,]

j <- 3
for (k in round(seq(from = 1, to = (nrow(mat1)-(nrow(mat1)/j)+1), by = nrow(mat1)/j), 0))  {
  plots <- list()
  mat2 <- mat1[mat1$Gene.names %in% mat1$Gene.names[k:(k+13)],]
  for (i in 1:nrow(mat2)) {
    title <- as.character(mat2$Gene.names[i])
    mat <- mat2[i,2:(ncol(mat2))]
    gtab <- melt(mat)
    vec <- ifelse(grepl("Tconv", gtab$variable), "Tconv", "Treg")
    gtab <- cbind(gtab, "CellType" = vec)
    g <- ggplot(gtab, aes(x = CellType, y = log2(value), col = CellType)) + geom_boxplot(width = 0.8) + geom_point(size = 1) + theme_bw() + xlab("") + ylab("LFQ value per replicate") + ggtitle(title) + scale_colour_manual(values = c("black", "royalblue4")) + theme(legend.position="none", axis.text=element_text(size=8), axis.title=element_text(size=11,face="bold"))
    plots[[length(plots)+1]] <- g
  }
  na <- paste0("Figures\\2PepOrMore\\MultipleBoxPlotsDown", n, ".pdf")
  pdf(na, useDingbats=FALSE, 8.27, 11.69) # width and height in inches.
  multiplot(plotlist = plots, cols = 4)
  dev.off()
  n <- n+1
}

```

I save the individual boxplots in `Figures/2PepOrMore/MultipleBoxplotsUp.pdf` and `Figures/2PepOrMore/MultipleBoxplotsDown.pdf`.

I create a table for the paper:

```{r, echo = F}
tablePaper <- df[,c(9,10,2,3,4)]
write.table(tablePaper, file = "OutputTables/2PepOrMore/Table1Paper.txt", sep = "\t", row.names = F)
```

## Heatmap of raw log2(LFQ) values:

This heatmap presents the genes with an absolute log2-transformed fold change above or equal to 1 and a p-value under or equal to 0.05.
This heatmap is ordered in function of the sum of the LFQ intensities in the Treg, decreasing. It is named `Figures/2PepOrMore/TileMap.pdf`.

```{r, echo=F}
df <- export[,c(1:18, 24)]
#df <- data.frame(df, "GeneName" = IDMapper[match(df$Majority.protein.IDs, IDMapper[,1]),2])
sub <- df[,grepl("Treg", names(df))]
s <- rowSums(sub)
df <- data.frame(df, "IntensityTreg"=s)
df <- df[((df$FC>=1) | (df$FC <= -1)) & (df$adjpval <= 0.05),]
df <- df[!is.na(df$Majority.protein.IDs),]
df$Gene.names <- factor(as.character(df$Gene.names), levels = df$Gene.names[order(df$IntensityTreg)])
df2 <- melt(df[,c(1:19)])
df2 <- df2[df2$variable != "pval",]
df2 <- df2[df2$variable != "adjpval",]
df2[is.na(df2)] <- NA

dftop <- df2[df2$Majority.protein.IDs %in% export$Majority.protein.IDs[export$FC >= 1],]
dfbottom <- df2[df2$Majority.protein.IDs %in% export$Majority.protein.IDs[export$FC <= -1],]
# colours = c("darkblue", "blue", "deepskyblue", "white", "orange", "red", "darkred")
colours = c("blue", "gray87", "red")
#colours = c("deepskyblue", "black", "red")
# colours = c("green4", "yellow", "red")
#colours = c("darkgreen", "yellow", "red")

top <- ggplot(dftop[dftop$variable != "FC",], aes(y = Gene.names, x = variable)) + theme_classic() + geom_tile(aes(fill = value)) + scale_fill_gradientn(colours = colours) + theme(axis.title.x = element_blank()) + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank(), axis.text.y = element_text(size = 5)) + geom_vline(xintercept = 7.5, size = 1) + ggtitle("log2(LFQ intensities)")

bottom <- ggplot(dfbottom[dfbottom$variable != "FC",], aes(y = Gene.names, x = variable)) + theme_classic() + geom_tile(aes(fill = value))  + scale_fill_gradientn(colours = colours) + theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_text(size = 5)) + geom_vline(xintercept = 7.5, size = 1)


plots <- list(top, bottom) # to align axis
grobs <- list()
widths <- list()
for (i in 1:length(plots)){
    grobs[[i]] <- ggplotGrob(plots[[i]])
    widths[[i]] <- grobs[[i]]$widths[2:5]
}
maxwidth <- do.call(grid::unit.pmax, widths)
for (i in 1:length(grobs)){
     grobs[[i]]$widths[2:5] <- as.list(maxwidth)
}

grid.arrange(top, bottom, ncol=1, grobs = grobs, heights=c(nrow(dftop), nrow(dfbottom)*1.7))

pdf("Figures\\2PepOrMore\\TileMap.pdf", useDingbats=FALSE, 8.27, 11.69) # width and height in inches (A4 portrait).
grid.arrange(top, bottom, ncol=1, grobs = grobs, heights=c(nrow(dftop), nrow(dfbottom)*1.7))
dev.off()
```

I make the same figure with only the 40 most regulated proteins up and down (corresponds to the proteins in the barcharts). This is called `Figures/2PepOrMore/TileMapRegulated.pdf`.

```{r, echo = F}
genes <- unique(c(as.character(p1$GeneName), as.character(p2$GeneName)))

top <- ggplot(dftop[dftop$variable != "FC" & dftop$Gene.names %in% genes,], aes(y = Gene.names, x = variable)) + theme_classic() + geom_tile(aes(fill = value)) + scale_fill_gradientn(colours = colours) + theme(axis.title.x = element_blank()) + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank(), axis.text.y = element_text(size = 8)) + geom_vline(xintercept = 7.5, size = 1) + ggtitle("log2(LFQ intensities)")

bottom <- ggplot(dfbottom[dfbottom$variable != "FC" & dfbottom$Gene.names %in% genes,], aes(y = Gene.names, x = variable)) + theme_classic() + geom_tile(aes(fill = value))  + scale_fill_gradientn(colours = colours) + theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_text(size = 8)) + geom_vline(xintercept = 7.5, size = 1)

grid.arrange(top, bottom, ncol=1, heights=c(40, 46))

pdf("Figures\\2PepOrMore\\TileMapRegulated.pdf", useDingbats=FALSE, 8.27, 11.69) # width and height in inches (A4 portrait).
grid.arrange(top, bottom, ncol=1, heights=c(40, 46))
dev.off()
```
There are `r nrow(sign)` proteins with a FC <= -1 or >= 1 (2-fold change between the 2 cell populations) and an adjusted p-value under or equal to 0.05. These are saved in the table `outputTable\SignificativeHits.txt`.

All the figures generated in this document are saved in `Figures\`.

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