---
title: "72_manual_clustering_cbaiv2.0"
author: "Aspen Coyle"
date: "Last compiled on `r format(Sys.time(), '%Y-%m-%d')`"
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 not filtered by taxa.

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 location of kallisto libraries
libpath <- "../output/kallisto_libraries/cbai_transcriptomev2.0/"

# 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/all_genes/"

```
### Assemble TPM counts for each crab

Loop through all crabs and all libraries within each crab, assembling the TPM counts from the abundance.tsv files to create files containing counts for each crab over time

```{r}
# In this, we're trying to create two types of files. First, a single file that contains all TPMs for all libraries. And second, a file for each crab that contains all TPMs for that particular crab. We'll do this by creating the larger file first, and then subsetting from that


# Create an empty list
all_TPMcts <- list()

for (i in 1:length(crabIDs)) {
  crabfiles <- paste0(libpath, "id", crabIDs[[i]], "/abundance.tsv")
  
  # Read in first kallisto library to start data frame
  TPMcts <- read.delim(file = crabfiles[1],
                       header = TRUE,
                       sep = "\t")
  
  # Eliminate all columns except transcript ID and TPM
  TPMcts <- TPMcts %>%
    select(target_id, tpm)
  
  # Rename columns for consistency and to ID TPM counts
  colnames(TPMcts) <- c("Transcript_ID",
                        crabIDs[[i]][1])
  
  # Loop through remaining libraries and create file of TPMs over time
  for (j in 2:length(crabIDs[[i]])) {
    kallisto_output <- read.delim(file = crabfiles[j],
                                  header = TRUE,
                                  sep = "\t")
    # Select only transcript ID and TPM columns
    kallisto_output <- kallisto_output %>%
      select(target_id, tpm)
    
    # Rename column names to give ID to count column
    colnames(kallisto_output) <- c("Transcript_ID",
                                    crabIDs[[i]][j])
    # Add TPM value to table of DEGs with full join
    TPMcts <- full_join(TPMcts, kallisto_output, 
                        by = "Transcript_ID")
  }
  
  # Write results to file
  write.table(TPMcts, file = paste0(TPM_outpath,
                                    "Crab_", 
                                    crabs[i],
                                    "_TPMcts.txt"),
              sep = "\t",
              row.names = FALSE)
  
  all_TPMcts[[i]] <- TPMcts
  
}

# Now that we've got a list with all TPM counts, put them in a matrix
all_TPMdf <- do.call("cbind", all_TPMcts)

# Remove duplicated columns (the Transcript_ID columns)
duplicated.cols <- duplicated(t(all_TPMdf))
all_TPMdf <- all_TPMdf[, !duplicated.cols]

# Write results to file
write.table(all_TPMdf, file = paste0(TPM_outpath, "all_indiv_library_TPMcts.txt"),
            sep = "\t",
            row.names = FALSE)

```




# 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 <- 10

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

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

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("LTH", "LHL", "HTL", "HTL2", "LHL2", "LTH2", "HLH"),  # Crab A
               c("LTH", "HLH", "LHL", "LHL2", "HTL", "HTL2"),  # Crab B
               c("HTL", "LTH", "HLH", "LHL", "HTL2", "LTH2"),  # Crab C
               c("HLH", "LHL", "HTL", "LTH", "HTL2", "HTL3", "LTH2"),  # Crab D
               c("LTH", "LHL", "LTH2", "HTL", "HLH", "HTL2", "LTH3"),  # Crab E
               c("LTH", "HTL", "LTH2", "HTL2", "LHL", "HLH", "HTL3"),  # Crab F
               c("MIX", "LH", "HL", "LH2", "MIX2", "MIX3"),  # Crab G
               c("MIX", "LH", "MIX2", "HL", "MIX3"),  # Crab H
               c("LH", "HL", "MIX", "LH2", "HL2", "MIX2"))  # 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)) {
  print(crabs[i])
  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"))

  }
```


