---
title: "08.2.1-Peve-sRNAseq-ShortStack-FastA-extraction"
author: "Sam White"
date: "2024-05-22"
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
bibliography: references.bib
link-citations: true
---

```{r setup, include=FALSE}
library(knitr)
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
)
```

Extract mature miRNAs identified with matches to miRBase by ShortStack in [08.2-Peve-sRNAseq-ShortStack-31bp-fastp-merged.Rmd](./08.2-Peve-sRNAseq-ShortStack-31bp-fastp-merged.Rmd) to FastA.


------------------------------------------------------------------------


# Create a Bash variables file

This allows usage of Bash variables across R Markdown chunks.

```{r save-bash-variables-to-rvars-file, engine='bash', eval=TRUE}
{
echo "#### Assign Variables ####"
echo ""

echo "# Trimmed FastQ naming pattern"

echo "# Data directories"
echo 'export deep_dive_dir=/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive'
echo 'export shortstack_dir="${deep_dive_dir}/E-Peve/output/08.2-Peve-sRNAseq-ShortStack-31bp-fastp-merged/ShortStack_out"'
echo 'export output_dir_top=${deep_dive_dir}/E-Peve/output/08.2.1-Peve-sRNAseq-ShortStack-FastA-extraction'
echo ""

echo "# Input/Output files"
echo 'export output_fasta="mature-miRBase-matches.fasta"'
echo 'export shortstack_fasta="mir.fasta"'
echo 'export shortstack_fasta_index="mir.fasta.fai"'
echo 'export shortstack_fixed_fasta="mir-coords-fixed.fasta"'
echo 'export shortstack_fixed_fasta_index="mir-coords-fixed.fasta.fai"'

echo 'export shortstack_results_file="Results.txt"'
echo 'export regions="mature-miRBase-regions.txt"'

echo ""

echo "# Set number of CPUs to use"
echo 'export threads=40'
echo ""

echo "# Programs"
echo 'export samtools=/home/shared/samtools-1.12/samtools'


} > .bashvars

cat .bashvars
```

# Examine `Results.txt`

## Head
```{r head-results.txt, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

head "${shortstack_dir}/${shortstack_results_file}" | column -t
```

## Columns of interest

Column 1: Region of miRNA match

Column 20: ShortStack miRNA? Y/N

Column 21: Match to miRBase? NA or miRBase match

```{r cols-of-interest, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

awk '{print $1"\t"$20"\t"$21}' "${shortstack_dir}/${shortstack_results_file}" | head | column -t
```

## miRNAs of interest
```{r miRNAs-of-interest, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

awk '$20 == "Y" && $21 != "NA" {print $1"\t"$20"\t"$21}' "${shortstack_dir}/${shortstack_results_file}" | head | column -t

echo ""
echo "------------------------------------------"
echo ""
echo "Number of miRNAs matching miRBase:"

awk '$20 == "Y" && $21 != "NA" {print $1"\t"$20"\t"$21}' "${shortstack_dir}/${shortstack_results_file}" | wc -l 
```

# Examine ShortStack miRNA FastA

## Head FastA

```{r head-FastA, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

grep "^>" "${shortstack_dir}/${shortstack_fasta}" | head
```

# Fix FastA description starting coordinates

Needed, due to [bug in code](https://github.com/MikeAxtell/ShortStack/issues/153#issuecomment-2122897486) (GitHub Issue) which incorrectly calculates the starting coordinates in the FastA output. All other files where start/stop coordinates are conveyed are correct.

The incorrect starting coordinates cause an issue in downstream manipulation, because the FastA headers need to match the ShortStack results file.

```{r fix-FastA-coordinates, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

awk '
/^>/ {
    # Split the line into main parts based on "::" delimiter
    split($0, main_parts, "::")
    
    # Extract the coordinate part and strand information separately
    coordinates_strand = main_parts[2]
    split(coordinates_strand, coord_parts, "[:-]")
    
    # Determine if the strand information is present and extract it
    strand = ""
    if (substr(coordinates_strand, length(coordinates_strand)) ~ /[\(\)\-\+]/) {
        strand = substr(coordinates_strand, length(coordinates_strand) - 1)
        coordinates_strand = substr(coordinates_strand, 1, length(coordinates_strand) - 2)
        split(coordinates_strand, coord_parts, "[:-]")
    }
    
    # Increment the starting coordinate by 1
    new_start = coord_parts[2] + 1
    
    # Reconstruct the description line with the new starting coordinate
    new_description = main_parts[1] "::" coord_parts[1] ":" new_start "-" coord_parts[3] strand
    
    # Print the modified description line
    print new_description
    
    # Skip to the next line to process the sequence line
    next
}

# For sequence lines, print them as-is
{
    print
}
' "${shortstack_dir}/${shortstack_fasta}" \
> "${shortstack_dir}/${shortstack_fixed_fasta}"

diff "${shortstack_dir}/${shortstack_fasta}" \
"${shortstack_dir}/${shortstack_fixed_fasta}" \
| head


```

# Create regions file for use with `samtools`

## Make FastA index

```{r FastA-index, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

${samtools} faidx "${shortstack_dir}/${shortstack_fixed_fasta}"


head "${shortstack_dir}/${shortstack_fixed_fasta_index}"
```

## Construct regions of miRBase matches for FastA index


```{r contrsuct-regions-file, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

# Make output directory, if it doesn't exist
mkdir --parents "${output_dir_top}"

{
   awk '$20 == "Y" && $21 != "NA" {print $2}' "${shortstack_dir}/${shortstack_results_file}" \
   | grep --fixed-strings --file - "${shortstack_dir}/${shortstack_fixed_fasta_index}" \
   | awk '{print $1}'
} \
> "${output_dir_top}/${regions}"

head "${output_dir_top}/${regions}"

```



# Extract FastAs

```{r extract-FastAs, engine='bash', eval=TRUE}
# Load bash variables into memory
source .bashvars

${samtools} faidx "${shortstack_dir}/${shortstack_fixed_fasta}" \
--region-file "${output_dir_top}/${regions}" \
> "${output_dir_top}/${output_fasta}"

head "${output_dir_top}/${output_fasta}"

echo ""
echo ""
echo ""
echo "Number of FastA sequences:"
grep "^>" --count "${output_dir_top}/${output_fasta}"

```