---
title: "08-E5-Deep-Dive-Expression-WGBS"
author: "Zoe Dellaert"
date: "2024-10-01"
always_allow_html: true
output: 
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: true 
---

## E5 Deep Dive Expression WGBS

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

# Create a Bash variables file

This allows usage of Bash variables (e.g. paths to common directories) across R Markdown chunks.
```{bash bash-vars}
{
echo "#### Assign Variables ####"
echo ""

echo "# Data directories"
echo 'export repo_dir=/work/pi_hputnam_uri_edu/deep-dive-expression'
echo 'export code_dir=${repo_dir}/D-Apul/code'
echo 'export genome_dir=${repo_dir}/D-Apul/data/genome'
echo 'export output_dir_top=${repo_dir}/D-Apul/output/08-Apul-WGBS'
echo 'export sim_data_dir=${repo_dir}/D-Apul/data/08-Apul-WGBS/sim-data'
echo ""

echo "# Paths to programs"
echo 'export sherman=/work/pi_hputnam_uri_edu/pgrams/Sherman/Sherman'
echo ""

echo "# Programs associative array"
echo "declare -A programs_array"
echo "programs_array=("
echo '[sherman]="${sherman}"'
echo ")"
} > .bashvars

cat .bashvars
```

## Genome

```{bash, eval=FALSE}
cd ../data/genome

wget -O Apulcra-genome.fa "https://osf.io/download/kn96u/"
```

## simulate WGBS data for *Acropora pulchra* on Unity

#### Install Sherman on unity in shared directory and make executable

```{bash install-sherman, eval=FALSE}
cd /work/pi_hputnam_uri_edu/pgrams

git clone https://github.com/FelixKrueger/Sherman.git

cd Sherman
chmod +x Sherman
```

```{bash run-sherman, eval=FALSE}
source .bashvars

${programs_array[sherman]} \
--genome ${genome_dir} \
--length 150 \
--number_of_seqs 1000000 \
--paired_end \
--output_dir ${sim_data_dir}

gzip ${sim_data_dir}/*.fastq
```

### Run nfcore-methylseq based on [Danielle's Notebook](https://github.com/daniellembecker/DanielleBecker_Lab_Notebook/blob/master/_posts/2023-03-13-P.meandrina-WGBS-Workflow-Host.md)

Create input csv file containing the simulated data 

```{bash}
source .bashvars

echo "sample,fastq_1,fastq_2" > ${sim_data_dir}/input.csv
echo "SIM_1,${sim_data_dir}/simulated_1.fastq.gz,${sim_data_dir}/simulated_2.fastq.gz" >> ${sim_data_dir}/input.csv

# Optional: Display the file content to verify
cat ${sim_data_dir}/input.csv
```

```{bash}
# Create the job script
cat << 'EOF' > methylseq.sh
#!/usr/bin/env bash
#SBATCH --export=NONE
#SBATCH --nodes=1 --ntasks-per-node=24
#SBATCH --signal=2
#SBATCH --no-requeue
#SBATCH --mem=500GB
#SBATCH --error="%x_error.%j"
#SBATCH --output="%x_output.%j"
#SBATCH -t 24:00:00
#SBATCH --mail-type=BEGIN,END,FAIL

# Load modules
module load nextflow/24.04.2

# Source the bashvars file
source .bashvars

# Navigate to the output directory
cd ${output_dir_top}

# Print variables from bashvars
echo "Simulated data path: ${sim_data_dir}"

# Run nextflow methylseq

nextflow run nf-core/methylseq -profile singularity \
--aligner bismark \
--fasta ${genome_dir}/Apulcra-genome.fa \
--save_reference \
--input ${sim_data_dir}/input.csv \
--zymo \
--cytosine_report \
--relax_mismatches \
--unmapped \
--outdir ${output_dir_top}/sim_Apul_methylseq
EOF

chmod +x methylseq.sh

# Submit the job to the cluster , go to terminal to do this - it gets confusing in Rmd...
# sbatch methylseq.sh
```

### code for running nf-core/methylseq on Andromeda since it is maybe not installed on unity...

```{bash, eval=FALSE}
mkdir scripts
nano scripts/methylseq.sh
```

```{bash, eval=FALSE}
#!/bin/bash
#SBATCH -t 120:00:00
#SBATCH --nodes=1 --ntasks-per-node=1
#SBATCH --cpus-per-task=10
#SBATCH --export=NONE
#SBATCH --mem=500GB
#SBATCH --mail-type=BEGIN,END,FAIL #email you when job starts, stops and/or fails
#SBATCH --mail-user=zdellaert@uri.edu #your email to send notifications
#SBATCH -D /data/putnamlab/zdellaert/deep-dive-exp-WGBS/

# load modules needed

module load Nextflow/22.10.1

# run nextflow methylseq

nextflow run nf-core/methylseq -profile singularity \
--aligner bismark \
--fasta /data/putnamlab/REFS/Apul/apul.hifiasm.s55_pa.p_ctg.fa.k32.w100.z1000.ntLink.5rounds.masked.fa \
--save_reference \
--input 'input.csv' \
--zymo \
--trim_poly_g \
--cytosine_report \
--relax_mismatches \
--unmapped \
--outdir sim_Apul_WGBS
```

```{bash, eval=FALSE}
sbatch scripts/methylseq.sh
```