#### Obsolete: Used when bar was 5 total counts of the gene per crab, not 1+ library with 30+ counts for any crab

# Create list of labels for each individual crab
clusts <- list(c("LTH", "LTH2", "LHL", "HTL", "HTL2", "HLH", "MIX"),  # Crab A
               c("LTH", "HLH", "LHL", "HTL", "LHL2", "HTL2", "HTL3"),  # Crab B
               c("LTH", "HLH", "LHL", "HTL", "HTL2", "LHL2", "LTH2"),  # Crab C
               c("HTL", "LHL", "HTL2", "HLH", "LHL2", "LTH", "LTH2"),  # Crab D
               c("HTL", "LHL", "LTH", "LTH2", "HLH", "HTL2", "LHL2"),  # Crab E
               c("LTH", "HTL", "LTH2", "HTL2", "LHL", "HLH"),  # Crab F
               c("LH", "HL", "HL2", "LH2", "MIX", "LL"),  # Crab G
               c("MIX", "HL", "MIX", "LH", "MIX", "MIX"),  # Crab H
               c("LH", "HL", "MIX", "MIX"))  # Crab I


# Create heat maps for each crab (manual cluster numbers)

This was the main method as of 5/25/2021, but is now obsolete. It involves setting cluster numbers for each crab individually, rather than looping through and then naming only at the end

IMPORTANT: When running this, you will need to determine the number of clusters to split, and then designate the clusters as low/med/high expression on each day (shorthand: LMH)

### Crab A

```{r CrabA}
# Set path to output folder that will contain all heatmaps and text files for Crab A
crab_out <- paste0(heatmap_output, "Crab_A/")

# Read in data
crabdat <- read.delim(file = paste0(TPM_outpath,
                                    "Crab_", 
                                    crabs[1],
                                    "_TPMcts.txt"),
                   row.names = 1)

# Drop any genes with less than 5 counts across all 3 samples
keep <- rowSums(crabdat) >= 5
f.crabdat <- crabdat[keep, ]

print(paste("# of genes remaining after pre-filtering:", nrow(f.crabdat)))

print(paste("# of genes dropped:", nrow(crabdat) - nrow(f.crabdat), sep = " "))
  
# Construct heat map
pheatmap(f.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"))

pheatmap(f.crabdat, cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "row",
           main = "gene counts by day",
           fontsize = 9)
dev.off()


# SPECIFY NUMBER OF CLUSTERS
clusters <- 3


# Extract gene clusters
out <- pheatmap(f.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, clusters) # split dendrogram into the number of groups specified above

# Loop through gene clusters, printing heatmap for each
for (i in 1:clusters) {
  clust <- f.crabdat[which(lbl==i) %>% names(), ]
  # Write results to table
  write.table(clust, file = paste0(crab_out,
                                   "cluster_",
                                   i, ".txt"),
              sep = "\t",
              row.names = TRUE)
  # Create heat map for each cluster
  png(paste0(crab_out, "cluster_", i, "_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 ", i),
           fontsize = 9))
  dev.off()
}

# Look at cluster heatmaps, designate each day within clusters as low, medium, or high. Order corresponds to order of clusters (ex: if Cluster 2 is high/low/medium, the second object in the vector should be "HLM")
clustorder <- c("LLH", "MHL", "HML")

# See if any duplicates are present in our clusters
any(duplicated(clustorder))

# 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 cluster order
file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*\\.txt")),
            paste0(crab_out, "cluster_",
                   clustorder, ".txt"))
```

## Crab B

