---
title: "12-piRNA-comparizon"
author: "Javier Rodriguez-Casariego"
date: "`r Sys.Date()`"
output: html
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, cache = TRUE)
knitr::opts_knit$set(root.dir = '/Users/jarcasariego/Dropbox/Other_projects/E5_work/deep-dive/DEF-cross-species/')
```

```{r}
#load libraries
library(Biostrings) # Biostrings_2.62.0 
require(ggplot2)    # ggplot2_3.3.6
library(ggpubr)
library(ggseqlogo)  # install.packages("ggseqlogo")
library(reshape2)   # reshape2_1.4.4
library(data.table)
library(gghighlight) #install.packages("gghighlight")
library(plotrix)
library(dplyr)
library(tidyr)
library(tibble)
library(tidyverse)

```

# putative piRNA reads size distribution
```{r}
my.outprefix <- paste(Sys.Date(),"all_species",sep="_")

#load list of file names
APUL_L = list.files(path = "../D-Apul/output/18-Apul-piRNA/0_piRNA_pipeline_proTRAC/piRNA_length_data", pattern=".txt", recursive = T, full.names = T)

PMEA_L = list.files(path = "../F-Pmea/output/18-Pmea-piRNA/0_piRNA_pipeline_proTRAC/piRNA_length_data", pattern=".txt", recursive = T, full.names = T)

PEVE_L = list.files(path = "../E-Peve/output/18-Peve-piRNA/0_piRNA_pipeline_proTRAC/piRNA_length_data", pattern=".txt", recursive = T, full.names = T)

L = c(APUL_L,PMEA_L,PEVE_L)
names <- L %>%
  substring(.,75) %>%
  gsub("_piRNA_length.txt", "", .)

DFs = lapply(L, function(x) {
  data <- read.table(paste0(x), comment.char = "")
  colnames(data) <- c("freq", "length")
  return(data)
})

combo = Reduce(function(...) merge(..., by = "length", all=T), DFs)
colnames(combo) <- c("length", names)
combo <- combo %>%
  mutate(ACR = rowSums(combo[,2:6])) %>%
  mutate(POC = rowSums(combo[,7:11])) %>%
  mutate(POR = rowSums(combo[,12:14])) %>%
  select(length,ACR,POC,POR) %>%
  pivot_longer(c("ACR", "POC", "POR"), 
  values_to = "sums", names_to = "group")

group.labs <- c("A. pulchra", "P. tuahiniensis", "P. evermanni")
names(group.labs) <- c("ACR","POC","POR")

length_plot <- ggplot(combo) + geom_bar(aes(x=length, y=sums, fill=group), stat="identity") +
 scale_x_continuous(name = "Length (bp)", breaks = c(24,25,26,27,28,29,30,31,32)) +
  scale_fill_manual(values = c("#408EC6", "#7A2048", "#1E2761")) +
  theme_bw() +
  facet_wrap(~group, ncol = 1, labeller = labeller(group = group.labs), strip.position="left") +
  theme(axis.ticks.y.left = element_blank(),
        axis.text.y.left = element_blank(),
        axis.title.y = element_blank(),
        axis.line.x = element_line(),
        panel.grid = element_blank(), 
        panel.border = element_blank(),
        strip.background = element_blank(),
        legend.position = 'none',
        strip.text = element_text(face = "italic", size = 9))

length_plot

ggsave(filename = paste0("output/12-piwiRNA_comparizon/", my.outprefix,"piRNA_length_plot.png"), length_plot, width = 3, height = 3)

```

# putative piRNA ping pong signatures
```{r}
#Graph comparing ping pong per million bootstrapped reads across species for PPmeter output with all mapped putative piRNAs

my.outprefix <- paste(Sys.Date(),"ppr_mbr_analysis",sep="_")

#load list of file names
APUL_L = list.files(path = "../D-Apul/output/18-Apul-piRNA/1_PPmeter/ppr_Apul", pattern=".txt", recursive = T, full.names = T)

PMEA_L = list.files(path = "../F-Pmea/output/18-Pmea-piRNA/1_PPmeter/ppr_Pmea", pattern=".txt", recursive = T, full.names = T)

PEVE_L = list.files(path = "../E-Peve/output/18-Peve-piRNA/1_PPmeter/ppr_Peve", pattern=".txt", recursive = T, full.names = T)

L = c(APUL_L,PMEA_L,PEVE_L)
names <- L %>%
  substring(.,51) %>%
  gsub("_ppr.txt", "", .)

#this function keeps only the instances of 10 bp overlaps per million read pairs
DFs = lapply(L, function(x) {
  data <- read.table(x, skip = 2, nrows = 100)
  data <- data[,c(1,34)]
  colnames(data) <- c("rep", "score")
  score <- median(data$score)
  return(score)
})
L[[1]]
data <- read.table("../D-Apul/output/18-Apul-piRNA/1_PPmeter/ppr_Apul/sRNA-ACR-140_ppr.txt", skip = 104)
data <- data[,c(1,33)]

#run for loop to get median ppm per replicate
med_ppm <- data.frame()
for(i in 1:length(L)) {
  total <- cbind(names[i], DFs[[i]])
  med_ppm <- rbind(med_ppm, total)
}
med_ppm <- med_ppm %>%
  mutate(group = substr(gsub("sRNA-", "", med_ppm$V1), 1, 3))%>%
  select(group, V2)

