---
title: "Ptuh Gene expression lncRNA correlations"
output: 
  github_document: null
date: "2025-02-04"
editor_options: 
  chunk_output_type: console
---

This script conducts correlation network analyses for gene expression and lncRNA.  

# Set up 

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

Load libraries 
```{r, warning=FALSE, message=FALSE}
library(tidyverse)
library(ggplot2)
library(DESeq2)
library(igraph)
library(psych)
library(tidygraph)
library(ggraph)
library(WGCNA)
library(edgeR)
library(reshape2)
library(ggcorrplot)
library(corrplot)
```

# Pocillopora tuahiniensis 

## Load and format data 

```{r}
ptuh_mrna<-read_csv(file="F-Ptuh/output/06.2-Ptuh-Hisat/Ptuh-gene_count_matrix.csv")%>%
  rename("mrna_id"=gene_id, "sample47"=`RNA-POC-47`, "sample48"=`RNA-POC-48`, "sample50"=`RNA-POC-50`, "sample53"=`RNA-POC-53`, "sample57"=`RNA-POC-57`)

ptuh_lncrna<-read_table(file="F-Ptuh/output/08-Ptuh-lncRNA-matrix/Peve-lncRNA-counts.txt", skip=1)%>%
  select(!c("Chr", "Start", "End", "Strand", "Length"))%>%
  rename("lncrna_id"=Geneid, 
         "sample47"=`../output/06-Ptuh-Hisat/RNA-POC-47.sorted.bam`, 
         "sample48"=`../output/06-Ptuh-Hisat/RNA-POC-48.sorted.bam`, 
         "sample50"=`../output/06-Ptuh-Hisat/RNA-POC-50.sorted.bam`, 
         "sample53"=`../output/06-Ptuh-Hisat/RNA-POC-53.sorted.bam`, 
         "sample57"=`../output/06-Ptuh-Hisat/RNA-POC-57.sorted.bam`)

ptuh_mrna<-as.data.frame(ptuh_mrna)
ptuh_lncrna<-as.data.frame(ptuh_lncrna)

str(ptuh_mrna)
str(ptuh_lncrna)
```

Remove mRNA and lncRNA with 0 in all samples. 

```{r}
# Remove any genes with 0 counts across samples 
ptuh_mrna<-ptuh_mrna %>%
     mutate(Total = rowSums(.[, 2:6]))%>%
    filter(!Total==0)%>%
    dplyr::select(!Total)

ptuh_lncrna<-ptuh_lncrna %>%
     mutate(Total = rowSums(.[, 2:6]))%>%
    filter(!Total==0)%>%
    dplyr::select(!Total)
```

We now have 10382 lncRNAs and 26508 mRNAs.  

Set row names. 

```{r}
rownames(ptuh_mrna)<-ptuh_mrna$mrna_id
ptuh_mrna<-ptuh_mrna%>%select(!mrna_id)
```

```{r}
rownames(ptuh_lncrna)<-ptuh_lncrna$lncrna_id
ptuh_lncrna<-ptuh_lncrna%>%select(!lncrna_id)
```

Make sample metadata. 

```{r}
ptuh_sample_info <- data.frame(
  species = c("Pocillopora tuahiniensis"),
  sample = colnames(ptuh_mrna)
)

rownames(ptuh_sample_info)<-ptuh_sample_info$sample
```

## Filter data 

Filter data to only keep those expressed in all samples. 
```{r}
dge_mRNA <- DGEList(counts = ptuh_mrna)
dge_lncRNA <- DGEList(counts = ptuh_lncrna)

# Calculate CPM
cpm_mRNA <- cpm(dge_mRNA)
cpm_lncRNA <- cpm(dge_lncRNA)

# Filtering criteria: CPM > 1 in all samples (100%)
min_samples <- ncol(ptuh_mrna) * 1

# Apply the filter
keep_mRNA <- rowSums(cpm_mRNA > 1) >= min_samples
keep_lncRNA <- rowSums(cpm_lncRNA > 1) >= min_samples

# Filtered count matrices
filtered_mRNA_counts <- ptuh_mrna[keep_mRNA, ]
filtered_lncRNA_counts <- ptuh_lncrna[keep_lncRNA, ]

# Check the number of features retained
cat("mRNAs retained:", sum(keep_mRNA), "\n")
cat("lncRNAs retained:", sum(keep_lncRNA), "\n")
```

This reduced mRNAs to 14797 and lncRNAs to 2712

## Normalize data 

