Hi All, I'm running the methylSeq pipeline with EM-seq data derived from a targeted panel. I've set the filtering parameters but noticed this only produces a filtered bedGraph and not a filtered coverage (

.cov

) file. This isnt ideal as many of the downstream packages for WGBS data analysis prefer coverage files as input. I'm wondering if running something like:

for f in *.cov.gz; do
    echo "Processing $f"
    zcat "$f" | \
    bedtools intersect -a - \
                      -b /Path_to_BED/covered_targets.bed \
                      -wa | \
    gzip > "${f%.cov.gz}.filtered.cov.gz"
done

is suitable to filter the

.cov

files with reference to a targeted regions BED file. Seems like it is but I'm looking for a second opinion. If anyone has done anything similar or has advice on the approach to filtering the

.cov

files, I'd appreciate their input. Thanks in advance!

Hi Ryan, I haven’t down this myself but your approach looks sound. Maybe the workflow should be extended to produce both filtered bedGraph and coverage files (as this generally is the more useful file…) You could test this very easily by importing the original and filtered .cov.gz files into SeqMonk (it has an option Import -> Bismark (.cov)), and checking it out (you could also import your file

covered_targets.bed

as annotations and then run a quantification, but I am pretty sure it will be obvious very quickly (<2 mins). Alternatively, you could verify this on the command line:

# Compare line counts
echo "Original: $(zcat original.cov.gz | wc -l)"
echo "Filtered: $(zcat original.filtered.cov.gz | wc -l)"

# Verify all filtered sites fall within targets
zcat original.filtered.cov.gz | \
    bedtools intersect -a - -b covered_targets.bed -v | \
    wc -l  # should be 0

All the best!

Hi Felix, Thanks for the response and PR. I think the adjustment will help a lot of people out. I wanted to flag a small issue I encountered when comparing the

.targeted.bedGraph

against the

.filtered.cov.gz

outputs which I ran as a sanity check. The filtered

.cov

file retains more CpGs than the targeted bedGraph, and I think this comes down to the coordinate difference between the two files. Assuming the target reference is a bunch of probes (as in my case having used the Twist Human Methylome panel), the 0-based start coordinate can fall outside the target boundary, causing sites to be dropped during bedGraph filtering that are otherwise retained after filtering the coverage file. For example: bedGraph interval: 1 10788 10789 Target reference: 1 10789 10909 cg16269199,cg26928153,cpg_inter missing_from_bed_vs_cov.bed 1 10788 10789 (see code below on this) Would it be worth adding 1 to the bedGraph start position prior to filtering to avoid losing these CpGs? I'm not sure how widespread this issue is, but it should be a straightforward fix unless I'm missing something. Otherwise, thanks for your help.

COV_FILE="Sample_1.filtered.cov.gz"
BED_FILE="Sample_1.targeted.bedGraph"

bedtools intersect \
    -a <(zcat "${COV_FILE}" | awk '{print $1"\t"$2-1"\t"$2}') \
    -b "${BED_FILE}" \
    -v > missing_from_bed.bed
echo "Missing from bed: $(wc -l < missing_from_bed.bed)"
Missing from bed: 19235 (This was the difference in CpG number between my filtered files).

Thanks for the reply. Claude thought the following: The issue with adding 1 to start:

10788 10789

→

10789 10789

creates a zero-width interval in BED convention (

[10789, 10789)

= empty), which would have undefined behavior in bedtools. Better approach — extend the end by 1 to represent the full CpG dinucleotide:

10788 10789

→

10788 10790

This makes the interval cover both the C and the G of the CpG, which correctly overlaps with a target starting at 10789. You’d do the intersection on the extended intervals, then restore the original bedGraph coordinates. In practice, this could be done in the

bedtools_intersect.config

by adding an

awk

pre-processing step, or more cleanly by using

bedtools slop -r 1

before the intersection. In the nf-core context, the cleanest option would probably be to add

ext.args = '-wa'

combined with a pre-processing

awk

in the process config, or to add a

BEDTOOLS_SLOP

step before

BEDTOOLS_INTERSECT

. Alternatively, the simplest one-liner fix: in the bedtools_intersect config, you could shift the logic to use

-a

with the bedGraph but apply

-loj

(left outer join) style filtering... though that changes semantics. The most robust and transparent fix would be: Extend bedGraph entries by 1bp on the right before intersection (

awk '{OFS="\t"; $3=$3+1; print}'

) Run

bedtools intersect

1. Restore original coordinates (

awk '{OFS="\t"; $3=$3-1; print}'

) This way the filtering logic explicitly accounts for the CpG dinucleotide spanning 2 bases. Shall we implement that?

Yeah, whatever is simplest! Thanks again, Ryan

Changes have been pushed to https://github.com/nf-core/methylseq/pull/604, could you give that branch a go to see if it fixes your the issue?

Hi Felix, I ran the PR as it is and it successfully produced the filtered

.cov

files, but the bedGraph files were empty. I left a comment on GitHub regarding line 37. Changing:

- awk 'BEGIN{OFS="\\t"} {\$3=\$3+1; print}' ${bedgraph} > extended.bedGraph
+ zcat ${bedgraph} | awk 'BEGIN{OFS="\t"} {$3=$3+1; print}' > extended.bedGraph

resolved the issue when I re-ran the PR with this adjustment. Eveything else seems to run fine so it should be good to go 👍

good catch, this is now fixed

Re-ran the updated pull request and the outputs all look good. Nothing is empty and the numbers match up. Should work well for the next version of nf-core/methylSeq.