---
title: "17-Apul-ShortStack-miRdeep2-comparison"
author: "Kathleen Durkin"
date: "2024-05-17"
output: 
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
  github_document:
    toc: true
    number_sections: true
  html_document:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
---

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

## intersectBed

Examine our input files (intersectBed accepts .bed and .gff files)

```{r view-OG-input-files, engine='bash', eval=TRUE}
head -5 ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/mirna_results_03_04_2024_t_13_00_39/novel_mature_03_04_2024_t_13_00_39_score-50_to_na.bed
echo ""
head -5 ../output/13.2.1-Apul-sRNAseq-ShortStack-31bp-fastp-merged-cnidarian_miRBase/ShortStack_out/Results.gff3
```

We need to get two input files that contain only mature miRNAs and have the same sequence naming conventions for matching purposes. That means the miRdeep2 mature miRNAs file need to be reformatted, and the ShortStack full results file needs to be filtered.

```{r prep-input-files, engine='bash', eval=TRUE}
# remove header lines, shorten chromosome ref names to match those used in ShortStack output, and ensure final file is tab-delimited instead of space-delimited
tail -n +5 ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/mirna_results_03_04_2024_t_13_00_39/novel_mature_03_04_2024_t_13_00_39_score-50_to_na.bed | \
awk 'BEGIN { FS = "\t" } { sub(/\.1.*/, ".1", $1); print }' | \
tr ' ' '\t' > ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/mirna_results_03_04_2024_t_13_00_39/novel_mature_03_04_2024_t_13_00_39_score-50_to_na_formatted.bed

# filter full results to obtain a gff file of only the mature miRNAs
awk -F'\t' '$3 == "mature_miRNA"' ../output/13.2.1-Apul-sRNAseq-ShortStack-31bp-fastp-merged-cnidarian_miRBase/ShortStack_out/Results.gff3 > ../output/13.2.1-Apul-sRNAseq-ShortStack-31bp-fastp-merged-cnidarian_miRBase/ShortStack_out/Results_mature.gff3
```

Check the files

```{r view-prepped-input-files, engine='bash', eval=TRUE}
head -5 ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/mirna_results_03_04_2024_t_13_00_39/novel_mature_03_04_2024_t_13_00_39_score-50_to_na_formatted.bed
echo ""
head -5 ../output/13.2.1-Apul-sRNAseq-ShortStack-31bp-fastp-merged-cnidarian_miRBase/ShortStack_out/Results_mature.gff3
```

Looks good! Now we can input into intersectBed.

```{r intersectBed, engine='bash', eval=TRUE}
# intersectBed to ID sequences in the miRdeep2 mature miRNA output that match mature miRNAs ID'd by ShortStack
# -wa and -wb ensure we recieve full annotations from both input files in the output
/home/shared/bedtools2/bin/intersectBed \
-wa \
-wb \
-a ../output/13.2.1-Apul-sRNAseq-ShortStack-31bp-fastp-merged-cnidarian_miRBase/ShortStack_out/Results_mature.gff3 \
-b ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/mirna_results_03_04_2024_t_13_00_39/novel_mature_03_04_2024_t_13_00_39_score-50_to_na_formatted.bed \
&> ../output/17-Apul-ShortStack-miRdeep2-comparison/ShortStack_miRdeep2_mature_intersect.bed

```

## Results

```{r count-intersects, engine='bash', eval=TRUE}
echo "Number of ShortStack mature miRNAs:"
wc -l < ../output/13.2.1-Apul-sRNAseq-ShortStack-31bp-fastp-merged-cnidarian_miRBase/ShortStack_out/Results_mature.gff3
echo ""
echo "Number of ShortStack mature miRNAs also identified by miRdeep2:"
wc -l < ../output/17-Apul-ShortStack-miRdeep2-comparison/ShortStack_miRdeep2_mature_intersect.bed
```

While looking through the output file I noticed that two of the intersects originate from the same cluster... not really sure what that's about...

```{r, engine='bash', eval=TRUE}
head -21 ../output/17-Apul-ShortStack-miRdeep2-comparison/ShortStack_miRdeep2_mature_intersect.bed | tail -2
```