```{r}
dds_mrna <- DESeqDataSetFromMatrix(countData = filtered_mRNA_counts, 
                              colData = ptuh_sample_info, 
                              design = ~ 1)

# Variance Stabilizing Transformation
vsd_mrna <- assay(vst(dds_mrna, blind = TRUE))
```

```{r}
dds_lncrna <- DESeqDataSetFromMatrix(countData = filtered_lncRNA_counts, 
                              colData = ptuh_sample_info, 
                              design = ~ 1)

# Variance Stabilizing Transformation
vsd_lncrna <- assay(vst(dds_lncrna, blind = TRUE))
```

## Conduct module correlations with WGCNA 

Set soft threshold.  

```{r}
options(stringsAsFactors = FALSE)
enableWGCNAThreads()  # Enable multi-threading
allowWGCNAThreads(nThreads = 2)

# Combine mRNA and lncRNA datasets
datExpr <- t(rbind(vsd_mrna, vsd_lncrna))

sum(is.na(datExpr))  # Should be 0
sum(!is.finite(as.matrix(datExpr)))  # Should be 0

# Remove genes/samples with missing or infinite values
datExpr <- datExpr[complete.cases(datExpr), ]
datExpr <- datExpr[, colSums(is.na(datExpr)) == 0]

# # Choose a set of soft-thresholding powers
powers <- c(seq(from = 1, to=19, by=2), c(21:50)) #Create a string of numbers from 1 through 10, and even numbers from 10 through 20
# 
# # Call the network topology analysis function
sft <-pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

sft$powerEstimate
```

Plot the results.  
```{r}
sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# # # Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
      xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
 text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# # # this line corresponds to using an R^2 cut-off
 abline(h=0.9,col="red")
# # # Mean connectivity as a function of the soft-thresholding power
 plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
 text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
```

Selected power will be 30 - no value crossed threshold, selecting limit of tested values.   

### Generate network 

```{r}
selected_power<-30

# Network construction (adjust power based on sft output)
net = blockwiseModules(datExpr, power = selected_power,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = FALSE, verbose = 3)
```

View modules.  
```{r}
moduleEigengenes = moduleEigengenes(datExpr, colors = net$colors)$eigengenes

length(table(net$unmergedColors))
length(table(net$colors))

MEs<-net$MEs
moduleLabels<-net$colors
```

There are 22 modules after merging similar modules from original 306 modules. 

Determine whether mRNA and or lncRNA are present in each module. 

```{r}
# Get gene names and corresponding module colors
gene_module_info <- data.frame(
  Gene = colnames(datExpr),
  Module = moduleLabels
)

# Check structure
head(gene_module_info)

#Add ME to all the names
gene_module_info$Module <- paste0("ME", gene_module_info$Module)

```

Classify modules based on the proportion of the module comprised by mRNAs.   
 
```{r}
# Function to calculate the proportion of mRNAs (genes with "FUN" in ID)
calculate_mRNA_proportion <- function(genes) {
  total_genes <- length(genes)
  mRNA_count <- sum(grepl("gene", genes))
  
  # Proportion of mRNAs
  proportion_mRNA <- mRNA_count / total_genes
  return(proportion_mRNA)
}

# Apply the function to each module
module_mRNA_proportion <- tapply(gene_module_info$Gene, 
                                 gene_module_info$Module, 
                                 calculate_mRNA_proportion)

# View the proportions
module_mRNA_proportion

#Add ME to all the names
#names(module_mRNA_proportion) <- paste0("ME", names(module_mRNA_proportion))
```

### Run correlation between modules. 

```{r}
cor_matrix = cor(moduleEigengenes)
```

Compute correlations with Spearman correlation and BH p-value adjustment.   
```{r}
# Compute Spearman correlation between mRNA and lncRNA
ptuh_cor_results <- corr.test(moduleEigengenes, method = "spearman", adjust = "BH")

# Extract correlation values and p-values
ptuh_cor_matrix <- ptuh_cor_results$r  # Correlation coefficients
ptuh_p_matrix <- ptuh_cor_results$p  # Adjusted p-values
```

Construct network. 
```{r}
# Set correlation and significance thresholds
cor_threshold <- 0.6  # Adjust based on desired stringency
p_threshold <- 0.05

# Convert correlation matrix into an edge list
ptuh_significant_edges <- which(abs(ptuh_cor_matrix) > cor_threshold & ptuh_p_matrix < p_threshold, arr.ind = TRUE)

ptuh_edge_list <- data.frame(
  mRNA = rownames(ptuh_cor_matrix)[ptuh_significant_edges[,1]],
  lncRNA = colnames(ptuh_cor_matrix)[ptuh_significant_edges[,2]],
  correlation = ptuh_cor_matrix[ptuh_significant_edges]
)

# Construct network graph
ptuh_network <- graph_from_data_frame(ptuh_edge_list, directed = FALSE)

module_mRNA_proportion<-as.data.frame(module_mRNA_proportion)

V(ptuh_network)$prop_mrna <- module_mRNA_proportion$module_mRNA_proportion[match(V(ptuh_network)$name, rownames(module_mRNA_proportion))]
```

