---
title: "04-Apul-RNA-sRNA-WGCNA"
author: "Kathleen Durkin"
date: "2024-08-26"
always_allow_html: true
output: 
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: true 
---

```{r setup, include=FALSE}
library(knitr)
knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = TRUE,         # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  comment = ""         # Prevents appending '##' to beginning of lines in code output
)
```

### Install and load packages

```{r load_libraries, inlcude = TRUE}
library(tidyverse)
library(ggplot2)
library(WGCNA)
library(magrittr)
library(genefilter)
```

Used Jennifer Chang's Bioinformatic Workflow [post on WGCNA analysis](https://bioinformaticsworkbook.org/tutorials/wgcna.html#gsc.tab=0) as a very helpful guide.

**NOTE: WGCNA is recommended for use with at *least* 15 samples to draw meaningful conclusions. We only have 5 samples from each species in our deep-dive dataset, so this code is primarily intended to be a trial run, developing an RNA-sRNA WGCNA coexpression pipeline in preparation for our *time series data*.**

## Load variance stabilized count data

Load in count matrices. The transformed RNA counts were generated using kallisto. The sRNA counts were generated using ShortStack. These counts were then normalized and variance stabilized using the DESeq2 package. For miRNAs, miRNAs were variance stabilized together with all sRNAs, then isolated based on ShortStack's miRNA annotation. See details on these transformations in `03-Apul-RN-summary` and `03.1-Apul-sRNA-summary`. 

```{r load-RNA-counts}
# RNA variance stabilized counts data
Apul_counts_RNA_vsd <- read.table("../output/03-Apul-RNA-summary/Apul_counts_RNA_variancestabilized.txt") 
Apul_counts_RNA_vsd_q75 <- read.table("../output/03-Apul-RNA-summary/Apul_counts_RNA_variancestabilized_q75.txt") 
Apul_counts_RNA_vsd_q95 <- read.table("../output/03-Apul-RNA-summary/Apul_counts_RNA_variancestabilized_q95.txt") 

# miRNA variance stabilized counts data
Apul_counts_miRNA_vsd <- read.table("../output/03.1-Apul-sRNA-summary/Apul_counts_miRNA_variancestabilized.txt") 
```



# Merge RNA and miRNA variance stabilized counts

Let's combine them into one dataset to feed into WGCNA.

```{r combine-counts-sets}
# Merge
Apul_counts_WGCNA_RNAq75 <- bind_rows(Apul_counts_miRNA_vsd_q75, Apul_counts_RNA_vsd)

# Convert from data frame to matrix
Apul_counts_WGCNA_RNAq75 <- as.matrix(Apul_counts_WGCNA_RNAq75)

# Transpose the normalized count data to meet WGCNA required input format
Apul_counts_WGCNA_RNAq75 = t(Apul_counts_WGCNA_RNAq75)
```


# WGCNA

Now we're ready to run WGCNA!

```{r WGCNA, cache=TRUE}
allowWGCNAThreads()          # allow multi-threading (optional)
# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to = 20, by = 2))

# Call the network topology analysis function
sft_Apul = pickSoftThreshold(
  Apul_counts_WGCNA_RNAq75,             # <= Input data
  #blockSize = 30,
  powerVector = powers,
  verbose = 5
  )

par(mfrow = c(1,2));
cex1 = 0.9;

plot(sft_Apul$fitIndices[, 1],
     -sign(sft_Apul$fitIndices[, 3]) * sft_Apul$fitIndices[, 2],
     xlab = "Soft Threshold (power)",
     ylab = "Scale Free Topology Model Fit, signed R^2",
     main = paste("Scale independence")
)
text(sft_Apul$fitIndices[, 1],
     -sign(sft_Apul$fitIndices[, 3]) * sft_Apul$fitIndices[, 2],
     labels = powers, cex = cex1, col = "red"
)
abline(h = 0.90, col = "red")
plot(sft_Apul$fitIndices[, 1],
     sft_Apul$fitIndices[, 5],
     xlab = "Soft Threshold (power)",
     ylab = "Mean Connectivity",
     type = "n",
     main = paste("Mean connectivity")
)
text(sft_Apul$fitIndices[, 1],
     sft_Apul$fitIndices[, 5],
     labels = powers,
     cex = cex1, col = "red")
```

I'm a little weirded out by the double curve in the Scale Independence plot, but this is a prelim test so I'm not going to worry about it right now. We want to pick a soft threshold power near the curve of the plot, so maybe 6, 7, 8, or 9? Let's pick 7 for now, but maybe experiment with other powers later.

