---
title: "06.2-cod-RNAseq-alignment-genome"
author: "Kathleen Durkin"
date: "2024-04-16"
always_allow_html: true
output: 
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: 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
)
```

```{r load-packages, eval=TRUE}
library(tidyverse)
library(dplyr)
library(magrittr)
library(knitr)
library(ggplot2)
library(plotly)
```

Code for aligning RNAseq data to reference genome, to be used on [Pacific cod RNAseq data](https://shedurkin.github.io/Roberts-LabNotebook/posts/projects/pacific_cod/2023_12_13_pacific_cod.html). 

- Raw reads found [here](https://owl.fish.washington.edu/nightingales/G_macrocephalus/30-943133806/)
- Trimmed reads: 
- Genome downloaded from [NCBI](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_031168955.1/), stored [here](https://owl.fish.washington.edu/halfshell/genomic-databank/GCF_031168955.1_ASM3116895v1_rna.fna) as a part of lab [genomic resources](https://robertslab.github.io/resources/Genomic-Resources/#gadus-macrocephalus-pacific-cod)
- Genome GTF downloaded from [NCBI](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_031168955.1/)

# Create a Bash variables file

This allows usage of Bash variables (e.g. paths to common directories) across R Markdown chunks.
```{r save-bash-variables-to-rvars-file, engine='bash', eval=TRUE}
{
echo "#### Assign Variables ####"
echo ""

echo "# Data directories"
echo 'export cod_dir=/home/shared/8TB_HDD_02/shedurkin/project-cod-temperature'
echo 'export output_dir_top=${cod_dir}/output/06.2-cod-RNAseq-alignment-genome'
echo 'export genome_fasta_dir=${cod_dir}/data'
echo 'export trimmed_reads_dir=${cod_dir}/output/05-cod-RNAseq-trimming/trimmed-reads'
echo ""


echo "# Input/Output files"
echo 'export genome_fasta_name="GCF_031168955.1_ASM3116895v1_genomic"'
echo 'export genome_fasta="${genome_fasta_dir}/${genome_fasta_name}"'
echo 'export genome_gtf_name="genomic"'
echo 'export genome_gtf="${genome_fasta_dir}/${genome_gtf_name}"'
echo 'export hisat2_exon_name="G_macrocephalus_exon"'
echo 'export hisat2_exon="${hisat2_output_dir}/${hisat_exon_name}"'
echo 'export hisat2_splice_sites_name="G_macrocephalus_splice_sites"'
echo 'export hisat2_splice_sites="${hisat2_output_dir}/${hisat2_splice_sites_name}"'
echo 'export hisat2_index_name="G_macrocephalus_Hisat2_index"'
echo 'export hisat2_index="${hisat2_output_dir}/${hisat2_index}"'


echo "# External data URLs and checksums"
echo 'export genome_fasta_url="https://owl.fish.washington.edu/halfshell/genomic-databank/GCF_031168955.1_ASM3116895v1_genomic.fna"'
echo 'export genome_fasta_checksum="5144890d4eceb0b258d92db3f35c681e"'
echo 'export genome_gtf_url="https://api.ncbi.nlm.nih.gov/datasets/v2alpha/genome/accession/GCF_031168955.1/download?include_annotation_type=GENOME_GTF"'
#echo 'export genome_gtf_checksum="173fb3c159e474391c5c4aa1f7230024"'


echo "# Paths to programs"
echo 'export hisat2_exons=/home/shared/hisat2-2.2.1/hisat2_extract_exons.py'
echo 'export hisat2_splice_sites=/home/shared/hisat2-2.2.1/hisat2_extract_splice_sites.py'
echo 'export hisat2_build=/home/shared/hisat2-2.2.1/hisat2-build'
echo ""


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


echo "# Programs associative array"
echo "declare -A programs_array"
echo "programs_array=("
echo '[hisat2_exons]="${hisat2_exons}" \'
echo '[hisat2_splice_sites]="${hisat2_splice_sites}" \'
echo '[hisat2_build]="${hisat2_build}" \'
echo '[trinity_abund_to_matrix]="${trinity_abund_to_matrix}" \'
echo ")"
} > .bashvars

cat .bashvars
```

# Align to reference genome (Hisat2)

## Retrieving the reference genome and gff

```{r download-genome-fasta, engine='bash'}
# Load bash variables into memory
source .bashvars

