---
title: "8_3_manual_clustering_immune_genes.Rmd"
author: "Aspen Coyle"
date: "8/18/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Introduction

This script is based off a script Laura Spencer wrote to take single-crab libraries and separate them into modules over time.

Inputs are libraries pseudoaligned (with kallisto) with cbai_transcriptomev2.0, which is unfiltered by taxa, and then filtered to only include sequences with GO terms linked to immune response.

Unlike in previous manual clustering scripts, we have no need to assemble counts of TPMs - we've already done that in an earlier script for immune-linked transcripts!

Original script is available in [Lab Github Discussion 1206](https://github.com/RobertsLab/resources/discussions/1206)

## Load libraries, install all required new libraries

```{r libraries, message=FALSE, warning=FALSE}
# Add all required libraries here
list.of.packages <- c("pheatmap", "tidyverse")
# Get names of all required packages that aren't installed
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
# Install all new packages
if(length(new.packages)) install.packages(new.packages)


# Load all required libraries
lapply(list.of.packages, FUN = function(X) {
  do.call("require", list(X))
})
```

## Set our inputs

```{r}
# Give IDs for all libraries for ambient-temperature crab. All these are infected crab.
crabA_libs <- c("178", "359", "463")
crabB_libs <- c("118", "349", "481")
crabC_libs <- c("132", "334", "485")
# Also giving IDs for all libraries for decreased-temperature crab. D and F are uninfected, E is infected. 
crabD_libs <- c("073", "221", "427")
crabE_libs <- c("151", "254", "445")
crabF_libs <- c("113", "222", "425")
# Give IDs for all libraries for elevated-temperature crab. Note: again, we have no libraries on Day 17 for these crabs
crabG_libs <- c("173", "272")
crabH_libs <- c("072", "294")
crabI_libs <- c("127", "280")

# Vector of crab IDs
crabs <- c("A", "B", "C", "D", "E", "F", "G", "H", "I")

# List of all crabs
crabIDs <- list(crabA_libs, crabB_libs, crabC_libs, crabD_libs, crabE_libs, crabF_libs, crabG_libs, crabH_libs, crabI_libs)

# Set path to TPM output folder
TPM_outpath <- "../output/TPM_counts/cbai_transcriptomev2.0/"

# Set path to general heatmap output folder
heatmap_output <- "../output/manual_clustering/cbai_transcriptomev2.0/immune_genes/"

```


# Create heat maps for each crab (setting specific cut height)

Each crab will be clustered  over time. These will be given one of five assignments based on their expression patterns. Multiple clusters may be given the same assignment.

- HTL = expression goes from high to low

- LTH = expression goes from low to high

- HLH = expression starts high, then goes low, then goes high again

- LHL = inverse of HLH

- MIX = mixed - no clear pattern of expression within the module

Crabs with only two time points (crabs G, H, and I) will be clustered as well. For these, five assignments will be given. Again, multiple clusters may be given the same assignment.

- LL = expression stays low

- HH = expression stays high

- LH = expression goes from low to high

- HL = expression goes from high to low

- MIX = mixed - no clear pattern of expression within the module

```{r}

# Set height of clusters for samples with 3 time points
clusterht <- 1.8 

# Set height of clusters for samples with 2 time points
twoday_clusterht <- 2.5

# Read in the initial data file
TPMcts <- read.delim(file = paste0(TPM_outpath,
                                   "immune_genes_all_libs_TPMcts.txt"
                                   ),
                     row.names = 1,
                     check.names = FALSE)

# Remove all genes without at least 1 samples with counts of 5+
TPMcts <- TPMcts %>%
  filter_all(any_vars(. > 5))

for (i in 1:length(crabs)) {
  # Set path to output folder that will contain all heatmaps and text files
  crab_out <- paste0(heatmap_output, "Crab_", crabs[i], "/")
  
  # Select the columns that match our particular crab
  crabdat <- TPMcts[, crabIDs[[i]]]
  
  # Remove all rows without any values
  crabdat <- crabdat[rowSums(crabdat[]) > 0, ]

  # When we just have two time points, we'll make a special heat map. 
  # It uses log2 counts, and doesn't scale by row
  if (ncol(crabdat) < 3) {
    # Construct heat map
    print(pheatmap(log2(crabdat+1), cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "none",
           main = "gene counts by day",
           fontsize = 9))
    
    # Plot again to save in our file
    png(paste0(crab_out, "heatmap.png"))
    
    print(pheatmap(log2(crabdat+1), cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "none",
           main = "gene counts by day",
           fontsize = 9))
    
    dev.off()
    
    out <- pheatmap(log2(crabdat+1), cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "none",
           main = "gene counts by day",
           fontsize = 9)
    
    hc <- out$tree_row
    lbl <- cutree(hc, h = twoday_clusterht)  # split dendrogram based on height specified above
    
      # Loop through gene clusters, printing heatmap for each
    for (j in 1:length(unique(lbl))) {
      clust <- crabdat[which(lbl==j) %>% names(), ]
      # Write results to table
      write.table(clust, file = paste0(crab_out,
                                       "cluster_",
                                       j, ".txt"),
                  sep = "\t",
                  row.names = TRUE)
      
      # Create heat map for each cluster. Again, like the graph above, we'll make it log2 scaled and not scale by row
      png(paste0(crab_out, "cluster_", j, "_heatmap.png"))
      
      print(pheatmap(log2(clust+1), cluster_rows = TRUE,
                         show_rownames = FALSE,
                     na.rm = TRUE,
                     cluster_cols = FALSE,
                     scale = "none",
                     main = paste0("gene counts by day, cluster ", j),
                     fontsize = 9))
      dev.off()
  }
    
  } else {
    # Construct heat map
    print(pheatmap(crabdat, cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "row",
           main = "gene counts by day",
           fontsize = 9))

    # Plot again to save in our file
    png(paste0(crab_out, "heatmap.png"))

    print(pheatmap(crabdat, cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "row",
           main = "gene counts by day",
           fontsize = 9))
    dev.off()
    
    out <- pheatmap(crabdat, cluster_rows = TRUE,
                show_rownames = FALSE,
                na.rm = TRUE,
                cluster_cols = FALSE,
                scale = "row",
                main = "gene counts by day",
                fontsize = 9)
    
    hc <- out$tree_row
    lbl <- cutree(hc, h = clusterht)  # split dendrogram based on height specified above
    
      # Loop through gene clusters, printing heatmap for each
    for (j in 1:length(unique(lbl))) {
      clust <- crabdat[which(lbl==j) %>% names(), ]
      # Write results to table
      write.table(clust, file = paste0(crab_out,
                                   "cluster_",
                                   j, ".txt"),
                  sep = "\t",
                  row.names = TRUE)
      # Create heat map for each cluster
      png(paste0(crab_out, "cluster_", j, "_heatmap.png"))
      print(pheatmap(clust, cluster_rows = TRUE,
                     show_rownames = FALSE,
                     na.rm = TRUE,
                     cluster_cols = FALSE,
                     scale = "row",
                     main = paste0("gene counts by day, cluster ", j),
                     fontsize = 9))
      dev.off()
  }
    
  }
  
}
```


### Assign labels manually to each crab

```{r}
# Create list of labels for each individual crab
clusts <- list(c("LHL", "HTL", "LTH"),  # Crab A          
               c("LHL", "LTH", "HLH", "HTL"),  # Crab B
               c("LTH", "HTL", "LHL"),  # Crab C
               c("LHL", "HTL"),  # Crab D
               c("LTH", "HTL", "HTL2", "LHL"),  # Crab E
               c("HTL", "LHL", "LTH"),  # Crab F
               c("MIX", "HL", "MIX2", "HL2", "HL3", "HL4"),  # Crab G
               c("MIX", "LL", "LH", "MIX2", "LH2"),  # Crab H
               c("HL", "MIX", "HL2", "HL3"))  # Crab I

# Ensure number of labels matches number of crab
length(clusts) == length(crabs)

# Ensure we have no duplicates in crab labels
for (i in 1:length(clusts)) {
  print(length(unique(clusts[[i]])) == length(clusts[[i]]))
}


for (i in 1:length(clusts)) {
  crab_out <- paste0(heatmap_output, "Crab_", crabs[i], "/")
  clustorder <- clusts[[i]]
  # Rename all heatmaps to include the cluster order
  file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*heatmap\\.png$")),
              paste0(crab_out, "cluster_", clustorder, "_heatmap.png"))
  # Rename all cluster filenames to include the cluster order
  file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*\\.txt")),
              paste0(crab_out, "cluster_",
                     clustorder, ".txt"))

  }
```