```{r CrabB}
# Set path to output folder that will contain all heatmaps and text files for Crab B
crab_out <- paste0(heatmap_output, "Crab_B/")

# Read in data
crabdat <- read.delim(file = paste0(TPM_outpath,
                                    "Crab_", 
                                    crabs[2],
                                    "_TPMcts.txt"),
                   row.names = 1)

# Drop any genes with less than 5 counts across all 3 samples
keep <- rowSums(crabdat) >= 5
f.crabdat <- crabdat[keep, ]

print(paste("# of genes remaining after pre-filtering:", nrow(f.crabdat)))

print(paste("# of genes dropped:", nrow(crabdat) - nrow(f.crabdat), sep = " "))
  
# Construct heat map
pheatmap(f.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"))

pheatmap(f.crabdat, cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "row",
           main = "gene counts by day",
           fontsize = 9)
dev.off()

# SPECIFY NUMBER OF CLUSTERS
clusters <- 5

# Extract gene clusters
out <- pheatmap(f.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, clusters) # split dendrogram into the number of groups specified above

# Loop through gene clusters, printing heatmap for each
for (i in 1:clusters) {
  clust <- f.crabdat[which(lbl==i) %>% names(), ]
  # Write results to table
  write.table(clust, file = paste0(crab_out,
                                   "cluster_",
                                   i, ".txt"),
              sep = "\t",
              row.names = TRUE)
  # Create heat map for each cluster
  png(paste0(crab_out, "cluster_", i, "_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 ", i),
           fontsize = 9))
  dev.off()
}

# Look at cluster heatmaps, designate each day within clusters as low, medium, or high. Order corresponds to order of clusters (ex: if Cluster 2 is high/low/medium, the second object in the vector should be "HLM")
clustorder <- c("LMH", "LHM", "HLH", "LHL", "HLL")

# See if any duplicates are present in our clusters
any(duplicated(clustorder))

# 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 cluster order
file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*\\.txt")),
            paste0(crab_out, "cluster_",
                   clustorder, ".txt"))
```


## Crab C

```{r CrabC}
# Set path to output folder that will contain all heatmaps and text files for Crab C
crab_out <- paste0(heatmap_output, "Crab_C/")

# Read in data
crabdat <- read.delim(file = paste0(TPM_outpath,
                                    "Crab_", 
                                    crabs[3],
                                    "_TPMcts.txt"),
                   row.names = 1)

# Drop any genes with less than 5 counts across all 3 samples
keep <- rowSums(crabdat) >= 5
f.crabdat <- crabdat[keep, ]

print(paste("# of genes remaining after pre-filtering:", nrow(f.crabdat)))

print(paste("# of genes dropped:", nrow(crabdat) - nrow(f.crabdat), sep = " "))
  
# Construct heat map
pheatmap(f.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"))

pheatmap(f.crabdat, cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "row",
           main = "gene counts by day",
           fontsize = 9)
dev.off()

# SPECIFY NUMBER OF CLUSTERS
clusters <- 8

# Extract gene clusters
out <- pheatmap(f.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, clusters) # split dendrogram into the number of groups specified above

# Loop through gene clusters, printing heatmap for each
for (i in 1:clusters) {
  clust <- f.crabdat[which(lbl==i) %>% names(), ]
  # Write results to table
  write.table(clust, file = paste0(crab_out,
                                   "cluster_",
                                   i, ".txt"),
              sep = "\t",
              row.names = TRUE)
  # Create heat map for each cluster
  png(paste0(crab_out, "cluster_", i, "_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 ", i),
           fontsize = 9))
  dev.off()
}

# Look at cluster heatmaps, designate each day within clusters as low, medium, or high. Order corresponds to order of clusters (ex: if Cluster 2 is high/low/medium, the second object in the vector should be "HLM")
clustorder <- c("LLH", "MLH", "HLM", "LHM", "MML", "HML", "LMH", "MHL")

# See if any duplicates are present in our clusters
any(duplicated(clustorder))

# 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 cluster order
file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*\\.txt")),
            paste0(crab_out, "cluster_",
                   clustorder, ".txt"))
```



## Crab E

