---
title: "27-Apul-mRNA-miRNA-interactions-topGO"
author: "Kathleen Durkin"
date: "2025-04-28"
always_allow_html: true
output: 
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: true 
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
---

This script will use topGO to analyze functional enrichment of miRNA targets for Apul

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

library(dplyr)
library(tidyr)
library(ggplot2)
library(topGO)
library(tidyverse)
```

Code used below was created by `Jill Ashey`, modified for use with A.pulchra datasets by `Kathleen Durkin`

# Format topGO files

## Read in and format annotation files

```{r}
# Read in Apul annotations
annot_locations <- read.delim("https://raw.githubusercontent.com/urol-e5/deep-dive-expression/refs/heads/main/D-Apul/output/02-Apul-reference-annotation/Apulcra-genome-mRNA-IDmapping-2024_12_12.tab")
# Remove unneeded columns 
annot_locations <- annot_locations %>% dplyr::select(-X, -V13)
# Ensure there are no duplicate rows
annot_locations <- annot_locations %>% distinct()

head(annot_locations)
# Looks good!
```

This file shows each gene as it's genomic location. We want to use gene IDs to associate genes, so add gene IDs to this annotation table

Read in file that associates each mRNA genomic location with corresponding gene ID

```{r}
mRNA_FUNids <- read.table("../output/05-Apul-annotate-UTRs/Apul-mRNA-FUNids.txt", header=FALSE, col.names=c("location", "type", "mRNA_ID", "gene_ID", "product"), sep="\t")

# Remove unwanted text from parent column
mRNA_FUNids$gene_ID <- gsub("Parent=", "", mRNA_FUNids$gene_ID)
# Only need to keep mRNA location and gene ID
mRNA_FUNids <- mRNA_FUNids %>% dplyr::select(location, gene_ID)
```

join with annotation file

```{r}
# join
annot <- left_join(annot_locations, mRNA_FUNids, by = c("V1" = "location"))

# ensure there are no duplicate rows
annot <- annot %>% distinct()
```

## Set up gene2GO object

Want to isolate a list of GO terms per gene

```{r}
gene2go <- annot %>% filter(!is.na(Gene.Ontology.IDs)) %>% dplyr::select(gene_ID, Gene.Ontology.IDs)
gene2go <- gene2go %>% dplyr::rename(GO.ID = Gene.Ontology.IDs)

gene2go_list <- setNames(
  strsplit(as.character(gene2go$GO.ID), ";"), 
  gene2go$gene_ID
)
```

Note: I think this means genes that had a Uniprot ID but no GO terms are excluded from this analysis

## Define reference set

Define reference set of genes. This should be all genes *found in our samples*, NOT all genes in the A.pulchra genome. Some genes (e.g., reproduction pathways) may not be found/expected in our samples for valid biological reasons.

```{r}
# Read in counts matrix
Apul_counts <- read.csv("../output/02.20-D-Apul-RNAseq-alignment-HiSat2/apul-gene_count_matrix.csv")
# Exclude genes with all 0 counts
Apul_counts <- Apul_counts[rowSums(Apul_counts[, 2:6]) != 0, ]

# Select gene IDs of the genes present in our samples
all_genes <- Apul_counts$gene_id
length(all_genes)
```

So we have a reference set of 30094 genes present in our samples.

## Read in PCC/miranda data

This is a table of all putative miRNA-mRNA binding predicted by miRanda, plus Pearsons correlation coefficients for coexpression of each putative binding pair. It only includes interactions with significant PCC (pval < 0.05)

```{r}
data <- read.csv("../output/14.1-Apul-miRNA-mRNA-coexpression-additional_inputs/miRanda-PCC-significant-mRNA_3UTR_5UTR.csv") %>% filter(!is.na(mirna))
nrow(data)
head(data)
```

For the purposes of functional annotation and enrichment, I don't care if an miRNA may bind to several locations of a given gene, and I don't want those duplicates being overrepresented during enrichment analysis.
```{r}
# Simplify target genes to gene IDs only (no binding coordinates)
data <- data %>% dplyr::select(-Target, -Score, -Energy_Kcal_Mol, -Subject_Aln, -Query_Aln, -Al_Len, -Subject_Identity, -Query_Identity) %>%
  distinct()
nrow(data)
```

# FA of all miRNA targets

Functional annotation of all putative miRNA targets

```{r}
cor_bind_FA <- left_join(data, annot, by = c("mRNA" = "gene_ID")) %>% distinct()

