Hello helpful methylation folk! I am assessing the mean read depth of bam files generated using this pipeline with samtools depth and then looking at the man depth of the bedfiles produced by MethylDackel and I see disconcordance. (Even when looking specifically at the cpg sites in my bedfile for the specific sample) Can anyone confirm that this is expected behaviour and why this happens? Note this is for EMseq (Twist) data and I am using the following MethylDackel command:

MethylDackel extract \
        --minDepth 10 \
        -l "$BED_FILE" \
        --OT "$ot_trim" \
        --OB "$ob_trim" \
        "$REFERENCE" \
        "$original_bam" \
        -o "$output_prefix"

Samtools depth command:

samtools depth -a -b "${CpG_BED}" "${BAM_FILE}"

I am no expert in MethylDackel but I wonder if 1) there is a mapQ / baseQ filter in there or 2) the trimming also removes nucleotides form the coverage and not just the methylation calling.

Hi @Eduard Casas, Adding the -s flag with samtools depth I got ~*35% reduction* in average depth, meaning roughly a third of coverage at CpG sites comes from overlapping read pairs - but there was still a discordance with my bedgraph files... To assess affect of quality filtering this I set the -p (phred) and -q (MAPQ) option to 0 instead of default 5 and 10 respectively. This hardly changed anything. Another difference between bam and bedgraph is trimming (--OT/OB). Trimming appears to have a small effect on the mean by 5.5% BUT then I noticed the bam depth counts were almost double that of bedgraphs.. The reason for this is at a CpG, MethylDackel only counts: • Reads on the forward strand looking at the forward C • Reads on the reverse strand looking at the reverse C (the G position on forward) Samtools depth, however, counts all reads overlapping that position regardless of strand, hence the number is double that of the Methydackle values.

# Count forward strand reads at CpG positions (-F 16 excludes reverse strand)
   TMP_FWD="${OUTPUT_DIR}/tmp_fwd_${i}.txt"
   samtools view -F 16 -c -L "${CpG_BED}" "${BAM_FILE}" > "${TMP_FWD}" &
   FWD_PID=$!

   # Count reverse strand reads at CpG positions (-f 16 includes only reverse strand)
   TMP_REV="${OUTPUT_DIR}/tmp_rev_${i}.txt"
   samtools view -f 16 -c -L "${CpG_BED}" "${BAM_FILE}" > "${TMP_REV}" &
   REV_PID=$!