wget \
--directory-prefix ${genome_fasta_dir} \
--recursive \
--no-check-certificate \
--continue \
--no-host-directories \
--no-directories \
--no-parent \
--quiet \
--execute robots=off \
--accept "${genome_fasta_name}.fna" ${genome_fasta_url}
```


v NOT CURRENTLY WORKING v, had to download locally and then upload to server for use
```{r download-genome-gtf, engine='bash'}
# Load bash variables into memory
source .bashvars

cd ../data

curl -O "https://api.ncbi.nlm.nih.gov/datasets/v2alpha/genome/accession/GCF_031168955.1/download?include_annotation_type=GENOME_GTF"

#curl -O "${genome_gtf_url}"
```


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

ls -lh "${genome_fasta_dir}"
```


## Verify genome FastA MD5 checksum

```{r verify-genome-fasta-checksum, engine='bash'}
# Load bash variables into memory
source .bashvars

cd "${genome_fasta_dir}"

# Checksums file contains other files, so this just looks for the sRNAseq files.
md5sum --check <<< "${genome_fasta_checksum}  ${genome_fasta_name}.fna"

```

```{r verify-genome-gff-checksum, engine='bash'}
# Load bash variables into memory
source .bashvars

cd "${genome_fasta_dir}"

# Checksums file contains other files, so this just looks for the sRNAseq files.
md5sum --check <<< "${genome_gtf_checksum}  ${genome_gtf_name}.gtf"

```

## Building Index

```{r create-exons-tab-file, engine='bash'}
# Load bash variables into memory
source .bashvars

# Create Hisat2 exons tab file
/home/shared/hisat2-2.2.1/hisat2_extract_exons.py \
../data/genomic.gtf \
> ../output/06.2-cod-RNAseq-alignment-genome/hisat2/G_macrocephalus_exon.tab

# Create Hisat2 exons tab file
#"${programs_array[hisat2_exons]}" \
#"${genome_gtf}.gtf" \
#> "${hisat2_exon}.tab"
```

``` {r create-splice-sites-tab-file, engine='bash'}
# Create Hisat2 splice sites tab file
/home/shared/hisat2-2.2.1/hisat2_extract_splice_sites.py \
../data/genomic.gtf \
> ../output/06.2-cod-RNAseq-alignment-genome/hisat2/G_macrocephalus_splice_sites.tab

#"${programs_array[hisat2_splice_sites]}" \
#"${genome_gtf}.gtf" \
#> "${hisat2_splice_sites}.tab"

```

```{r build-hisat2-index, engine='bash'}
# Build Hisat2 reference index using splice sites and exons
/home/shared/hisat2-2.2.1/hisat2-build \
../data/GCF_031168955.1_ASM3116895v1_genomic.fna \
../output/06.2-cod-RNAseq-alignment-genome/hisat2/G_macrocephalus_index.idx \
--exon ../output/06.2-cod-RNAseq-alignment-genome/hisat2/G_macrocephalus_exon.tab \
--ss ../output/06.2-cod-RNAseq-alignment-genome/hisat2/G_macrocephalus_splice_sites.tab \
-p 20 \
2> ../output/06.2-cod-RNAseq-alignment-genome/hisat2/hisat2-build_stats.txt

#"${programs_array[hisat2_build]}" \
#"${genome_fasta}.fna" \
#"${hisat2_index}.idx" \
#--exon "${hisat2_exon}.tab" \
#--ss "${hisat2_splice_sites}.tab" \
#-p "${threads}" \
#2> "${output_dir_top}/hisat2-build_stats.txt"

```

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

ls -lh "${output_dir_top}"
```


## Alignment
```{r hisat2-alignment, engine='bash'}
# Load bash variables into memory
source .bashvars

## Sample Quantification

