---
title: "03-hisat2-summer2021_2022"
output: html_document
date: "2024-05-21"
---
Rmd to run `HISAT2` with _Pycnopodia helianthoides_ genome and _Pycnpoodia helianthoides_ coelomocyte RNAseq sequences. 

All trimmed reads from summer 2021 and 2022 live on Raven:       
Summer 2021     
`/home/shared/8TB_HDD_02/graceac9/data/pycno2021`

Summer 2022    
`/home/shared/8TB_HDD_02/graceac9/data/pycno2022` 

Get genome from NCBI into data directory on Raven:        
https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_032158295.1/       

```{r, engine='bash'}
cd ../data

/home/shared/datasets download genome accession GCA_032158295.1 --include gff3,rna,cds,protein,genome,seq-report
```

```{r, engine='bash'}
cd ../data
unzip ncbi_dataset.zip

```

Get gtf file on raven into data folder:    
I downloaded it originally from 

rsync /Users/graciecrandall/Downloads/augustus.hints.gtf graceac9@raven.fish.washington.edu:/home/shared/8TB_HDD_02/graceac9/GitHub/paper-pycno-sswd-2021-2022/data
graceac9@raven.fish.washington.edu's password: 





# 1. BUILD AN INDEX

Follow this code: https://htmlpreview.github.io/?https://github.com/urol-e5/deep-dive/blob/8fd4ad4546d1d95464952f0509406efd9e42ffa0/D-Apul/code/04-Apulcra-hisat.html 

```{bash}
pwd
```

From the gtf, get exon list
```{bash}
/home/shared/hisat2-2.2.1/hisat2_extract_exons.py \
../data/augustus.hints.gtf \
> ../analyses/03-hisat2/m_exon.tab
```

```{bash}
head ../analyses/03-hisat2/m_exon.tab
```

from the gtf, get the splice sites
```{bash}
#!/bin/bash

# This script will extract splice sites from the gtf file

# This is the command to extract splice sites from the gtf file
/home/shared/hisat2-2.2.1/hisat2_extract_splice_sites.py \
../data/augustus.hints.gtf \
> ../analyses/03-hisat2/m_splice_sites.tab

```

use the genome fasta to make an index for alignment
```{bash}
# build an index 
/home/shared/hisat2-2.2.1/hisat2-build \
../data/ncbi_dataset/data/GCA_032158295.1/GCA_032158295.1_ASM3215829v1_genomic.fna \
../data/Phel_genome_index \
--exon ../analyses/03-hisat2/m_exon.tab \
--ss ../analyses/03-hisat2/m_splice_sites.tab \
-p 40 \
../data/augustus.hints.gtf \
2> ../analyses/03-hisat2/hisat2-build_stats.txt
```

```{bash}
tail ../analyses/03-hisat2/hisat2-build_stats.txt
```

# 2. ALIGNMENT
old code that doens't have --dta function and that get the stats output in separate files per library 
```{bash}
#for file in ../data/2021_trimmed/*_R1*fq.gz; do
    # Remove the  part to get the base name
    base=$(basename "$file" _R1_001.fastq.gz.fastp-trim.20220810.fq.gz)
    # Construct the names of the pair of files
    file1=${base}_R1_001.fastq.gz.fastp-trim.20220810.fq.gz
    file2=${base}_R2_001.fastq.gz.fastp-trim.20220810.fq.gz
    # Run the hisat2 command
#    /home/shared/hisat2-2.2.1/hisat2 \
#    -x ../data/Phel_genome_index \
#    --dta \
#    -p 20 \
#    -1 ../data/2021_trimmed/$file1 \
#    -2 ../data/2021_trimmed/$file2 \
#    -S ../analyses/03-hisat2/${base}.sam \
#    2> ../analyses/03-hisat2/hisat.out
#done
```

this code chunk doesn't work because my hisat2 stats are individual files for each library 
```{bash}
#cat ../analyses/03-hisat2/hisat2_stats.txt
```

Convert sam to bam
```{bash}
#for samfile in ../analyses/03-hisat2/*.sam; do
 # bamfile="${samfile%.sam}.bam"
#  sorted_bamfile="${samfile%.sam}.sorted.bam"
 # /home/shared/samtools-1.12/samtools view -bS -@ 20 "$samfile" > "$bamfile"
  #/home/shared/samtools-1.12/samtools sort -@ 20 "$bamfile" -o "$sorted_bamfile"
#  /home/shared/samtools-1.12/samtools index -@ 20 "$sorted_bamfile"
#done
```

2024-05-25   
Should have done the --dta option for alignment based on [Steven's notebook post](https://sr320.github.io/tumbling-oysters/posts/sr320-12/)    

Re-running alignment with that:    
Steven's code: 
```{bash}
#find ../data/fastq/*R2_001.fastp-trim.20230519.fastq.gz \
#| xargs basename -s -S1-TP2_R2_001.fastp-trim.20230519.fastq.gz | xargs -I{} \
#/home/shared/hisat2-2.2.1/hisat2 \
#-x ../output/15-Apul-hisat/GCF_013753865.1_Amil_v2.1.index \
#--dta \
#-p 20 \
#-1 ../data/fastq/{}-S1-TP2_R1_001.fastp-trim.20230519.fastq.gz \
#-2 ../data/fastq/{}-S1-TP2_R2_001.fastp-trim.20230519.fastq.gz \
#-S ../output/15-Apul-hisat/{}.sam \
#2> ../output/15-Apul-hisat/hisat.out
```

Code for mine: 
```{bash}
find ../data/2021_trimmed/*.fq.gz \
| xargs basename -s -.fastq.gz.fastp-trim.20220810.fq.gz | xargs -I{} \
/home/shared/hisat2-2.2.1/hisat2 \
-x ../data/Phel_genome_index \
--dta \
-p 20 \
-1 ../data/2021_trimmed/{}_R1_001.fastq.gz.fastp-trim.20220810.fq.gz \
-2 ../data/2021_trimmed/{}_R2_001.fastq.gz.fastp-trim.20220810.fq.gz \
-S ../analyses/03-hisat2/{}.sam \
2> ../analyses/03-hisat2/hisat.out
```