nrow(cor_bind_FA)
nrow(cor_bind_FA[!is.na(cor_bind_FA$Gene.Ontology.IDs),])
```

Of the 71634 putative miRNA targets predicted by miRanda with significant PCC, 20294 have available annotations

```{r}
high0.5_cor_bind_FA <- cor_bind_FA[abs(cor_bind_FA$PCC.cor) > 0.5,]

nrow(high0.5_cor_bind_FA)
nrow(high0.5_cor_bind_FA[!is.na(high0.5_cor_bind_FA$Gene.Ontology.IDs),])
```

Of the 9480 putative miRNA targets predicted by miRanda that that have highly correlated expression (magnitude of correlation > 0.5), 2711 have available annotations.

```{r}
high0.6_cor_bind_FA <- cor_bind_FA[abs(cor_bind_FA$PCC.cor) > 0.6,]

nrow(high0.6_cor_bind_FA)
nrow(high0.6_cor_bind_FA[!is.na(high0.6_cor_bind_FA$Gene.Ontology.IDs),])
```

2090 have correlation of at least 0.6, and of those 601 are annotated

```{r}
high0.7_cor_bind_FA <- cor_bind_FA[abs(cor_bind_FA$PCC.cor) > 0.7,]

nrow(high0.7_cor_bind_FA)
nrow(high0.7_cor_bind_FA[!is.na(high0.7_cor_bind_FA$Gene.Ontology.IDs),])
```

314 have correlation of at least 0.7, and of those 86 are annotated


Save
```{r}
write.csv(cor_bind_FA, "../output/27-Apul-mRNA-miRNA-interactions-topGO/miRNA_targets_FA.csv")
write.csv(high0.5_cor_bind_FA, "../output/27-Apul-mRNA-miRNA-interactions-topGO/miRNA_high0.5_cor_targets_FA.csv")
```

# FE of specific miRNA's targets (all targets)

Create topGO function for use with miRNA names

```{r}

miRNA_topGO_FE <- function(miRNA.name, input_interactions) {

  #Isolate genes in our input module of interest
  interacting_genes <- input_interactions %>%
    filter(miRNA == miRNA.name) %>%
    pull(mRNA)
  
  if (length(interacting_genes) > 0) {
    # Create factor for all reference genes, where 1 represents module membership and 0 means the gene is not in module of interest
    gene_list <- factor(as.integer(all_genes %in% interacting_genes))
    names(gene_list) <- all_genes
    str(gene_list)

    ## Biological Process ##
    # Create topGO object
    GO_BP <- new("topGOdata", 
                ontology = "BP", # Biological Process
                allGenes = gene_list,
                annot = annFUN.gene2GO, 
                gene2GO = gene2go_list,
                geneSel=topDiffGenes)
    
    # Run GO enrichment test
    GO_BP_FE <- runTest(GO_BP, algorithm = "weight01", statistic = "fisher")
    # View the results
    GO_BP_results <- GenTable(GO_BP, Fisher = GO_BP_FE, orderBy = "Fisher",  topNodes = 100, numChar = 51)
    # Filter by significant results
    GO_BP_results$Fisher<-as.numeric(GO_BP_results$Fisher)
    GO_BP_results_sig<-GO_BP_results[GO_BP_results$Fisher<0.05,]
    
    
    ## Molecular Function ##
    # Create topGO object
    GO_MF <- new("topGOdata", 
                ontology = "MF", # Molecular Function
                allGenes = gene_list,
                annot = annFUN.gene2GO, 
                gene2GO = gene2go_list,
                geneSel=topDiffGenes)
    
    # Run GO enrichment test
    GO_MF_FE <- runTest(GO_MF, algorithm = "weight01", statistic = "fisher")
    # View the results
    GO_MF_results <- GenTable(GO_MF, Fisher = GO_MF_FE, orderBy = "Fisher",  topNodes = 100, numChar = 51)
    # Filter by significant results
    GO_MF_results$Fisher<-as.numeric(GO_MF_results$Fisher)
    GO_MF_results_sig<-GO_MF_results[GO_MF_results$Fisher<0.05,]
  
    # Return
    # Add type column only if results exist
    if (nrow(GO_BP_results_sig) > 0) {
      GO_BP_results_sig$type <- "Biological.Process"
    }
    if (nrow(GO_MF_results_sig) > 0) {
      GO_MF_results_sig$type <- "Molecular.Function"
    }
    GO_results <- rbind(GO_BP_results_sig, GO_MF_results_sig)
    print(GO_results)
  }
}