colnames(med_ppm) <- c("group", "score")
med_ppm$score <- as.numeric(med_ppm$score)

#divide ppr-mbr score by 10e5 for plotting
med_ppm$score <- med_ppm$score/100000

#order for plotting
med_ppm$group <- factor(med_ppm$group , levels=c("ACR", "POC", "POR"))

#get specific p-values for plotting

####### Tests
apul.peve <- wilcox.test (med_ppm$score[med_ppm$group == "ACR"],med_ppm$score[med_ppm$group == "POR"])
apul.pmea <- wilcox.test (med_ppm$score[med_ppm$group == "ACR"],med_ppm$score[med_ppm$group == "POC"])
pmea.peve <- wilcox.test (med_ppm$score[med_ppm$group == "POC"],med_ppm$score[med_ppm$group == "POR"])


my_comparisons <- list( c("POC", "POR"), c("ACR", "POC"), c("ACR", "POR"))
zscore_plt <- ggboxplot(med_ppm, x = "group", y = "score",
          color = "group", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons, method = "wilcox.test",label.y = c(2.5, 3.4, 3.7)) +
  scale_y_continuous(name = "median x10e5 ppr-mbr (100 bootstraps)") +
  theme(axis.title.x = element_blank(),
        legend.position = "none")
zscore_plt

ggsave(filename = paste0("output/12-piwiRNA_comparizon/", my.outprefix,"ppm_all_species.png"), zscore_plt, width = 5, height = 4)

#dev.off()

# plot ping pong overlap per nucleotide position 

DF_nuc = lapply(L, function(x) {
  data <- read.table(x, skip = 104)
  data <- data[1,c(1:30)] %>%
    transpose()
})

nuc_score <- data.frame()
for(i in 1:length(L)) {
  total <- cbind(names[i], DF_nuc[[i]])
  nuc_score <- rbind(nuc_score, total)
}

nuc_score <- nuc_score %>%
  mutate(group = substr(gsub("sRNA-", "", nuc_score$`names[i]`), 1, 3)) %>%
  mutate(score = nuc_score$V1) %>%
  mutate(position = rep(c(1:30),16)) %>%
  select(position, score, group)

overlap_plot <- ggplot(nuc_score, aes(fill = group)) + geom_bar(aes(x=position, y=score), stat="identity") +
 scale_x_continuous(name = "5' overlap (bp)", limits = c(1,23)) +
 scale_fill_manual(values = c("#408EC6", "#7A2048", "#1E2761")) +
  facet_wrap(~group, ncol = 1) +
  gghighlight(position == 10, calculate_per_facet = TRUE, unhighlighted_params = list(fill = NULL, alpha = 0.3 )) +
  theme_bw() +
  theme(axis.ticks.y.left = element_blank(),
        axis.text.y.left = element_blank(),
        axis.title.y = element_blank(),
        panel.grid = element_blank(), 
        panel.border = element_blank(), 
        axis.line.x = element_line(),
        strip.text.x = element_blank(),
        legend.position = 'none',
        )

overlap_plot

ggsave(filename = paste0("output/12-piwiRNA_comparizon/", my.outprefix,"piRNA_pp_overlap_plot.png"), overlap_plot, width = 3, height = 3)




```

# putative piRNA base composition (logo graphs)
```{r}
library(gridExtra)

# Nucleotide Distribution Plots

my.outprefix <- paste(Sys.Date(),"all_species",sep="_")

#load trimmed fasta sequences without headers
APUL_sequences <- read.table("../D-Apul/output/18-Apul-piRNA/0_piRNA_pipeline_proTRAC/total_fasta_APUL_replaced.fasta")
PMEA_sequences <- read.table("../F-Pmea/output/18-Pmea-piRNA/0_piRNA_pipeline_proTRAC/total_fasta_PMEA_replaced.fasta")
PEVE_sequences <- read.table("../E-Peve//output/18-Peve-piRNA/0_piRNA_pipeline_proTRAC/total_fasta_PEVE_replaced.fasta")

#make position weight matrix
pfm_APUL <- consensusMatrix(APUL_sequences$V1)
pfm_PMEA <- consensusMatrix(PMEA_sequences$V1)
pfm_PEVE <- consensusMatrix(PEVE_sequences$V1)

#remove large modified fasta file
rm(APUL_sequences, PMEA_sequences, PEVE_sequences)  

# create a list to facet the plot
seqs <- list(pfm_APUL,pfm_PMEA, pfm_PEVE)
names(seqs) <- c("APUL", "PTHU","PEVE")

#make logo plot
logo_plot <- ggplot() + 
  geom_logo(seqs, method = 'bit', seq_type = "rna",) + 
  theme_logo() +
  facet_grid(vars(seq_group), switch = "y") +
  theme_bw() + 
  scale_x_continuous(name = "position (bp)", limits = c(0.5,20)) +
  theme(axis.ticks.y.left = element_blank(),
        axis.text.y.left = element_blank(),
        axis.title.y = element_blank(),
        panel.grid = element_blank(), 
        panel.border = element_blank(), 
        axis.line.x = element_line(),
        strip.text = element_blank(),
        strip.background = element_blank()
        )
logo_plot


ggsave(filename = paste0("output/12-piwiRNA_comparizon/", my.outprefix,"piRNA_nucleotide_plot.png"), logo_plot, width = 3, height = 3)
```

# Ping-pong z scores for putative piRNAs
```{r}
# Graph z-scores per species