Plot network. 
```{r}
# Visualize network
ggraph(ptuh_network, layout = "fr") +  # Force-directed layout
  geom_edge_link(aes(edge_alpha = correlation), show.legend = TRUE, width=3) +
  geom_node_point(aes(colour=prop_mrna), size = 5) +
  scale_colour_gradient(name="Prop. mRNA", low = "purple", high = "cyan3")+
  geom_node_text(aes(label = name), repel = TRUE, size = 4) +
  theme_void() +
  labs(title = "P. tuahiniensis mRNA-lncRNA Network")
```

### Plot eigengene patterns and proportions of mRNA and lncRNAs

```{r}
module_mRNA_proportion$module_lncRNA_proportion<-1-module_mRNA_proportion$module_mRNA_proportion
```

View total size of modules 
```{r}
module_sizes <- table(moduleLabels)
module_sizes<-as.data.frame(module_sizes)
module_sizes$module<-paste0("ME", module_sizes$moduleLabels)
```

Plot a stacked bar plot. 

```{r}
stack_data<-module_mRNA_proportion
stack_data$module<-rownames(stack_data)
stack_data$size<-module_sizes$Freq[match(stack_data$module, module_sizes$module)]

stack_data$module <- factor(stack_data$module, 
                             levels = rev(stack_data$module[order(stack_data$size)]))

stack_data<-stack_data%>%
  mutate(mRNAs=module_mRNA_proportion*size)%>%
  mutate(lncRNAs=module_lncRNA_proportion*size)%>%
  select(!c(module_mRNA_proportion, module_lncRNA_proportion, size))


# Reshape the data for ggplot (long format)

stack_long <- melt(stack_data[, c("module", "mRNAs", "lncRNAs")], 
                    id.vars = "module",
                    variable.name = "RNA_Type", 
                    value.name = "Count")

ggplot(stack_long, aes(x = module, y = Count, fill = RNA_Type)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("mRNAs" = "skyblue", 
                               "lncRNAs" = "salmon"), 
                    labels = c("mRNA", "lncRNA"))+
  theme_classic() +
  labs(title = "P. tuahiniensis module mRNA and lncRNA components",
       x = "Module",
       y = "Counts",
       fill = "RNA Type") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

Next plot eigengene expression of each module across samples. 

```{r}
#convert wide format to long format for plotting  
head(moduleEigengenes)

plot_MEs <- moduleEigengenes 

plot_MEs$sample<-rownames(plot_MEs)

plot_MEs<-plot_MEs%>%
  pivot_longer(
    cols = where(is.numeric),  # Select only numeric columns
    names_to = "Module",       # Name for the new column containing the column names
    values_to = "Mean"         # Name for the new column containing the values
  )

expression_plots<-plot_MEs%>%
  group_by(Module) %>%
  
  ggplot(aes(x=sample, y=Mean)) +
  facet_wrap(~ Module)+
    geom_hline(yintercept = 0, linetype="dashed", color = "grey")+
  geom_point()+
  ggtitle("Pocillopora tuahiniensis mRNA-lncRNA modules")+
  geom_line(aes(group=1))+
  theme_classic()+ 
  theme(axis.text.x=element_text(angle = 90, hjust=1)); expression_plots
```

### Plot correlations of specific modules 

Which modules were significantly correlated? Show correlation matrix. 

```{r}
# Compute Spearman correlation between mRNA and lncRNA
ptuh_cor_results 

# Extract correlation values and p-values
ptuh_cor_matrix  # Correlation coefficients
ptuh_p_matrix  # Adjusted p-values
```

```{r}
ggcorrplot(ptuh_cor_results$r, 
           type = "lower", # Only plot lower triangle
           p.mat = ptuh_p_matrix, 
           sig.level = 0.05,  # Show significant correlations
           insig = "blank",  # Remove insignificant correlations
           lab = TRUE,  # Show correlation coefficients
           lab_size = 4,  # Label size
           colors = c("blue", "white", "red"),  # Color gradient
           title = "P. tuahiniensis Module Correlation Matrix")
```


## Next steps

Examine lncRNAs and mRNAs in co regulated networks 
MEs 

We could also do functional enrichment on each individual module to look at these different groups of coexpression.  