# Hisat2 alignments
find ../output/05-cod-RNAseq-trimming/trimmed-reads/*.gz \
| xargs basename -s .flexbar_trim.R_1.fastq.gz | xargs -I{} \
/home/shared/hisat2-2.2.1/hisat2 \
-x ../output/06.2-cod-RNAseq-alignment-genome/hisat2/G_macrocephalus_index.idx \
-p 20 \
-1 ../output/05-cod-RNAseq-trimming/trimmed-reads/{}.flexbar_trim.R_1.fastq.gz \
-2 ../output/05-cod-RNAseq-trimming/trimmed-reads/{}.flexbar_trim.R_2.fastq.gz \
-S ../output/06.2-cod-RNAseq-alignment-genome/hisat2/{}.sam
&> ../output/06.2-cod-RNAseq-alignment-genome/hisat2/{}_hisat2.log



# # Hisat2 alignments
# find ${trimmed_reads_dir}/*.gz \
# | xargs basename -s .flexbar_trim.R_1.fastq.gz | xargs -I{} \
# "${programs_array[hisat2]}" \
# -x "${hisat2_index}.idx" \
# -p 20 \
# -1 ${trimmed_reads_dir}/{}.flexbar_trim.R_1.fastq.gz \
# -2 ${trimmed_reads_dir}/{}.flexbar_trim.R_2.fastq.gz \
# -S ${output_dir_top}/{}.sam
# 
```

```{r sort-sam-files, engine='bash'}
# Sort SAM files, convert to BAM, and index

for samfile in ../output/06.2-cod-RNAseq-alignment-genome/hisat2/*.sam; do
  bamfile="${samfile%.sam}.bam"
  sorted_bamfile="${samfile%.sam}.sorted.bam"
  
  # Check if the output file already exists
  if [[ ! -e "$sorted_bamfile" ]]; then
    # Convert SAM to BAM
    /home/shared/samtools-1.12/samtools view -bS -@ 10 "$samfile" > "$bamfile"
  
    # Sort BAM
    /home/shared/samtools-1.12/samtools sort -@ 10 "$bamfile" -o "$sorted_bamfile"
  
    # Index sorted BAM
    /home/shared/samtools-1.12/samtools index -@ 10 "$sorted_bamfile"
  fi
done
```

```{r count-output-files, engine='bash', eval=TRUE}
# Count the number of samples for which we have sorted bam files -- I should have 79 (one for each input sample)
find ../output/06.2-cod-RNAseq-alignment-genome/hisat2/ -type f -name "*.sorted.bam" | wc -l
```

```{r delete-excess-files, engine='bash'}
# Delete unneccessary index files
rm ../output/06.2-cod-RNAseq-alignment-genome/hisat2/*.ht2

# Delete unneeded SAM files
rm ../output/06.2-cod-RNAseq-alignment-genome/hisat2/*.sam


# # Sort SAM files, convert to BAM, and index
# ${programs_array[samtools_view]} \
# -@ "${threads}" \
# -Su "${sample_name}".sam \
# | ${programs_array[samtools_sort]} - \
# -@ "${threads}" \
# -o "${sample_name}".sorted.bam
# ${programs_array[samtools_index]} "${sample_name}".sorted.bam
# 
# 
# # Delete unneccessary index files
# rm "${genome_index_name}"*.ht2
# 
# # Delete unneeded SAM files
# rm ./*.sam

```

```{r check-output-dir, engine='bash', eval=TRUE}
# View the output files
ls -lh ../output/06.2-cod-RNAseq-alignment-genome/hisat2/
```

# Read Summarization

Will be summarizing reads using [featureCounts](https://www.rdocumentation.org/packages/Rsubread/versions/1.22.2/topics/featureCounts) in the [Rsubread](https://bioconductor.org/packages/release/bioc/vignettes/Rsubread/inst/doc/SubreadUsersGuide.pdf) Bioconductor package

## Exons

```{r featureCounts-exon, engine='bash'}

/home/shared/subread-2.0.5-Linux-x86_64/bin/featureCounts \
-p --countReadPairs \
-T 5 \
-t exon \
-g gene_id \
-a ../data/genomic.gtf \
-o ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-exon/featureCounts_exon_matrix.txt \
../output/06.2-cod-RNAseq-alignment-genome/hisat2/*.sorted.bam \
&> ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-exon/featureCounts_exon.log

```

Save a version of the featureCounts count file with no header (aka remove line 1, which contains the program and command info)
```{r remove-header-feaTureCounts-exon, engine='bash'}
sed '1d' ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-exon/featureCounts_exon_matrix.txt  > ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-exon/featureCounts_exon_matrix_noheader.txt
```

```{r featureCounts-exon-multiqc, engine='bash'}

/home/sam/programs/mambaforge/bin/multiqc \
../output/06.2-cod-RNAseq-alignment-genome/featureCounts-exon/featureCounts_exon_matrix.txt.summary \
-o ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-exon

# View directory contents
ls -lh ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-exon

```

I also want to include the treatment/tank info when plotting alignment rates across samples
```{r plot-alignment-rates-exon, eval=TRUE}
# Load multiqc stats
featureCounts_exon_multiqc <- read.csv("../output/06.2-cod-RNAseq-alignment-genome/featureCounts-exon/multiqc_data/multiqc_featureCounts.txt", sep = '\t')
# Adjust sample name formatting (to prep for join)
featureCounts_exon_multiqc$Sample <- paste("sample_", featureCounts_exon_multiqc$Sample, sep = "")
# Load experimental data
cod_sample_info_OG <- read.csv("../data/DESeq2_Sample_Information.csv")

featureCounts_exon_multiqc_plustreatment <- left_join(cod_sample_info_OG, featureCounts_exon_multiqc, by = c("sample_name" = "Sample")) %>% 
  na.omit()
featureCounts_exon_multiqc_plustreatment <- featureCounts_exon_multiqc_plustreatment[order(featureCounts_exon_multiqc_plustreatment$sample_number),]

ggplot(featureCounts_exon_multiqc_plustreatment,
       aes(x=reorder(sample_name, sample_number), y=percent_assigned, fill=as.factor(temp_treatment))) +
  geom_bar(stat="identity") +
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size=7))

ggplot(featureCounts_exon_multiqc_plustreatment,
       aes(x=reorder(sample_name, sample_number), y=Total, fill=as.factor(temp_treatment))) +
  geom_bar(stat="identity") +
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size=7))

# sample149 is kind of throwing off the visualization, so lets remove and redo
ggplot(featureCounts_exon_multiqc_plustreatment[featureCounts_exon_multiqc_plustreatment$sample_name != "sample_149", ],
       aes(x=reorder(sample_name, sample_number), y=Total, fill=as.factor(temp_treatment))) +
  geom_bar(stat="identity") +
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size=7))

```

## Genes

```{r featureCounts-gene, engine='bash'}

/home/shared/subread-2.0.5-Linux-x86_64/bin/featureCounts \
-p --countReadPairs \
-T 5 \
-t gene \
-g gene_id \
-a ../data/genomic.gtf \
-o ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-gene/featureCounts_gene_matrix.txt \
../output/06.2-cod-RNAseq-alignment-genome/hisat2/*.sorted.bam \
&> ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-gene/featureCounts_gene.log

```

Save a version of the featureCounts count file with no header (aka remove line 1, which contains the program and command info)
```{r remove-header-feaTureCounts-gene, engine='bash'}
sed '1d' ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-gene/featureCounts_gene_matrix.txt  > ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-gene/featureCounts_gene_matrix_noheader.txt
```

```{r featureCounts-gene-multiqc, engine='bash'}

/home/sam/programs/mambaforge/bin/multiqc \
../output/06.2-cod-RNAseq-alignment-genome/featureCounts-gene/featureCounts_gene_matrix.txt.summary \
-o ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-gene

# View directory contents
ls -lh ../output/06.2-cod-RNAseq-alignment-genome/featureCounts-gene

```

I also want to include the treatment/tank info when plotting alignment rates across samples
```{r plot-alignment-rates-gene, eval=TRUE}
# Load multiqc stats
featureCounts_gene_multiqc <- read.csv("../output/06.2-cod-RNAseq-alignment-genome/featureCounts-gene/multiqc_data/multiqc_featureCounts.txt", sep = '\t')
# Adjust sample name formatting (to prep for join)
featureCounts_gene_multiqc$Sample <- paste("sample_", featureCounts_gene_multiqc$Sample, sep = "")

featureCounts_gene_multiqc_plustreatment <- left_join(cod_sample_info_OG, featureCounts_gene_multiqc, by = c("sample_name" = "Sample")) %>% 
  na.omit()
featureCounts_gene_multiqc_plustreatment <- featureCounts_gene_multiqc_plustreatment[order(featureCounts_gene_multiqc_plustreatment$sample_number),]

ggplot(featureCounts_gene_multiqc_plustreatment,
       aes(x=reorder(sample_name, sample_number), y=percent_assigned, fill=as.factor(temp_treatment))) +
  geom_bar(stat="identity") +
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size=7))

ggplot(featureCounts_gene_multiqc_plustreatment,
       aes(x=reorder(sample_name, sample_number), y=Total, fill=as.factor(temp_treatment))) +
  geom_bar(stat="identity") +
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size=7))

# sample149 is kind of throwing off the visualization, so lets remove and redo
ggplot(featureCounts_gene_multiqc_plustreatment[featureCounts_gene_multiqc_plustreatment$sample_name != "sample_149", ],
       aes(x=reorder(sample_name, sample_number), y=Total, fill=as.factor(temp_treatment))) +
  geom_bar(stat="identity") +
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size=7))

```