APUL_Z <- read.table("output/12-piwiRNA_comparizon/pp_z-scores/APUL_pp_z-scores.txt") %>%
  dplyr::mutate(sp = "ACR",
                z_score = V2) %>%
  select(sp,z_score) 
  
PMEA_Z <- read.table("output/12-piwiRNA_comparizon/pp_z-scores/PMEA_pp_z-scores.txt") %>%
  dplyr::mutate(sp = "POC",
                z_score = V2) %>%
  select(sp,z_score) 
PEVE_Z <- read.table("output/12-piwiRNA_comparizon/pp_z-scores/PEVE_pp_z-scores.txt") %>%
  dplyr::mutate(sp = "POR",
                z_score = V2) %>%
  select(sp,z_score) 

z_group <- rbind(PMEA_Z, APUL_Z, PEVE_Z) %>%
  dplyr::group_by(sp) %>%
  summarise_each(funs(mean,sd,std.error)) %>%
  mutate(mean = round(mean, 1))
z_group$sp <- factor(z_group$sp, levels = c("ACR", "POC", "POR"), ordered = TRUE)
stat_res <- data_frame(y = c(rep("mean",3)), group1 = c("ACR","ACR","POC"), group2 = c("POC", "POR", "POR"), p.value = c("*** (Z=27.36)","*** (Z=22.24)","*** (Z=4.93)"))

