---
title: "02-gallo-align"
format: html
editor: visual
---

```{bash}

# Directories and programs
bismark_dir="/home/shared/Bismark-0.24.0"
bowtie2_dir="/home/shared/bowtie2-2.4.4-linux-x86_64/"
genome_folder="../data/"
reads_dir="../data/Raw_WGBS_Mch/"

 ${bismark_dir}/bismark \
--path_to_bowtie ${bowtie2_dir} \
-genome ${genome_folder} \
-p 4 \
-score_min L,0,-0.6 \
--non_directional \
-1 ${reads_dir}LCo_BSr1_R1.fastq \
-2 ${reads_dir}LCo_BSr1_R2.fastq \
-o ../output/02-gallo-align/gallo-out
```

```{bash}
# Directories and programs
bismark_dir="/home/shared/Bismark-0.24.0"
bowtie2_dir="/home/shared/bowtie2-2.4.4-linux-x86_64/"
genome_folder="../data/"
reads_dir="../data/Raw_WGBS_Mch/"
output_base="../output/02-gallo-align/gallo-out2"

# Score_min parameters to test
score_min_params=(
    "L,0,-0.4"
    "L,0,-0.6"
    "L,0,-0.8"
    "L,0,-1.0"
    "L,-1,-0.6"
)

# Loop through parameters
for score in "${score_min_params[@]}"; do
    score_tag=$(echo $score | sed 's/,/_/g')
    output_dir="${output_base}_${score_tag}"
    log_file="${output_dir}/bismark_run.log"

    mkdir -p "$output_dir"

    echo "Running Bismark with score_min=$score" | tee "$log_file"

    ${bismark_dir}/bismark \
        --path_to_bowtie ${bowtie2_dir} \
        -genome ${genome_folder} \
        -p 4 \
        --score_min "$score" \
        -u 100000 \
        -1 ${reads_dir}LCo_BSr1_R1.fastq \
        -2 ${reads_dir}LCo_BSr1_R2.fastq \
        -o "$output_dir" \
        >> "$log_file" 2>&1
done
```

```{bash}
#!/bin/bash

# Output base
base_dir="../output/02-gallo-align"

# Header for CSV
echo "score_min,total_reads,unique_alignments,ambiguous_alignments,unmapped_reads,mapping_efficiency,cpg_meth,chg_meth,chh_meth,unknown_meth" > ../output/02-gallo-align/merged_alignment_summary.csv

# Loop through each score_min folder
for dir in ${base_dir}/gallo-out2_*; do
    # Find report file
    report=$(find "$dir" -type f -name "*report.txt" | head -n 1)

    if [[ -f "$report" ]]; then
        score_min=$(basename "$dir" | sed 's/gallo-out2_//')

        total_reads=$(grep "Sequence pairs analysed" "$report" | awk '{print $4}')
        unique=$(grep "Number of paired-end alignments with a unique best hit" "$report" | awk '{print $10}')
        ambiguous=$(grep "Number of paired-end alignments with multiple best hits" "$report" | awk '{print $10}')
        unmapped=$(grep "Number of read pairs that could not be aligned" "$report" | awk '{print $11}')
        efficiency=$(grep "Mapping efficiency" "$report" | awk '{print $3}' | sed 's/%//')

        # Methylation levels
        cpg=$(grep "C methylated in CpG context" "$report" | awk '{print $6}' | sed 's/%//')
        chg=$(grep "C methylated in CHG context" "$report" | awk '{print $6}' | sed 's/%//')
        chh=$(grep "C methylated in CHH context" "$report" | awk '{print $6}' | sed 's/%//')
        unknown=$(grep "C methylated in unknown context" "$report" | awk '{print $8}' | sed 's/%//')

        echo "${score_min},${total_reads},${unique},${ambiguous},${unmapped},${efficiency},${cpg},${chg},${chh},${unknown}" >> ../output/02-gallo-align/merged_alignment_summary.csv
    else
        echo "⚠️ No report found in $dir"
    fi
done
```

```{r}

library(tidyverse)  # includes readr, dplyr, ggplot2, etc.
summary_df <- read_csv("../output/02-gallo-align/merged_alignment_summary.csv")
glimpse(summary_df)
View(summary_df)
```



```{bash}

# Directories and programs
bismark_dir="/home/shared/Bismark-0.24.0"
bowtie2_dir="/home/shared/bowtie2-2.4.4-linux-x86_64/"
genome_folder="../data/"
reads_dir="../data/Raw_WGBS_Mch/"

cd ../output/02-gallo-align/gallo-out/

find *.bam | \
xargs basename -s .bam | \
xargs -I{} ${bismark_dir}/deduplicate_bismark \
--bam \
--paired \
{}.bam 

```


```{bash}
# Directories and programs
bismark_dir="/home/shared/Bismark-0.24.0"

cd ../output/02-gallo-align/gallo-out/

${bismark_dir}/bismark_methylation_extractor \
--bedGraph \
--counts \
--comprehensive \
--merge_non_CpG \
--multicore 42 \
--buffer_size 75% \
*deduplicated.bam
```

```{bash}

# Directories and programs
bismark_dir="/home/shared/Bismark-0.24.0"
bowtie2_dir="/home/shared/bowtie2-2.4.4-linux-x86_64/"
genome_folder="../data/"
reads_dir="../data/Raw_WGBS_Mch/"

cd ../output/02-gallo-align/gallo-out/

find *deduplicated.bismark.cov.gz \
| xargs basename -s _trimmed_bismark_bt2.deduplicated.bismark.cov.gz \
| xargs -I{} ${bismark_dir}/coverage2cytosine \
--genome_folder ${genome_folder} \
-o {} \
--merge_CpG \
--zero_based \
{}_bismark_bt2.deduplicated.bismark.cov.gz
```