miRNA_topGO_FE("Cluster_9532", cor_bind_FA)
```

Loop through all miRNA and run functional enrichment on the miRNA's targets (all predicted targets)

```{r, cache=TRUE}

interacting_miRNAs <- unique(cor_bind_FA$miRNA) %>% na.omit
results_all_targets <- NULL  # initialize empty df

for(miRNA in interacting_miRNAs) {
  
  # Run topGO enrichment function
  miRNA_results <- miRNA_topGO_FE(miRNA, cor_bind_FA)
  
  # Only keep results if not empty
  if (nrow(miRNA_results) > 0) {
    
    # Add the miRNA source column
    miRNA_results$miRNA <- miRNA

    # Bind to the accumulating results data frame
    results_all_targets <- rbind(results_all_targets, miRNA_results)
  }
}

head(results_all_targets)

```

Save results

```{r}
write.csv(results_all_targets, "../output/27-Apul-mRNA-miRNA-interactions-topGO/miRNA_all_targets_topGO_FE.csv")
```

# FE of specific miRNA's targets (high 0.5 cor targets)

Loop through all miRNA and run functional enrichment on the miRNA's highly correlated targets (PCC magnitude > 0.5)

```{r, cache=TRUE}

interacting_miRNAs_high0.5 <- unique(high0.5_cor_bind_FA$miRNA) %>% na.omit
results_high0.5_cor_targets <- NULL  # initialize empty df

for(miRNA in interacting_miRNAs_high0.5) {
  
  # Run topGO enrichment function
  miRNA_results <- miRNA_topGO_FE(miRNA, high0.5_cor_bind_FA)
  
  # Only keep results if not empty
  if (nrow(miRNA_results) > 0) {
    
    # Add the miRNA source column
    miRNA_results$miRNA <- miRNA

    # Bind to the accumulating results data frame
    results_high0.5_cor_targets <- rbind(results_high0.5_cor_targets, miRNA_results)
  }
}

head(results_high0.5_cor_targets)

```

Save results

```{r}
write.csv(results_high0.5_cor_targets, "../output/27-Apul-mRNA-miRNA-interactions-topGO/miRNA_high0.5_cor_targets_topGO_FE.csv")
```


# FE of specific miRNA's targets (high 0.6 cor targets)

Loop through all miRNA and run functional enrichment on the miRNA's highly correlated targets (PCC magnitude > 0.6)

```{r, cache=TRUE}

interacting_miRNAs_high0.6 <- unique(high0.6_cor_bind_FA$miRNA) %>% na.omit
results_high0.6_cor_targets <- NULL  # initialize empty df

for(miRNA in interacting_miRNAs_high0.6) {
  
  # Run topGO enrichment function
  miRNA_results <- miRNA_topGO_FE(miRNA, high0.6_cor_bind_FA)
  
  # Only keep results if not empty
  if (nrow(miRNA_results) > 0) {
    
    # Add the miRNA source column
    miRNA_results$miRNA <- miRNA

    # Bind to the accumulating results data frame
    results_high0.6_cor_targets <- rbind(results_high0.6_cor_targets, miRNA_results)
  }
}

head(results_high0.6_cor_targets)

```

Save results

```{r}
write.csv(results_high0.6_cor_targets, "../output/27-Apul-mRNA-miRNA-interactions-topGO/miRNA_high0.6_cor_targets_topGO_FE.csv")
```


# FE of specific miRNA's targets (high 0.7 cor targets)

Loop through all miRNA and run functional enrichment on the miRNA's highly correlated targets (PCC magnitude > 0.7)

```{r, cache=TRUE}

interacting_miRNAs_high0.7 <- unique(high0.7_cor_bind_FA$miRNA) %>% na.omit
results_high0.7_cor_targets <- NULL  # initialize empty df

for(miRNA in interacting_miRNAs_high0.7) {
  
  # Run topGO enrichment function
  miRNA_results <- miRNA_topGO_FE(miRNA, high0.7_cor_bind_FA)
  
  # Only keep results if not empty
  if (nrow(miRNA_results) > 0) {
    
    # Add the miRNA source column
    miRNA_results$miRNA <- miRNA

    # Bind to the accumulating results data frame
    results_high0.7_cor_targets <- rbind(results_high0.7_cor_targets, miRNA_results)
  }
}

head(results_high0.7_cor_targets)

```

Save results

```{r}
write.csv(results_high0.7_cor_targets, "../output/27-Apul-mRNA-miRNA-interactions-topGO/miRNA_high0.7_cor_targets_topGO_FE.csv")
```