It looks like they have pretty much identical entries (cluster, coordinates, etc.), except the end of the miRdeep locus name (NC_058072.1_Acropora_millepora_isolate_JS-1_chromosome_7\_Amil_v2.1_whole_genome_shotgun_sequence\_***295614*** vs NC_058072.1_Acropora_millepora_isolate_JS-1_chromosome_7\_Amil_v2.1_whole_genome_shotgun_sequence\_***295613***) and the miRdeep2 "score" value assigned to them (this is the column following the miRdeep2 locus name, 12623.6 vs 1.1). We can also check these two loci in the full miRdeep2 output.
```{r check-mirdeep-output-dupes, engine='bash', eval=TRUE}

# View full mirdeep2 output for these two loci
awk -F'\t' '$1 == "NC_058072.1_Acropora_millepora_isolate_JS-1_chromosome_7_Amil_v2.1_whole_genome_shotgun_sequence_295614"' ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/parsable-result_03_04_2024_t_13_00_39.csv
echo""
awk -F'\t' '$1 == "NC_058072.1_Acropora_millepora_isolate_JS-1_chromosome_7_Amil_v2.1_whole_genome_shotgun_sequence_295613"' ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/parsable-result_03_04_2024_t_13_00_39.csv

echo""
echo""

echo "mature miRNA sequences for these two loci:"
awk -F'\t' '$1 == "NC_058072.1_Acropora_millepora_isolate_JS-1_chromosome_7_Amil_v2.1_whole_genome_shotgun_sequence_295614"' ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/parsable-result_03_04_2024_t_13_00_39.csv | awk '{print $13}'

awk -F'\t' '$1 == "NC_058072.1_Acropora_millepora_isolate_JS-1_chromosome_7_Amil_v2.1_whole_genome_shotgun_sequence_295613"' ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/parsable-result_03_04_2024_t_13_00_39.csv | awk '{print $13}'

echo ""
echo ""

echo "precursor miRNA sequences for these two loci:"
awk -F'\t' '$1 == "NC_058072.1_Acropora_millepora_isolate_JS-1_chromosome_7_Amil_v2.1_whole_genome_shotgun_sequence_295614"' ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/parsable-result_03_04_2024_t_13_00_39.csv | awk '{print $15}'

awk -F'\t' '$1 == "NC_058072.1_Acropora_millepora_isolate_JS-1_chromosome_7_Amil_v2.1_whole_genome_shotgun_sequence_295613"' ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/parsable-result_03_04_2024_t_13_00_39.csv | awk '{print $15}'
```
Interesting... The two loci have identical mature miRNAs, but different precursors!


Let's set that aside for now and do some quick investigation of the miRdeep2 evaluation criteria for all of these ShortStack/miRdeep2 shared miRNAs. This could give us an idea of what thresholds may be appropriate for filtering the miRdeep2 output. 
```{r check-mirdeep2-output-all-shared, engine='bash', eval=TRUE}
mirdeepIDs=$(awk '{print $13}' ../output/17-Apul-ShortStack-miRdeep2-comparison/ShortStack_miRdeep2_mature_intersect.bed)

head -1 ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/parsable-result_03_04_2024_t_13_00_39.csv > ../output/17-Apul-ShortStack-miRdeep2-comparison/intersect_miRdeep2_stats.txt
while IFS= read -r id; do
  # Use awk to process fileA and match column 1 with the current ID
  awk -F'\t' -v id="$id" '$1 == id {print}' ../output/11.1-Apul-sRNAseq-miRdeep2-31bp-fastp-merged/parsable-result_03_04_2024_t_13_00_39.csv
done <<< "$mirdeepIDs" >> ../output/17-Apul-ShortStack-miRdeep2-comparison/intersect_miRdeep2_stats.txt
```

```{r view-criteria, eval=TRUE}
intersect_miRdeep2_stats <- read.csv("../output/17-Apul-ShortStack-miRdeep2-comparison/intersect_miRdeep2_stats.txt", sep="\t")

intersect_miRdeep2_stats %>% 
  select(miRDeep2.score, significant.randfold.p.value) %>% 
  arrange(desc(miRDeep2.score)) %>%
  kable()
```

The miRdeep2 score is "the log-odds score assigned to the hairpin" and is essentially a probability that the locus is a miRNA (presumably based on the hairpin structure) on, with higher values indicating higher probability. MiRdeep2's default score threshold for miRNA classification is 0, but coral miRNA papers we've seen use myriad thresholds (e.g., 4, 10).

The randfold is also an evaluation of ncrna secondary structure, and we want a significant randfold value (this would indicate high likelihood of miRNA precursory structure)

For some reason there are no miRdeep2 miRNAs with significant randfold p-values for A. pulchra, even though the other species have a bunch of loci with significant p-values. The miRdeep2 score is interesting though! All but 3 of the shared miRNAs have scores >10 -- maybe that could be a good cutoff?
