To create pseudo-replicates for statistical analysis, splitting the POD5 files is the best approach. This ensures that each replicate is processed independently from basecalling to alignment, which better simulates technical variability and avoids biases introduced by splitting already-aligned BAM files (which might have mapping dependencies).
Here is the strategy:
- Split POD5s: For each sample, split the original
pod5_passdirectory into two subsets (Replicate 1 and Replicate 2).To create pseudo-replicates for statistical analysis, splitting the POD5 files is the best approach. This ensures that each replicate is processed independently from basecalling to alignment, which better simulates technical variability and avoids biases introduced by splitting already-aligned BAM files (which might have mapping dependencies).
Here is the strategy:
- Split POD5s: For each sample, split the original
pod5_passdirectory into two subsets (Replicate 1 and Replicate 2). We will use a simple script to randomly assign ~50% of reads to each replicate. - Adapt
generate_mapped_modbam.sh: Modify the script to process these split POD5 directories. - Update Samplesheets: Create new samplesheets that list the pseudo-replicates.
Step 1: Script to Split POD5 Files
First, let’s create a helper script to split your POD5 files. Save this as split_pod5.sh.
#!/bin/bash
# split_pod5.sh - Splits POD5 files into two pseudo-replicates
set -euo pipefail
# Configuration
BASE_DIR="/home/jhuang/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans"
POD5_BASE="${BASE_DIR}/X101SC26036392-Z01-J004/Release-X101SC26036392-Z01-J004-20260625_01/Data-X101SC26036392-Z01-J004"
# List of sample folders relative to POD5_BASE
SAMPLES=("S2_Light" "S2_Dark" "T" "O" "WT_Trans" "O_Trans" "WT")
for SAMPLE in "${SAMPLES[@]}"; do
echo "Processing ${SAMPLE}..."
# Define source directory
SRC_DIR="${POD5_BASE}/${SAMPLE}/1732_2B_PBK77125_17f5982a/${SAMPLE}_pod5_pass"
# Define output directories for replicates
REP1_DIR="${SRC_DIR}_rep1"
REP2_DIR="${SRC_DIR}_rep2"
# Create output dirs if they don't exist
mkdir -p "${REP1_DIR}" "${REP2_DIR}"
# Get all pod5 files
FILES=("${SRC_DIR}"/*.pod5)
TOTAL=${#FILES[@]}
if [ "$TOTAL" -eq 0 ]; then
echo "No POD5 files found for ${SAMPLE}. Skipping."
continue
fi
echo "Found ${TOTAL} POD5 files. Splitting..."
# Shuffle and split
# Note: shuf is part of coreutils on Linux
SHUFFLED=($(shuf -e "${FILES[@]}"))
MID=$((TOTAL / 2))
# Copy first half to rep1, second half to rep2
for i in "${!SHUFFLED[@]}"; do
FILE="${SHUFFLED[$i]}"
BASENAME=$(basename "${FILE}")
if [ "$i" -lt "$MID" ]; then
cp "${FILE}" "${REP1_DIR}/${BASENAME}"
else
cp "${FILE}" "${REP2_DIR}/${BASENAME}"
fi
done
echo "✅ ${SAMPLE} split complete: ${MID} files in rep1, $((TOTAL - MID)) files in rep2"
done
echo "🎉 All POD5 files split successfully!"Run this script once:
chmod +x split_pod5.sh
./split_pod5.shStep 2: Adapted generate_mapped_modbam_rep.sh
This updated script now processes the _rep1 and _rep2 directories created above. It generates separate BAM files for each replicate.
#!/bin/bash
#===============================================================================
# generate_mapped_modbam_pseudo_reps_rep.sh - Local Dorado generation for pseudo-replicates
#===============================================================================
set -euo pipefail
# === Configuration ===
DORADO="/home/jhuang/Tools/dorado-2.0.0-linux-x64/bin/dorado" # Your local Dorado path
BASE_DIR="/home/jhuang/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans"
OUTDIR="${BASE_DIR}"
# Reference genome paths
REF_S2_LIGHT="${BASE_DIR}/S2_Light-trycycler-medaka_polished_genome.fa"
REF_S2_DARK="${BASE_DIR}/S2_Dark-trycycler-medaka_polished_genome.fa"
REF_T="${BASE_DIR}/T-trycycler-medaka_polished_genome.fa"
REF_O="${BASE_DIR}/O-trycycler-medaka_polished_genome.fa"
REF_WT_TRANS="${BASE_DIR}/WT_Trans-trycycler-medaka_polished_genome.fa"
REF_O_TRANS="${BASE_DIR}/O_Trans-trycycler-medaka_polished_genome.fa"
REF_WT="${BASE_DIR}/WT-trycycler-medaka_polished_genome.fa"
# POD5 data paths (Updated to point to split replicas)
POD5_BASE="${BASE_DIR}/X101SC26036392-Z01-J004/Release-X101SC26036392-Z01-J004-20260625_01/Data-X101SC26036392-Z01-J004"
# Helper function to get pod5 path
get_pod5_path() {
local SAMPLE=$1
local REP=$2
echo "${POD5_BASE}/${SAMPLE}/1732_2B_PBK77125_17f5982a/${SAMPLE}_pod5_pass_${REP}"
}
# Dorado model
MODEL="dna_r10.4.1_e8.2_400bps_sup@v5.0.0"
# Ensure all reference genomes are indexed
echo "📚 Indexing reference genomes..."
for REF in "${REF_S2_LIGHT}" "${REF_S2_DARK}" "${REF_T}" "${REF_O}" "${REF_WT_TRANS}" "${REF_O_TRANS}" "${REF_WT}"; do
if [ ! -f "${REF}.fai" ]; then
echo " Indexing: $(basename ${REF})"
samtools faidx "${REF}"
fi
done
# === Function: Generate modBAM ===
generate_modbam() {
local SAMPLE_NAME=$1
local REF=$2
local POD5_DIR=$3
local MOD_BASES=$4
echo "🚀 Generating ${SAMPLE_NAME} ${MOD_BASES} modBAM (aligned)..."
# Check if pod5 dir exists
if [ ! -d "${POD5_DIR}" ]; then
echo "❌ Error: POD5 directory not found: ${POD5_DIR}"
return 1
fi
"${DORADO}" basecaller \
--modified-bases "${MOD_BASES}" \
--emit-moves \
--device cuda:0 \
--reference "${REF}" \
"${MODEL}" \
"${POD5_DIR}" | samtools view -b - > "${OUTDIR}/${SAMPLE_NAME}_${MOD_BASES//\//_}_mapped.mod.bam"
}
# === Process all samples and replicates ===
SAMPLES=("S2_Light" "S2_Dark" "T" "O" "WT_Trans" "O_Trans" "WT")
REPS=("rep1" "rep2")
MOD_TYPES=("6mA" "4mC_5mC")
# Associative array for references
declare -A REFS
REFS["S2_Light"]="${REF_S2_LIGHT}"
REFS["S2_Dark"]="${REF_S2_DARK}"
REFS["T"]="${REF_T}"
REFS["O"]="${REF_O}"
REFS["WT_Trans"]="${REF_WT_TRANS}"
REFS["O_Trans"]="${REF_O_TRANS}"
REFS["WT"]="${REF_WT}"
for SAMPLE in "${SAMPLES[@]}"; do
REF="${REFS[$SAMPLE]}"
for REP in "${REPS[@]}"; do
POD5_DIR=$(get_pod5_path "${SAMPLE}" "${REP}")
SAMPLE_TAG="${SAMPLE}_${REP}"
for MOD in "${MOD_TYPES[@]}"; do
generate_modbam "${SAMPLE_TAG}" "${REF}" "${POD5_DIR}" "${MOD}"
done
done
done
# === Verify Output ===
echo ""
echo "🔍 Verifying BAM files..."
for BAM in "${OUTDIR}"/*_rep[12]_*_mapped.mod.bam; do
if [ -f "${BAM}" ]; then
if samtools quickcheck "${BAM}" 2>/dev/null; then
READS=$(samtools view -c "${BAM}")
echo "✅ $(basename "${BAM}"): ${READS} reads"
else
echo "❌ $(basename "${BAM}"): Corrupted or invalid format"
fi
fi
done
echo ""
echo "🎉 Pseudo-replicate modBAM generation complete!"
echo "📁 Output directory: ${OUTDIR}"Step 3: Updated Samplesheets
You need to update your CSV files to include the pseudo-replicates. The group column should remain the same (e.g., “WT”) so that downstream tools know they belong to the same biological condition, but the sample name must be unique (e.g., “WT_rep1”).
samplesheet_6mA_rep.csv
group,sample,path,ref,method
WT,WT_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_rep1_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT-trycycler-medaka_polished_genome.fa,ont
WT,WT_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_rep2_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT-trycycler-medaka_polished_genome.fa,ont
T,T_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/T_rep1_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/T-trycycler-medaka_polished_genome.fa,ont
T,T_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/T_rep2_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/T-trycycler-medaka_polished_genome.fa,ont
O,O_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_rep1_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O-trycycler-medaka_polished_genome.fa,ont
O,O_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_rep2_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O-trycycler-medaka_polished_genome.fa,ont
WT_Trans,WT_Trans_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_Trans_rep1_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_Trans-trycycler-medaka_polished_genome.fa,ont
WT_Trans,WT_Trans_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_Trans_rep2_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_Trans-trycycler-medaka_polished_genome.fa,ont
O_Trans,O_Trans_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_Trans_rep1_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_Trans-trycycler-medaka_polished_genome.fa,ont
O_Trans,O_Trans_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_Trans_rep2_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_Trans-trycycler-medaka_polished_genome.fa,ont
S2_Light,S2_Light_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Light_rep1_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Light-trycycler-medaka_polished_genome.fa,ont
S2_Light,S2_Light_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Light_rep2_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Light-trycycler-medaka_polished_genome.fa,ont
S2_Dark,S2_Dark_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Dark_rep1_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Dark-trycycler-medaka_polished_genome.fa,ont
S2_Dark,S2_Dark_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Dark_rep2_6mA_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Dark-trycycler-medaka_polished_genome.fa,ontsamplesheet_4mC_5mC_rep.csv
group,sample,path,ref,method
WT,WT_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_rep1_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT-trycycler-medaka_polished_genome.fa,ont
WT,WT_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_rep2_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT-trycycler-medaka_polished_genome.fa,ont
T,T_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/T_rep1_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/T-trycycler-medaka_polished_genome.fa,ont
T,T_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/T_rep2_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/T-trycycler-medaka_polished_genome.fa,ont
O,O_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_rep1_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O-trycycler-medaka_polished_genome.fa,ont
O,O_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_rep2_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O-trycycler-medaka_polished_genome.fa,ont
WT_Trans,WT_Trans_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_Trans_rep1_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_Trans-trycycler-medaka_polished_genome.fa,ont
WT_Trans,WT_Trans_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_Trans_rep2_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/WT_Trans-trycycler-medaka_polished_genome.fa,ont
O_Trans,O_Trans_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_Trans_rep1_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_Trans-trycycler-medaka_polished_genome.fa,ont
O_Trans,O_Trans_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_Trans_rep2_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/O_Trans-trycycler-medaka_polished_genome.fa,ont
S2_Light,S2_Light_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Light_rep1_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Light-trycycler-medaka_polished_genome.fa,ont
S2_Light,S2_Light_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Light_rep2_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Light-trycycler-medaka_polished_genome.fa,ont
S2_Dark,S2_Dark_rep1,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Dark_rep1_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Dark-trycycler-medaka_polished_genome.fa,ont
S2_Dark,S2_Dark_rep2,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Dark_rep2_4mC_5mC_mapped.mod.bam,/mnt/md1/DATA/Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans/S2_Dark-trycycler-medaka_polished_genome.fa,ontNext Steps
- Run
split_pod5.shto create the replicated POD5 directories. - Run the new
generate_mapped_modbam.shto create the BAM files. Note: This will take roughly twice as long as before since you are processing all reads again. - Use the new CSV files with
nf-core/methylong. - When running
modkit pileupmanually afterwards, ensure you update the loop to iterate over the new replicate names (e.g.,WT_rep1,WT_rep2, etc.) if you wish to generate per-replicate BED files, or keep the original logic ifnf-core/methylonghandles the aggregation correctly. Usually, for differential methylation analysis later, having per-replicate BEDs or counts is beneficial.
