4

I am interested in identifying indels in whole genome bisulfite sequencing data (76bp paired end). Currently, I do this by setting the -rfg and -rdg affine gap penalty scores for bowtie2 to more permissive values than the default 5+3N and mapping using the bisulfite sequencing alignment wrapper, bismark.

My question is, what values of -rfg and -rdg will allow me to identify longest possible indels without sacrificing alignment quality? Is it better to set the affine penalty to zero with a high penalty for initially opening an indel (ex. 8+0N)? Or is it better to keep the initial penalty low and having a nonzero penalty for extension (ex. 2+1N)?

Ben D.
  • 397
  • 1
  • 10

1 Answers1

3

You'd be best off by starting with -rfg and -rdg as is and reverting bismark's change of --score-min back to the default for bowtie2. That alone will allow for much longer indels. If that still doesn't suffice, then I'd play around more with --score-min before messing with the gap open/extend penalties. If you do need to play with those, then increase the gap open penalty and decrease the gap extension to 1. Do check that the resulting alignments aren't nonsense though!

Devon Ryan
  • 19,602
  • 2
  • 29
  • 60
  • With the bismark default --score-min, -rdg, and -rfg settings, it appears that the maximum sized indel identifiable with my read length is 3bp. I changed both -rdg and -rfg to a gap open/extension of (1,1) and the maximum length I can detect is 14bp. I will try modifying --score-min only. It's a fine line to walk, detecting the longest possible insertions while keeping the alignment clean! – Ben D. Aug 04 '17 at 03:38