```{r CrabE}
# Set path to output folder that will contain all heatmaps and text files for Crab E
crab_out <- paste0(heatmap_output, "Crab_E/")

# Read in data
crabdat <- read.delim(file = paste0(TPM_outpath,
                                    "Crab_", 
                                    crabs[4],
                                    "_TPMcts.txt"),
                   row.names = 1)

# Drop any genes with less than 5 counts across all 3 samples
keep <- rowSums(crabdat) >= 5
f.crabdat <- crabdat[keep, ]

print(paste("# of genes remaining after pre-filtering:", nrow(f.crabdat)))

print(paste("# of genes dropped:", nrow(crabdat) - nrow(f.crabdat), sep = " "))
  
# Construct heat map
pheatmap(f.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"))

pheatmap(f.crabdat, cluster_rows = TRUE,
           show_rownames = FALSE,
           na.rm = TRUE,
           cluster_cols = FALSE,
           scale = "row",
           main = "gene counts by day",
           fontsize = 9)
dev.off()

# SPECIFY NUMBER OF CLUSTERS
clusters <- 6

# Extract gene clusters
out <- pheatmap(f.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, clusters) # split dendrogram into the number of groups specified above

# Loop through gene clusters, printing heatmap for each
for (i in 1:clusters) {
  clust <- f.crabdat[which(lbl==i) %>% names(), ]
  # Write results to table
  write.table(clust, file = paste0(crab_out,
                                   "cluster_",
                                   i, ".txt"),
              sep = "\t",
              row.names = TRUE)
  # Create heat map for each cluster
  png(paste0(crab_out, "cluster_", i, "_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 ", i),
           fontsize = 9))
  dev.off()
}

# Look at cluster heatmaps, designate each day within clusters as low, medium, or high. Order corresponds to order of clusters (ex: if Cluster 2 is high/low/medium, the second object in the vector should be "HLM")
clustorder <- c("LMM", "MLH", "HLM", "HML", "MHL", "MHM")

# See if any duplicates are present in our clusters
any(duplicated(clustorder))

# 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 cluster order
file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*\\.txt")),
            paste0(crab_out, "cluster_",
                   clustorder, ".txt"))
```

## Crab G

```{r CrabG}
# Set path to output folder that will contain all heatmaps and text files for Crab G
crab_out <- paste0(heatmap_output, "Crab_G/")

# Read in data
crabdat <- read.delim(file = paste0(TPM_outpath,
                                    "Crab_", 
                                    crabs[5],
                                    "_TPMcts.txt"),
                   row.names = 1)

# Drop any genes with less than 5 counts across all 3 samples
keep <- rowSums(crabdat) >= 5
f.crabdat <- crabdat[keep, ]

print(paste("# of genes remaining after pre-filtering:", nrow(f.crabdat)))

print(paste("# of genes dropped:", nrow(crabdat) - nrow(f.crabdat), sep = " "))
  
# Construct heat map
pheatmap(log2(f.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"))

pheatmap(log2(f.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()

# SPECIFY NUMBER OF CLUSTERS
clusters <- 2

# Extract gene clusters
out <- pheatmap(log2(f.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, clusters) # split dendrogram into the number of groups specified above

# Loop through gene clusters, printing heatmap for each
for (i in 1:clusters) {
  clust <- f.crabdat[which(lbl==i) %>% names(), ]
  # Write results to table
  write.table(clust, file = paste0(crab_out,
                                   "cluster_",
                                   i, ".txt"),
              sep = "\t",
              row.names = TRUE)
  # Create heat map for each cluster
  png(paste0(crab_out, "cluster_", i, "_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 ", i),
           fontsize = 9))
  dev.off()
}

# Look at cluster heatmaps, designate each day within clusters as low, medium, or high. Order corresponds to order of clusters (ex: if Cluster 2 is high/low/medium, the second object in the vector should be "HLM")
clustorder <- c("LH", "HL")

# See if any duplicates are present in our clusters
any(duplicated(clustorder))

# 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 cluster order
file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*\\.txt")),
            paste0(crab_out, "cluster_",
                   clustorder, ".txt"))
```

## Crab H

