Generating Pseudo-Replicates for ONT Methylation Analysis: Splitting POD5 Files vs. BAM Subsampling for nf-core/methylong (Data_Tam_Methylation_2026_WT_T_O_T_Trans_O_Trans)

Leave a reply

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:

  1. Split POD5s: For each sample, split the original pod5_pass directory 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:

  1. Split POD5s: For each sample, split the original pod5_pass directory into two subsets (Replicate 1 and Replicate 2). We will use a simple script to randomly assign ~50% of reads to each replicate.
  2. Adapt generate_mapped_modbam.sh: Modify the script to process these split POD5 directories.
  3. 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.sh

Step 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,ont

samplesheet_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,ont

Next Steps

  1. Run split_pod5.sh to create the replicated POD5 directories.
  2. Run the new generate_mapped_modbam.sh to create the BAM files. Note: This will take roughly twice as long as before since you are processing all reads again.
  3. Use the new CSV files with nf-core/methylong.
  4. When running modkit pileup manually 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 if nf-core/methylong handles the aggregation correctly. Usually, for differential methylation analysis later, having per-replicate BEDs or counts is beneficial.

Leave a Reply

Your email address will not be published. Required fields are marked *