```{r, cache=TRUE}
picked_power = 7
temp_cor <- cor       
cor <- WGCNA::cor         # Force it to use WGCNA cor function (fix a namespace conflict issue)
netwk_Apul <- blockwiseModules(Apul_counts_WGCNA_RNAq75,                # <= input here

                          # == Adjacency Function ==
                          power = picked_power,                # <= power here
                          networkType = "signed",

                          # == Tree and Block Options ==
                          deepSplit = 2,
                          pamRespectsDendro = F,
                          # detectCutHeight = 0.75,
                          minModuleSize = 30,
                          maxBlockSize = 4000,

                          # == Module Adjustments ==
                          reassignThreshold = 0,
                          mergeCutHeight = 0.25,

                          # == TOM == Archive the run results in TOM file (saves time)
                          saveTOMs = T,
                          saveTOMFileBase = "ER",

                          # == Output Options
                          numericLabels = T,
                          verbose = 3)


cor <- temp_cor     # Return cor function to original namespace

```

```{r move-WGCNA-output, engine='bash'}
# Check if there are any files starting with "ER" in the current directory
if ls ER-block* 1> /dev/null 2>&1; then
    # Move the files if they exist
    mv ER-block* ../output/04-Apul-RNA-sRNA-WGCNA
else
    echo "No files starting with 'ER' found."
fi
```

Take a look
```{r}
# Convert labels to colors for plotting
mergedColors = labels2colors(netwk_Apul$colors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(
  netwk_Apul$dendrograms[[1]],
  mergedColors[netwk_Apul$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE,
  hang = 0.03,
  addGuide = TRUE,
  guideHang = 0.05 )
```

```{r}
module_df_Apul <- data.frame(
  gene_id = names(netwk_Apul$colors),
  colors = labels2colors(netwk_Apul$colors)
)

module_df_Apul[1:5,]
```

```{r}
# Get Module Eigengenes per cluster
MEs0_Apul <- moduleEigengenes(Apul_counts_WGCNA_RNAq75, mergedColors)$eigengenes

# # Reorder modules so similar modules are next to each other
# MEs0_Apul <- orderMEs(MEs0_Apul)
# module_order_Apul = names(MEs0_Apul) %>% gsub("ME","", .)

# Add treatment names
MEs0_Apul$sample = row.names(MEs0_Apul)

# tidy & plot data
mME_Apul = MEs0_Apul %>%
  pivot_longer(-sample) %>%
  mutate(
    name = gsub("ME", "", name),
    # name = factor(name, levels = module_order)
  )

mME_Apul %>% ggplot(., aes(x=sample, y=name, fill=value)) +
  geom_tile() +
  theme_bw() +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    midpoint = 0,
    limit = c(-1,1)) +
  theme(axis.text.x = element_text(angle=90)) +
  labs(title = "Module-trait Relationships", y = "Modules", fill="corr")
```

```{r}
# Check which modules include miRNAs
module_df_Apul %>%
  filter(grepl("mir",gene_id)) %>%
  pull(colors) %>%
  unique()

module_df_Apul %>%
  filter(grepl("mir",gene_id)) %>%
  arrange(colors) %>%
  head(n=38)
```

```{r}
# pick out a few modules of interest. Let's do some of the modules that contain previously described miRNAs, since we have the best idea of their function
modules_of_interest = c("grey", "lightyellow")
# Define the colors corresponding to the modules of interest
module_colors <- c("grey" = "black", "lightyellow" = "yellow")

# Pull out list of genes in that module
submod_Apul = module_df_Apul %>%
  subset(colors %in% modules_of_interest)

row.names(module_df_Apul) = module_df_Apul$gene_id

subexpr_Apul = t(Apul_counts_WGCNA_RNAq75)[submod_Apul$gene_id,]

submod_df_Apul = data.frame(subexpr_Apul) %>%
  mutate(
    gene_id = row.names(.)
  ) %>%
  pivot_longer(-gene_id) %>%
  mutate(
    module = module_df_Apul[gene_id,]$colors
  )

submod_df_Apul %>% ggplot(., aes(x=name, y=value, group=gene_id)) +
  geom_line(aes(color = module),
            alpha = 0.2) +
  scale_color_manual(values = module_colors) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90)
  ) +
  facet_grid(rows = vars(module)) +
  labs(x = "treatment",
       y = "normalized expression")
```