```{r CrabH}
# Set path to output folder that will contain all heatmaps and text files for Crab H
crab_out <- paste0(heatmap_output, "Crab_H/")

# Read in data
crabdat <- read.delim(file = paste0(TPM_outpath,
                                    "Crab_", 
                                    crabs[6],
                                    "_TPMcts.txt"),
                   row.names = 1)

# Drop any genes with less than 5 counts across all 3 samples
keep <- rowSums(crabdat) >= 5
f.crabdat <- crabdat[keep, ]

print(paste("# of genes remaining after pre-filtering:", nrow(f.crabdat)))

print(paste("# of genes dropped:", nrow(crabdat) - nrow(f.crabdat), sep = " "))
  
# Construct heat map
pheatmap(log2(f.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"))

pheatmap(log2(f.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()

# SPECIFY NUMBER OF CLUSTERS
clusters <- 2

# Extract gene clusters
out <- pheatmap(log2(f.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, clusters) # split dendrogram into the number of groups specified above

# Loop through gene clusters, printing heatmap for each
for (i in 1:clusters) {
  clust <- f.crabdat[which(lbl==i) %>% names(), ]
  # Write results to table
  write.table(clust, file = paste0(crab_out,
                                   "cluster_",
                                   i, ".txt"),
              sep = "\t",
              row.names = TRUE)
  # Create heat map for each cluster
  png(paste0(crab_out, "cluster_", i, "_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 ", i),
           fontsize = 9))
  dev.off()
}

# Look at cluster heatmaps, designate each day within clusters as low, medium, or high. Order corresponds to order of clusters (ex: if Cluster 2 is high/low/medium, the second object in the vector should be "HLM")
clustorder <- c("LH", "HL")

# See if any duplicates are present in our clusters
any(duplicated(clustorder))

# 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 cluster order
file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*\\.txt")),
            paste0(crab_out, "cluster_",
                   clustorder, ".txt"))
```

## Crab I

```{r CrabI}
# Set path to output folder that will contain all heatmaps and text files for Crab I
crab_out <- paste0(heatmap_output, "Crab_I/")

# Read in data
crabdat <- read.delim(file = paste0(TPM_outpath,
                                    "Crab_", 
                                    crabs[7],
                                    "_TPMcts.txt"),
                   row.names = 1)

# Drop any genes with less than 5 counts across all 3 samples
keep <- rowSums(crabdat) >= 5
f.crabdat <- crabdat[keep, ]

print(paste("# of genes remaining after pre-filtering:", nrow(f.crabdat)))

print(paste("# of genes dropped:", nrow(crabdat) - nrow(f.crabdat), sep = " "))
  
# Construct heat map
pheatmap(log2(f.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"))

pheatmap(log2(f.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()

# SPECIFY NUMBER OF CLUSTERS
clusters <- 2

# Extract gene clusters
out <- pheatmap(log2(f.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, clusters) # split dendrogram into the number of groups specified above

# Loop through gene clusters, printing heatmap for each
for (i in 1:clusters) {
  clust <- f.crabdat[which(lbl==i) %>% names(), ]
  # Write results to table
  write.table(clust, file = paste0(crab_out,
                                   "cluster_",
                                   i, ".txt"),
              sep = "\t",
              row.names = TRUE)
  # Create heat map for each cluster
  png(paste0(crab_out, "cluster_", i, "_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 ", i),
           fontsize = 9))
  dev.off()
}

# Look at cluster heatmaps, designate each day within clusters as low, medium, or high. Order corresponds to order of clusters (ex: if Cluster 2 is high/low/medium, the second object in the vector should be "HLM")
clustorder <- c("LH", "HL")

# See if any duplicates are present in our clusters
any(duplicated(clustorder))

# 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 cluster order
file.rename(paste0(crab_out, list.files(crab_out, pattern = "^cluster_.*\\.txt")),
            paste0(crab_out, "cluster_",
                   clustorder, ".txt"))
```

