---
title: "Hierarchical Clustering Analysis of Proteomics Data"
author: "Mark Ziemann - Burnet Bioinformatics Group"
date: "`r Sys.Date()`"
output:
  html_document:
    theme: simplex
    toc: true
    code_folding: hide
    fig_width: 7
    fig_height: 7
    toc_float:
      collapsed: true
---

# Intro

Need to cluster the proteomics data.

to run this script, load it into Rstudio and hit "knit".

```{r,libs}

library(readxl)
library(limma)
library(gplots)
library(stringr)
library(RColorBrewer)
library(kableExtra)

```

# Read data

There were some cells with "Missing Value" and this is causing those coluns to be recognised
as characters.

I deleted those cells' contents and then loaded into R.
```{r,read_data1}

excel_sheets("1b.2_MS25-036_Dimple_Reruns._Perseus.xlsx")

xl <- read_excel("1b.2_MS25-036_Dimple_Reruns._Perseus.xlsx")

head(xl)

xl <- as.data.frame(read_excel("1b.2_MS25-036_Dimple_Reruns._Perseus.xlsx"))

str(xl)

```

Cleaning data by fixing the rownames.

```{r,clean1}

rownames(xl) <- paste(xl[,1],xl[,2])
xl[,1:2]=NULL

myrownames <- rownames(xl)

```

Remove the last 2 rows.

```{r,clean2}

xl <- xl[1:(nrow(xl) -2),]

```
# Log transformation and normalisation

```{r,norm1}

lxl <- log10(xl + 1)

nxl <- normalizeBetweenArrays(lxl, method = "quantile",
  targets = NULL, cyclic.method = "fast")

colSums(nxl,na.rm=TRUE)

```

There's a wide range of total intensity values, as expected for pollen.

# Basic visualisation

```{r,heatmap1}

colfunc <- colorRampPalette(c("blue", "white", "red"))

rownames(nxl) <- str_trunc(rownames(nxl), 20)

nxl2 <- nxl
nxl2[is.na(nxl2)] <- 0

heatmap.2(as.matrix(nxl2),trace="none",col=colfunc(25),scale="row",
  margins = c(10,6), cexRow=0.6, cexCol=0.8 )

heatmap.2(cor(as.matrix(nxl2),method="spearman"),trace="none",scale="none",
  main="Sample spearman cor", margins = c(10,10))

heatmap.2(cor(t(nxl2),method="spearman"),trace="none",scale="none",
  main="Protein spearman cor", margins = c(6,6))

```

# Hierarchical clustering

This will perform clustering with different tree cut thresholds.

The protein clusters will be saved to different files with names like "ProteinClusters_h".

```{r,hca1}

heights <- c(1.2,1.4,1.6,1.8,2,2.2)
dump <- lapply(heights,function(height) {
  cl<-as.dist(1-cor(t(nxl2), method="spearman"))
  hr <- hclust(cl , method="complete")
  mycl <- cutree(hr, h=max(hr$height/height))
  clusterCols <- colorRampPalette(c("red", "orange", "yellow","green",
    "blue","violet","black"))
  clusterCols <-  clusterCols(length(unique(mycl)))
  myClusterSideBar <- clusterCols[mycl]
  colfunc <- colorRampPalette(c("blue", "white", "red"))
  myfilename <- paste("ProteinClusters_h",height,".txt",sep="")
  write.table(mycl,file=myfilename,quote=F,sep="\t")

  heatmap.2(nxl2, main=paste("Protein Clustering h=",height),  Rowv=as.dendrogram(hr),
    dendrogram="both", scale="none", col = colfunc(25), trace="none",
    RowSideColors= myClusterSideBar, margins = c(9,7))
})

```

# Session information

For reproducibility.

```{r,session}

sessionInfo()

```
