Dear all, I have a question concerning the analysis of targeted bisulfite sequencing data. I am completely new to bisulfite sequencing analysis, so please forgive me if this is a stupid question. We performed bisulfite conversion and target region amplification of several loci on 3 chromosomes using 52 primer pairs in total, followed by paired-end sequencing on an Illumina MiSeq platform. Using the nf-core/methylseq workflow with default trimming parameters (nextflow run nf-core/methylseq --input samplesheet_final.csv --outdir methylseq_final --genome GRCh38 --skip_deduplication true -profile singularity -c slurm.config), alignment rates are very low, around 2%. I tried some things and found that alignment goes up to ~ 70% when trimming 80 bases off the 3’ end using --three_prime_clip_r1/r2, so there must be something at the 3’ end. When checking my sequences after trimming, I noticed that after the reverse primer sequences, there is some potentially junk sequence (beginning with AGACCAAGTCTCT) that is not removed. The automatically recognized and trimmed adapter sequence was the Illumina TruSeq adapter (AGATCGGAAGAGC), but the strange sequence I could find in a lot of reads. My question now is, does anyone have advice on how to proceed here? I thought about manually trimming with cutadapt using the primer sequences, but some reads have the primer sequences in the middle because the amplicons are designed to overlap. Is it possible to align the reads to only the target sequences (starting with the forward primer and ending with the reverse primer) so that the 3’ junk gets ignored? Any help would be very much appreciated! Please let me know if I should paste any workflow logs or you need more information in general. Best wishes, Luisa

Hi Luisa, this is probably less of a Nextflow question but more specific to your data. I’d be happy to help you out with this, but maybe you we can take this offline via email? Some things that would be helpful: • can you send a trimming report as this will show how the auto-detection performed (and how many reads it found. • AGATCGGAAGAGC is not any adapter sequence I know of, it could be biological or one of the amplicons? • could you attach a FastQC report (the entire html report), this might give us clues • amplicon sequencing sometimes gives weird strand alignments, e.g. OT+CTOT, OT+CTOB. Can you try out --non-directional mapping and send the strand stats (or mapping report?) • finally, if you could send a few test reads (e.g. 200K, gzipped, untrimmed) I could take a look

Sorry for coming to the wrong place with this issue, I sent you an email now! Thank you very much! :)

it’s not the wrong place per se, but a good starting point