---
title: "03-G4-Library1-MinION-Dorado-recall"
author: "Kathleen Durkin"
date: "2025-10-28"
always_allow_html: true
output: 
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: true 
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: 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
)
```


# 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 nanopore_dir=/home/shared/8TB_HDD_02/shedurkin/SIFP-nanopore'
echo 'export output_dir_top=${nanopore_dir}/D-Group4/output/03-G4-Library1-MinION-Dorado-recall'
echo 'export data_dir_top=${nanopore_dir}/D-Group4/data/03-G4-Library1-MinION-Dorado-recall'
echo 'export raw_pod5_dir=${data_dir_top}/pod5_pass'
echo 'export raw_pod5_url1="https://gannet.fish.washington.edu/kdurkin1/SIFP_2025/Group4_MinION/Library1/20250905_1158_MD-101223_FBD08455_b87a1f92/pod5_pass/"'
echo 'export raw_pod5_url2="https://gannet.fish.washington.edu/kdurkin1/SIFP_2025/Group4_MinION/Library1_wash3/20250907_1505_MD-101223_FBD08455_50a8f15e/pod5_pass/"'
echo 'export genome_dir=${nanopore_dir}/data/GCA_965233905.1_jaEunKnig1.1/'

echo ""

echo "# Paths to programs"
echo 'export modkit=/home/shared/dist_modkit_v0.5.1_8fa79e3/modkit'
echo 'export samtools=/home/shared/samtools-1.12/samtools'
echo 'export dorado=/home/shared/dorado-1.2.0-linux-x64/bin/dorado'
echo ""

echo "# Set FastQ filename patterns"
echo "export pod5_pattern='*.pod5'"
echo ""

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

echo "# Input/output files"
echo 'export raw_checksums=checksums.md5'
echo 'export trimmed_checksums=trimmed_fastq_checksums.md5'
echo ""
} > .bashvars

cat .bashvars
```

# Raw reads
## Download raw pod5 files

Reads are downloaded from: 
https://gannet.fish.washington.edu/kdurkin1/SIFP_2025/Group4_MinION/Library1/20250905_1158_MD-101223_FBD08455_b87a1f92/pod5_pass/ and https://gannet.fish.washington.edu/kdurkin1/SIFP_2025/Group4_MinION/Library1_wash3/20250907_1505_MD-101223_FBD08455_50a8f15e/pod5_pass/

Note that there are two directories of data for the Group 4 Library 1 MinION runs. One contains the bulk of the run, the second is everything after wash 3 (the run errored during wash 3)

Note that this directory contains multiple subdirectories, each representing one barcode (specimen) and containing the pod5 files associated with that barcode

The `--cut-dirs 7` command cuts the preceding directory structure (i.e. `nightingales/P_evermanni/30-789513166/`) so that we just end up with the reads.

```{r download-raw-reads, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

# pre-wash3
wget \
--directory-prefix ${raw_pod5_dir}/pre_wash3 \
--recursive \
--no-check-certificate \
--continue \
--cut-dirs 6 \
--no-host-directories \
--no-parent \
--quiet \
--accept ${pod5_pattern} ${raw_pod5_url1}

# post-wash3
wget \
--directory-prefix ${raw_pod5_dir}/post_wash3 \
--recursive \
--no-check-certificate \
--continue \
--cut-dirs 6 \
--no-host-directories \
--no-parent \
--quiet \
--accept ${pod5_pattern} ${raw_pod5_url2}
```

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

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

## Verify raw read checksums
```{r verify-raw-read-checksums, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

wget \
--directory-prefix ${raw_pod5_dir}/pre_wash3 \
--recursive \
--no-check-certificate \
--continue \
--cut-dirs 6 \
--no-host-directories \
--no-parent \
--quiet \
--accept 'checksums.md5' ${raw_pod5_url1}

cd "${raw_pod5_dir}/pre_wash3"

# Recursively verify checksums in all subdirectories
find $(pwd) -type d | while read -r DIR; do
    # Check if checksums.md5 exists in this directory
    if [[ -f "$DIR/checksums.md5" ]]; then
        echo "Verifying checksums in $DIR"
        (
            cd "$DIR" || exit 1
            md5sum -c checksums.md5
        )
        echo ""
    fi
done
```