z_score_plot <- ggplot(z_group) +
    geom_bar( aes(x=sp, y=mean, fill=sp), stat="identity", width=0.3) +
    geom_errorbar( aes(x=sp, ymin=mean-std.error, ymax=mean+std.error), width=0.1) +
  coord_flip() +
  theme_bw() +
  scale_x_discrete(limits=rev) +
  scale_fill_manual(values = c("#408EC6", "#7A2048", "#1E2761")) +
  scale_y_continuous(name = "Z-score", limits = c(0,100)) +
  theme(panel.grid = element_blank(), 
        panel.border = element_blank(), 
        axis.line.x = element_line(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(), 
        legend.position = 'none'
        ) +
  geom_hline(yintercept = 0, size = .75) +
  stat_pvalue_manual(data = stat_res, 
                     label = "p.value",
                     y.position = c(51, 88, 81), coord.flip = TRUE, label.size = 2, tip.length = 0.01);z_score_plot


```



```{r}

library(ggpubr)
library(cowplot)
final_plt <- plot_grid(length_plot,overlap_plot, logo_plot, z_score_plot, nrow = 1, align = "h", labels = c("A", "B", "C", "D"), label_x = 0, label_y = 1.15, hjust = -0.1, vjust = 5)
final_plt

ggsave(filename = "figures/piRNA_final_plot.png", final_plt, width = 12, height = 3)
```

# Overlaps

```{r}

library(ggtext) 
library(ggExtra)
library(gtable)
overlap <- read.csv("output/12-piwiRNA_comparizon/piRNA_feature_overlap.csv")
overlap$element <- factor(overlap$element, levels = c("TE", "gene","intergenic"), ordered = TRUE)

overlap <- overlap %>% 
  group_by(sp) %>% # Variable to be transformed
  mutate(perc = (count / sum(count)*100)) %>% 
  ungroup()
  
cp <- coord_polar(theta = "y", start = 0)
cp$is_free <- function() TRUE

over_plot <- ggplot(overlap, aes(x="", y=perc, alpha=element, fill = sp)) +
  scale_fill_manual(values = c("#408EC6", "#7A2048", "#1E2761"), labels = c("*A. pulchra*", "*P. tuahiniensis*", "*P. evermanni*")) +
  scale_alpha_discrete(range = c(.4,1), labels = c("Repetitive elements-derived piRNAs", "Gene-derived piRNAs", "Intergenic piRNAs")) +
  geom_bar(stat="identity", width=1) +
  cp +
  facet_wrap(~sp, scales = "free") +
  geom_text(aes(label = ifelse(perc > 15 , paste0(round(perc), "%"), "")), 
            position = position_stack(vjust = 0.5),
            show.legend = FALSE) +
  geom_text(aes(x = 1.7, label = ifelse(perc < 15 , paste0(round(perc), "%"), "")), 
            position = position_stack(vjust = 0.5),
            show.legend = FALSE) +
  theme_void() +
  theme(aspect.ratio = 1, 
        strip.text = element_blank(), 
        legend.title = element_blank(),
        legend.text = element_markdown(),
        legend.position = "bottom",
        legend.box = "vertical");over_plot;
over_plot

ggsave(filename = "figures/piRNA_overlap_plot.png", over_plot)


```

# clustered piRNAs characterization 
## Cluster piRNA composition graph
```{r}
APUL_clust <- read.table("output/12-piwiRNA_comparizon/clusters/APUL_cluster_summary.tab", header = T) %>%
  dplyr::mutate(sp = "ACR") 
PMEA_clust <- read.table("output/12-piwiRNA_comparizon/clusters/PMEA_cluster_summary.tab", header = T) %>%
  dplyr::mutate(sp = "POC") 
PEVE_clust <- read.table("output/12-piwiRNA_comparizon/clusters/PEVE_cluster_summary.tab", header = T) %>%
  dplyr::mutate(sp = "POR") 

clust_summary <- rbind(APUL_clust, PMEA_clust, PEVE_clust) %>%
  select(., sp, perc_unique_tags_in_clusters, perc_reads_in_clusters, avg_perc_1T_in_clusters) %>%
  group_by(sp) %>%
  summarise_all(list(mean)) %>%
  gather(key = vars, value = percent, !sp) %>%
  mutate(vars = factor(vars, level=c('avg_perc_1T_in_clusters','perc_unique_tags_in_clusters', 'perc_reads_in_clusters'), ordered = T))

clust_cov_plot <- ggplot(clust_summary, aes(x=sp, y=percent, fill = sp, alpha = vars)) +
  geom_col(position = "dodge") +
  coord_flip() +
  theme_classic() +
  scale_x_discrete(guide = 'none', limits=rev) +
  scale_fill_manual(values = c("#408EC6", "#7A2048", "#1E2761")) +
  scale_alpha_discrete(range = c(1,.4), labels = c('Average percent 1T in clusters', 'Average percent of unique sequences in clusters', 'Average percent of total reads in clusters')) + 
  ylab("Fraction (%)") +
  theme(legend.position = "right",
        legend.title = element_blank(),
        panel.grid = element_blank(), 
        panel.border = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.y = element_blank(), 
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        );clust_cov_plot  
ggsave(filename = "figures/cluster_charact_plot.png", clust_cov_plot, width = 6, height = 3)
```


## Cluster overlaps

```{r}

APUL_m_clust <- read.table("../D-Apul/output/18-Apul-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_APUL/APUL.merged.clusters.bed", header = F) %>%
  mutate(bp_overlap = V3-V2) %>%
  select(V1,V2,V3,bp_overlap) %>%
  mutate(sp = "ACR")
colnames(APUL_m_clust) <- c("chr","start","end","bp_overlap","sp")
PMEA_m_clust <- read.table("../F-Pmea/output/18-Pmea-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_PMEA/PMEA.merged.clusters.bed", header = F) %>%
  mutate(bp_overlap = V3-V2) %>%
  select(V1,V2,V3,bp_overlap) %>%
  mutate(sp = "POC")
colnames(PMEA_m_clust) <- c("chr","start","end","bp_overlap","sp")
PEVE_m_clust <- read.table("../E-Peve/output/18-Peve-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_PEVE/PEVE.merged.clusters.bed", header = F) %>%
  mutate(bp_overlap = V3-V2) %>%
  select(V1,V2,V3,bp_overlap) %>%
  mutate(sp = "POR")
colnames(PEVE_m_clust) <- c("chr","start","end","bp_overlap","sp")

APUL_gene <- read.table("../D-Apul/output/18-Apul-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_APUL/APUL.merged.clusters.gene.overlap.tab", header = F) %>%
  select(V1, V2, V3, V13, ) %>%
  mutate(sp = "ACR")
colnames(APUL_gene) <- c("chr","start","end","bp_overlap","sp")
PMEA_gene <- read.table("../F-Pmea/output/18-Pmea-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_PMEA/PMEA.merged.clusters.gene.overlap.tab", header = F) %>%
  select(V1, V2, V3, V13) %>%
  mutate(sp = "POC")
colnames(PMEA_gene) <- c("chr","start","end","bp_overlap","sp")
PEVE_gene <- read.table("../E-Peve/output/18-Peve-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_PEVE/PEVE.merged.clusters.gene.overlap.tab", header = F) %>%
  select(V1, V2, V3, V13) %>%
  mutate(sp = "POR")
colnames(PEVE_gene) <- c("chr","start","end","bp_overlap","sp")

APUL_TE <- read.table("../D-Apul/output/18-Apul-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_APUL/APUL_piRNA_clusters_TEs_intersect.txt", header = F) %>%
  select(V1, V2, V3, V9) %>%
  mutate(sp = "ACR")
colnames(APUL_TE) <- c("chr","start","end","bp_overlap","sp")
PMEA_TE <- read.table("../F-Pmea/output/18-Pmea-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_PMEA/PMEA_piRNA_cluster_TEs_intersect.txt", header = F) %>%
  select(V1, V2, V3, V9)%>%
  mutate(sp = "POC")
colnames(PMEA_TE) <- c("chr","start","end","bp_overlap","sp")
PEVE_TE <- read.table("../E-Peve/output/18-Peve-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_PEVE/PEVE_piRNA_cluster_TEs_intersect.txt", header = F) %>%
  select(V1, V2, V3, V9)%>%
  mutate(sp = "POR")
colnames(PEVE_TE) <- c("chr","start","end","bp_overlap","sp")

merged_clust_over <- rbind(APUL_gene,PMEA_gene,PEVE_gene) %>%
  select(bp_overlap, sp) %>%
  group_by(sp) %>%
  summarise(sum(bp_overlap)) %>%
  mutate(element = "gene") %>% rbind(rbind(APUL_TE,PMEA_TE,PEVE_TE) %>%
  select(bp_overlap, sp) %>%
  group_by(sp) %>%
  summarise(sum(bp_overlap)) %>%
  mutate(element = "TE")) %>% rbind(rbind(APUL_m_clust,PMEA_m_clust,PEVE_m_clust) %>%
  select(bp_overlap, sp) %>%
  group_by(sp) %>%
  summarise(sum(bp_overlap)) %>%
  mutate(element = "clust")) %>%
  spread(element, "sum(bp_overlap)") %>%
  mutate(intergenic = clust-(gene + TE)) %>%
  select(sp, gene, TE, intergenic) %>%
  gather(key = "element", value = "bp_count", 2:4) %>%
  group_by(sp) %>% 
  mutate(perc = (bp_count/sum(bp_count)*100)) %>%
  ungroup()

cp <- coord_polar(theta = "y", start = 0)
cp$is_free <- function() TRUE

clust_plot <- ggplot(merged_clust_over, aes(x="", y=perc, alpha=element, fill = sp)) +
  geom_bar(stat="identity", width=1) +
  cp +
  facet_wrap(~sp, scales = "free", ncol = 1) +
  geom_text(aes(label = ifelse(perc > 15 , paste0(round(perc), "%"), "")), 
            position = position_stack(vjust = 0.5),
            show.legend = FALSE) +
  geom_text(aes(x = 1.7, label = ifelse(perc < 15 , paste0(round(perc), "%"), "")), 
            position = position_stack(vjust = 0.5),
            show.legend = FALSE) +
  scale_fill_manual(values = c("#408EC6", "#7A2048", "#1E2761"), labels = c("*A. pulchra*", "*P. tuahiniensis*", "*P. evermanni*")) +
  scale_alpha_discrete(range = c(.4,1), labels = c("Gene-derived piRNAs", "Intergenic piRNAs", "TE-derived piRNAs")) +
  theme_void() +
  theme(aspect.ratio = 1,
        panel.spacing = unit(0, "lines"),
        strip.text = element_blank(), 
        legend.title = element_blank(),
        legend.text = element_markdown(),
        legend.position = "right",
        legend.box = "vertical");clust_plot

ggsave(filename = "figures/cluster_overlap_plot.png", clust_plot, width = 7, height = 5)
```

## GO enrichment in clusters
```{r}

library(topGO)
library(clusterProfiler)
library(pathview)

# Load annotation
GO_APUL <- read.delim("../data/Apul-rna-GO.tsv", sep = "\t", header = T)  %>%
  select(query, GeneOntologyIDs) ## Need to replace transcript names with gene names
GO_PMEA <- read.delim("../data/Pmea-rna-GO.tsv", sep = "\t", header = T)  %>%
  mutate(gene_name = query) %>%
  select(gene_name, GeneOntologyIDs)
write_tsv(GO_PMEA, file = "../data/Pmea_gene2go.tab")
GO_PEVE <- read.delim("../data/Peve-rna-GO.tsv", sep = "\t", header = T)  %>%
  select(query, GeneOntologyIDs) %>%
  mutate(gene_name = gsub("Parent=","", query)) %>%
  select(gene_name, GeneOntologyIDs)
write_tsv(GO_PEVE, file = "../data/Peve_gene2go.tab")

## Relate gene_ID and transcripts

#grep $'\tmRNA\t' ../genome/APUL/GCF_013753865.1_Amil_v2.1_genomic.gff |
#awk '
#BEGIN { FS = "\t"; OFS = "\t" } 
#{
#  if ($9 ~ /gene=/ && $9 ~ /ID=/) {
#    match($9, /gene=[^;]*/);
#    gene = substr($9, RSTART+5, RLENGTH-5);
#    match($9, /ID=[^;]*/);
#    transcript_id = substr($9, RSTART+3, RLENGTH-3);
#    print gene, transcript_id;
#  }
#}' > ../data/APUL_gene2transcript.tab

gene2trans <- read.table("../data/APUL_gene2transcript.tab", header = F) %>%
  mutate(query = gsub("rna-", "", V2)) %>%
  select(query, V1)
colnames(gene2trans) <- c("query", "gene_name")

GO_APUL <- GO_APUL %>%
  left_join(.,gene2trans) %>%
  select(gene_name, GeneOntologyIDs)
write_tsv(GO_APUL, file = "../data/Apul_gene2go.tab")

# Get overlapping genes
APUL_gene_IDs <- read.table("../D-Apul/output/18-Apul-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_APUL/APUL.merged.clusters.gene.overlap.tab", header = F) %>%
  select(V12) %>%
  separate(V12, into = paste0("col", 1:6), sep = ";", fill = "right") %>%
  select(starts_with("col")) %>%
  pivot_longer(cols = everything(), names_to = "key", values_to = "value") %>%
  filter(grepl("gene=", value)) %>%
  mutate(gene_name = sub("gene=", "", value)) %>%
  select(gene_name) 
apul_clust_genes <- as.character(APUL_gene_IDs$gene_name)
apul_all_genes <- as.character(GO_APUL$gene_name)
Apul_GeneList <- factor(as.integer(apul_all_genes %in% apul_clust_genes))
names(Apul_GeneList) <- apul_all_genes


PMEA_gene_IDs <- read.table("../F-Pmea/output/18-Pmea-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_PMEA/PMEA.merged.clusters.gene.overlap.tab", header = F) %>%
  select(V12) %>%
  mutate(gene_name = sub("ID=", "", V12)) %>%
  select(gene_name)

pmea_clust_genes <- as.character(PMEA_gene_IDs$gene_name)
pmea_all_genes <- as.character(GO_PMEA$gene_name)
Pmea_GeneList <- factor(as.integer(pmea_all_genes %in% pmea_clust_genes))
names(Pmea_GeneList) <- pmea_all_genes

PEVE_gene_IDs <- read.table("../E-Peve/output/18-Peve-piRNA/0_piRNA_pipeline_proTRAC/proTRAC_results_PEVE/PEVE.merged.clusters.gene.overlap.tab", header = F) %>%
  select(V12) %>%
  separate(V12, into = paste0("col", 1:6), sep = ";", fill = "right") %>%
  select(starts_with("col")) %>%
  pivot_longer(cols = everything(), names_to = "key", values_to = "value") %>%
  filter(grepl("ID=", value)) %>%
  mutate(gene_name = sub("ID=", "", value)) %>%
  select(gene_name)

peve_clust_genes <- as.character(PEVE_gene_IDs$gene_name)
peve_all_genes <- as.character(GO_PEVE$gene_name)
Peve_GeneList <- factor(as.integer(peve_all_genes %in% peve_clust_genes))
names(Peve_GeneList) <- peve_all_genes
```
 
# Run TopGO
```{r}
topDiffGenes <- function(allScore) {
return(allScore < 0.05)}

#GOnode<-10

APUL_gene2go<-readMappings("../data/Apul_gene2go.tab", IDsep=";", sep="\t")
PMEA_gene2go<-readMappings("../data/Pmea_gene2go.tab", IDsep=";", sep="\t")
PEVE_gene2go<-readMappings("../data/Peve_gene2go.tab", IDsep=";", sep="\t")

# Biological Process
GO_APUL_BP <-new("topGOdata", ontology="BP", gene2GO=APUL_gene2go, allGenes=Apul_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

GO_PMEA_BP <-new("topGOdata", ontology="BP", gene2GO=PMEA_gene2go, allGenes=Pmea_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

GO_PEVE_BP <-new("topGOdata", ontology="BP", gene2GO=PEVE_gene2go, allGenes=Peve_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

# Run tests 
GO_APUL_BP_FE <- runTest(GO_APUL_BP, algorithm="weight01", statistic="fisher")
GO_PMEA_BP_FE <- runTest(GO_PMEA_BP, algorithm="weight01", statistic="fisher")
GO_PEVE_BP_FE <- runTest(GO_PEVE_BP, algorithm="weight01", statistic="fisher")

##Table of Results
GO_APUL_BP_En <- GenTable(GO_APUL_BP, Fisher = GO_APUL_BP_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_APUL_BP_En$Fisher<-as.numeric(GO_APUL_BP_En$Fisher)

GO_PEVE_BP_En <- GenTable(GO_PEVE_BP, Fisher = GO_PEVE_BP_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_PEVE_BP_En$Fisher<-as.numeric(GO_PEVE_BP_En$Fisher)

GO_PMEA_BP_En <- GenTable(GO_PMEA_BP, Fisher = GO_PMEA_BP_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_PMEA_BP_En$Fisher<-as.numeric(GO_PMEA_BP_En$Fisher)

##Filter for Significant Terms
GO_APUL_BP_En<-GO_APUL_BP_En[GO_APUL_BP_En$Fisher<0.01,]
GO_PEVE_BP_En<-GO_PEVE_BP_En[GO_PEVE_BP_En$Fisher<0.01,]
GO_PMEA_BP_En<-GO_PMEA_BP_En[GO_PMEA_BP_En$Fisher<0.01,]

##Fix order of terms for plotting
GO_APUL_BP_En$Term<-factor(GO_APUL_BP_En$Term, levels = rev(GO_APUL_BP_En$Term))
GO_PEVE_BP_En$Term<-factor(GO_PEVE_BP_En$Term, levels = rev(GO_PEVE_BP_En$Term))
GO_PMEA_BP_En$Term<-factor(GO_PMEA_BP_En$Term, levels = rev(GO_PMEA_BP_En$Term))

##Add Contrast, Timepoint, and Direction columns
GO_APUL_BP_En$Species<-"*A. pulchra*"
GO_APUL_BP_En$ontology<-"Biological Process" 

GO_PMEA_BP_En$Species<-"*P. tuahiniensis*"
GO_PMEA_BP_En$ontology<-"Biological Process" 

GO_PEVE_BP_En$Species<-"*P. evermanni*"
GO_PEVE_BP_En$ontology<-"Biological Process"  

Reduce(intersect, list(GO_APUL_BP_En$GO.ID, GO_PEVE_BP_En$GO.ID))
Reduce(intersect, list(GO_APUL_BP_En$GO.ID, GO_PMEA_BP_En$GO.ID))
Reduce(intersect, list(GO_PMEA_BP_En$GO.ID, GO_PEVE_BP_En$GO.ID))

# Cellular components
GO_APUL_CC <-new("topGOdata", ontology="CC", gene2GO=APUL_gene2go, allGenes=Apul_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

GO_PMEA_CC <-new("topGOdata", ontology="CC", gene2GO=PMEA_gene2go, allGenes=Pmea_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

GO_PEVE_CC <-new("topGOdata", ontology="CC", gene2GO=PEVE_gene2go, allGenes=Peve_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

# Run tests 
GO_APUL_CC_FE <- runTest(GO_APUL_CC, algorithm="weight01", statistic="fisher")
GO_PMEA_CC_FE <- runTest(GO_PMEA_CC, algorithm="weight01", statistic="fisher")
GO_PEVE_CC_FE <- runTest(GO_PEVE_CC, algorithm="weight01", statistic="fisher")

##Table of Results
GO_APUL_CC_En <- GenTable(GO_APUL_CC, Fisher = GO_APUL_CC_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_APUL_CC_En$Fisher<-as.numeric(GO_APUL_CC_En$Fisher)

GO_PEVE_CC_En <- GenTable(GO_PEVE_CC, Fisher = GO_PEVE_CC_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_PEVE_CC_En$Fisher<-as.numeric(GO_PEVE_CC_En$Fisher)

GO_PMEA_CC_En <- GenTable(GO_PMEA_CC, Fisher = GO_PMEA_CC_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_PMEA_CC_En$Fisher<-as.numeric(GO_PMEA_CC_En$Fisher)

##Filter for Significant Terms
GO_APUL_CC_En<-GO_APUL_CC_En[GO_APUL_CC_En$Fisher<0.01,]
GO_PEVE_CC_En<-GO_PEVE_CC_En[GO_PEVE_CC_En$Fisher<0.01,]
GO_PMEA_CC_En<-GO_PMEA_CC_En[GO_PMEA_CC_En$Fisher<0.01,]

##Fix order of terms for plotting
GO_APUL_CC_En$Term<-factor(GO_APUL_CC_En$Term, levels = rev(GO_APUL_CC_En$Term))
GO_PEVE_CC_En$Term<-factor(GO_PEVE_CC_En$Term, levels = rev(GO_PEVE_CC_En$Term))
GO_PMEA_CC_En$Term<-factor(GO_PMEA_CC_En$Term, levels = rev(GO_PMEA_CC_En$Term))

##Add Contrast, Timepoint, and Direction columns
GO_APUL_CC_En$Species<-"*A. pulchra*"
GO_APUL_CC_En$ontology<-"Cellular Component" 

GO_PMEA_CC_En$Species<-"*P. tuahiniensis*"
GO_PMEA_CC_En$ontology<-"Cellular Component" 

GO_PEVE_CC_En$Species<-"*P. evermanni*"
GO_PEVE_CC_En$ontology<-"Cellular Component" 

Reduce(intersect, list(GO_APUL_CC_En$GO.ID, GO_PEVE_CC_En$GO.ID))
Reduce(intersect, list(GO_APUL_CC_En$GO.ID, GO_PMEA_CC_En$GO.ID))
Reduce(intersect, list(GO_PMEA_CC_En$GO.ID, GO_PEVE_CC_En$GO.ID))

# Molecular functions
## Create datasets
GO_APUL_MF <-new("topGOdata", ontology="MF", gene2GO=APUL_gene2go, allGenes=Apul_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

GO_PMEA_MF <-new("topGOdata", ontology="MF", gene2GO=PMEA_gene2go, allGenes=Pmea_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

GO_PEVE_MF <-new("topGOdata", ontology="MF", gene2GO=PEVE_gene2go, allGenes=Peve_GeneList, annot = annFUN.gene2GO, geneSel=topDiffGenes)

# Run tests 
GO_APUL_MF_FE <- runTest(GO_APUL_MF, algorithm="weight01", statistic="fisher")
GO_PMEA_MF_FE <- runTest(GO_PMEA_MF, algorithm="weight01", statistic="fisher")
GO_PEVE_MF_FE <- runTest(GO_PEVE_MF, algorithm="weight01", statistic="fisher")

##Table of Results
GO_APUL_MF_En <- GenTable(GO_APUL_MF, Fisher = GO_APUL_MF_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_APUL_MF_En$Fisher<-as.numeric(GO_APUL_MF_En$Fisher)

GO_PEVE_MF_En <- GenTable(GO_PEVE_MF, Fisher = GO_PEVE_MF_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_PEVE_MF_En$Fisher<-as.numeric(GO_PEVE_MF_En$Fisher)

GO_PMEA_MF_En <- GenTable(GO_PMEA_MF, Fisher = GO_PMEA_MF_FE,
  orderBy = "Fisher",  topNodes = 100, numChar = 51)
GO_PMEA_MF_En$Fisher<-as.numeric(GO_PMEA_MF_En$Fisher)

##Filter for Significant Terms
GO_APUL_MF_En<-GO_APUL_MF_En[GO_APUL_MF_En$Fisher<0.01,]
GO_PEVE_MF_En<-GO_PEVE_MF_En[GO_PEVE_MF_En$Fisher<0.01,]
GO_PMEA_MF_En<-GO_PMEA_MF_En[GO_PMEA_MF_En$Fisher<0.01,]

##Fix order of terms for plotting
GO_APUL_MF_En$Term<-factor(GO_APUL_MF_En$Term, levels = rev(GO_APUL_MF_En$Term))
GO_PEVE_MF_En$Term<-factor(GO_PEVE_MF_En$Term, levels = rev(GO_PEVE_MF_En$Term))
GO_PMEA_MF_En$Term<-factor(GO_PMEA_MF_En$Term, levels = rev(GO_PMEA_MF_En$Term))

##Add Contrast, Timepoint, and Direction columns
GO_APUL_MF_En$Species<-"*A. pulchra*"
GO_APUL_MF_En$ontology<-"Molecular Function" 

GO_PMEA_MF_En$Species<-"*P. tuahiniensis*"
GO_PMEA_MF_En$ontology<-"Molecular Function" 

GO_PEVE_MF_En$Species<-"*P. evermanni*"
GO_PEVE_MF_En$ontology<-"Molecular Function" 

Reduce(intersect, list(GO_APUL_MF_En$GO.ID, GO_PEVE_MF_En$GO.ID))
Reduce(intersect, list(GO_APUL_MF_En$GO.ID, GO_PMEA_MF_En$GO.ID))
Reduce(intersect, list(GO_PMEA_MF_En$GO.ID, GO_PEVE_MF_En$GO.ID))

all_GOs <- rbind(GO_APUL_BP_En, GO_PMEA_BP_En, GO_PEVE_BP_En,GO_APUL_MF_En, GO_PMEA_MF_En, GO_PEVE_MF_En,GO_APUL_CC_En, GO_PMEA_CC_En, GO_PEVE_CC_En) %>%
  mutate(prop.sig.genes = Significant/Annotated)
all_GOs$ontology <- as.factor(all_GOs$ontology)
all_GOs$Species <- as.factor(all_GOs$Species)

str(all_GOs)

```


```{r}
####Plot Feature Sizes
axis.title.sz=18
axis.txt.sz=14
leg.title.sz=15
leg.txt.sz=12
levels.sz=7
sig.sz=5
panel.lab.sz=25
point.sz=4
bar.sz=1
cap.sz=0.3



# Biological Processes 

GO_plot<-ggplot(all_GOs, aes(x = Term, y = -log10(Fisher), size = prop.sig.genes, fill = -log10(Fisher))) +
  #expand_limits(y = 1.5) +
  #ylim(1, 7.25) +
  # Add horizontal lines with a single aesthetic value
  geom_hline(yintercept = -log10(0.01), linetype = "longdash", colour = "black", linewidth = .6) +
  geom_hline(yintercept = -log10(0.001), linetype = "solid", colour = "black", linewidth = .6) +
  geom_point(shape = 21) + 
  scale_size(range = c(2, 12)) + 
  scale_fill_continuous(low = "#1AD3D1FF", high = "#4686FBFF") +
  xlab('') + 
  ylab('Enrichment score') +
  labs(caption = 'Cut-off lines at p=0.01 and p=0.001') +
  theme_bw() +
  facet_grid(vars(ontology), vars(Species), scales = "free", space = "free_y") + 
  theme(
    legend.position = 'right',
    legend.background = element_rect(),
    plot.title = element_text(size = axis.title.sz, face = 'bold', vjust = 1, colour = "black"),
    plot.subtitle = element_text(size = axis.txt.sz, vjust = 1, colour = "black"),
    plot.caption = element_text(size = leg.txt.sz, vjust = 1, face = 'italic', colour = "black"),
    axis.text.x = element_text(size = axis.txt.sz, hjust = 1.10, colour = "black"),
    axis.text.y = element_text(size = axis.txt.sz, vjust = 0.5, colour = "black"),
    axis.title.x = element_text(size = axis.title.sz, colour = "black"),
    axis.title.y = element_text(size = axis.title.sz, colour = "black"),
    axis.line = element_line(colour = 'black'),
    legend.key = element_blank(),
    legend.text = element_text(size = leg.txt.sz),
    title = element_text(size = leg.txt.sz),
    strip.background = element_blank(), 
    strip.text = element_markdown(),
    panel.border = element_rect(colour = "black", fill = "NA"), 
    panel.background = element_rect(fill = NA, color = "black") 
  ) +
  coord_flip();GO_plot

ggsave(filename = "figures/cluster_GOs_plot.png", GO_plot, width = 15, height = 20)
```
# Plot overlaps
```{r}
library(ggvenn)
library(cowplot)

GOs_BP_wide <- list(APUL=GO_APUL_BP_En$GO.ID,PTUA=GO_PMEA_BP_En$GO.ID,PEVE=GO_PEVE_BP_En$GO.ID)
GOs_CC_wide <- list(APUL=GO_APUL_CC_En$GO.ID,PTUA=GO_PMEA_CC_En$GO.ID,PEVE=GO_PEVE_CC_En$GO.ID)
GOs_MF_wide <- list(APUL=GO_APUL_MF_En$GO.ID,PTUA=GO_PMEA_MF_En$GO.ID,PEVE=GO_PEVE_MF_En$GO.ID)

gv1 <- ggvenn(GOs_BP_wide, fill_color = c("#408EC6", "#7A2048", "#1E2761"),
  stroke_size = 0.5, set_name_size = 0, show_percentage = F, stroke_color = F)
gv2 <- ggvenn(GOs_CC_wide, fill_color = c("#408EC6", "#7A2048", "#1E2761"),
  stroke_size = 0.5, set_name_size = 0, show_percentage = F, stroke_color = F 
  )
gv3 <- ggvenn(GOs_MF_wide, fill_color = c("#408EC6", "#7A2048", "#1E2761"),
  stroke_size = 0.5, set_name_size = 0, show_percentage = F, stroke_color = F 
  )

GO_venn <- plot_grid(gv1,gv2, gv3, nrow = 1, align = "h", labels = c("Biological Process", "Cellular Component", "Molecular Function"), vjust = 7, hjust = -0.2)

ggsave("figures/GO_overlaps_vennDiagram.png", GO_venn, height = 4.6, width = 8)
```