```{r verify-raw-read-checksums, engine='bash', eval=FALSE}
# Load bash variables into memory
source .bashvars

wget \
--directory-prefix ${raw_pod5_dir}/post_wash3 \
--recursive \
--no-check-certificate \
--continue \
--cut-dirs 6 \
--no-host-directories \
--no-parent \
--quiet \
--accept 'checksums.md5' ${raw_pod5_url2}

cd "${raw_pod5_dir}/pre_wash3"

# Recursively verify checksums in all subdirectories
find $(pwd) -type d | while read -r DIR; do
    # Check if checksums.md5 exists in this directory
    if [[ -f "$DIR/checksums.md5" ]]; then
        echo "Verifying checksums in $DIR"
        (
            cd "$DIR" || exit 1
            md5sum -c checksums.md5
        )
        echo ""
    fi
done
```

# Trimming and Demultiplexing

We want to demultiplex and trim adapter and barcode sequences, then basecall including modified bases

Notes:

- When selecting multiple modification models, only one modification per canonical base may be active at once. For example, sup,4mC,5mC is invalid as both modification models operate on the C canonical base context.

- Since basecalling is GPU-intensive, but alignment is CPU-intensive, including alignment during basecalling shouldn't have a significant impact on throughput

- Dorado uses `minimap2` for alignment by default, which seems to be more of a k-mer-based pseudoalignment

- When using alignment during basecalling, no reads are excluded following alignment. All reads are basecalled. Then, reads that do not align to the provided reference are included in the output SAM/BAM as *unmapped* reads, with the appropriate tags. If you want to exclude unmapped reads, you'll have to do that *after* basecalling

```{bash}
source .bashvars

# --trim 'all': will trim adapter, primer (if present) and barcode sequences
# --kit-name <kit-name>: allows simultaneous demultiplexing, with which the trimming does not interfere
# --reference: allows basecalling with alignment
# --modified-bases: allows for modified basecalling, space-separated. Note that, for cytosine modifications, there are options for calling modifications on *all* cytosines, and for calling modifications that occur *specifically* in a CG context. Cannot specify modifications that occur on same base (e.g., cannot simultaneously call 4mC and 5mCG)

# --trim: 'all', trim adapter and barcodes
# --kit-name: specifying the library prep i used, Native Barcoding V14 96
# --reference: aligning to Eunicea knighti genome (closest published relative of Eunicea tourneforti)
# -- modified-bases: 5mCG_5hmCG calls 5mC and 5hmC occuring in CG contexts; 6mA calles 6mA modifications.

$dorado basecaller \
hac \
-r ${raw_pod5_dir}/ \
--kit-name SQK-NBD114-96 \
--trim 'all' \
--reference ${genome_dir}/GCA_965233905.1_jaEunKnig1.1_genomic.fna \
--modified-bases 5mCG_5hmCG 6mA \
> ${output_dir_top}/FBD08455_pass_recalled.bam

```


Started above basecalling run 2025/10/28, ~7:30pm.
At 2025/10/29 12:00pm, ~15% complete






Separate out mapped reads, and separate by barcode

```{r, engine='bash', eval=TRUE}
source .bashvars

cd ${output_dir_top}

# Alignment summary
$samtools flagstat FBD08455_pass_recalled.bam > alignment_summary.txt
head alignment_summary.txt

echo ""

# Summarize by barcode
$samtools view FBD08455_pass_recalled.bam | \
awk '{for(i=12;i<=NF;i++) if($i ~ /^BC:Z:/) {print substr($i,6)}}' | \
sort | uniq -c | sort -nr

```

```{bash}
source .bashvars

cd ${output_dir_top}

# Separate out mapped reads
$samtools view -b -F 4 FBD08455_pass_recalled.bam > FBD08455_pass_recalled_mapped.bam

# Separate by barcode
samtools view -h FBD08455_pass_recalled_mapped.bam | awk '/^@/ || /BC:Z:SQK-NBD114-96_barcode27/' | samtools view -b > FBD08455_recalled_mapped_barcode27.bam
samtools view -h FBD08455_pass_recalled_mapped.bam | awk '/^@/ || /BC:Z:SQK-NBD114-96_barcode28/' | samtools view -b > FBD08455_recalled_mapped_barcode28.bam
samtools view -h FBD08455_pass_recalled_mapped.bam | awk '/^@/ || /BC:Z:SQK-NBD114-96_barcode29/' | samtools view -b > FBD08455_recalled_mapped_barcode29.bam


```

## modkit

```{r, engine='bash', eval=TRUE}
source .bashvars
$samtools view -c ${output_dir_top}/FBD08455_pass_recalled.bam
$samtools view -c ${output_dir_top}/FBD08455_pass_recalled_mapped.bam

$modkit summary ${output_dir_top}/FBD08455_pass_recalled.bam
$modkit summary ${output_dir_top}/FBD08455_pass_recalled_mapped.bam

```

