I'm looking at Reduced representation bisulfite sequencing (RRBS) data from ENCODE, and to align the FASTQ files I've used Bismark with Bowtie 1. When I load the resulting BAM file into R with Rsamtools and GenomicRanges, I get something like this:

> reads
GRanges object with 2 ranges and 5 metadata columns:
      seqnames         ranges strand |                     seq                                   XM          XR          XG        NM
         <Rle>      <IRanges>  <Rle> |          <DNAStringSet>                          <character> <character> <character> <integer>
  [1]     chrM [16523, 16558]      - | ATAAAACCTA...CCCTTAAATA .....h........h.......z.............          CT          GA         3
  [2]     chrM [16524, 16559]      + | TAAAGTTTAA...TTTTAAATAA .....hh.......hhh.h.z...hhhh........          CT          CT        11
  seqinfo: 25 sequences from an unspecified genome

For the read on the - strand (which is aligned to the "G->A"-converted reference), how should I compare methylation calls to those on the + strand, since they don't line up?

For example, the z in the methylation string of the - strand is one position ahead of the z in the + strand (which makes sense because of symmetric CpG methylation). But how should I determine whether these two methylation calls are "essentially the same" or not?


1 Answer 1


The strand the bismark reports is related to the strand from which the read originated, not necessarily how it's aligned. So, alignments on the + strand shouldn't have calls overlapping those on the - strand, since you can't have a C in the same place on the same strand. One should often see Z/z next to each other on opposite strand, like in your example, since these are CpG (so it should be a Z/z on the + strand and then a Z/z on the following base on the - strand). Thus, in your example you have two reads supporting unmethylation for the CpG as a whole (1 read for each of the Cs).

The most confusing thing about BSseq is that reads (or paired-end reads) are only ever informative for a single strand. This is actually true for all sequencing with Illumina instruments (single-stranded fragments are loaded, after all), but since strands are almost always complementary we can usually ignore this fact.

As an aside, you might find MethylDackel useful (full disclosure, I wrote it). It'll be much faster at extracting methylation calls than bismark and supports nice things like excluding regions of bias methylation and likely variant positions.

  • $\begingroup$ OK, that's what I wasn't sure about, because of the description given by the Babraham Institute about how Bismark works (pg. 5, for anyone who's curious). So not having the Z/z's overlap between reads on opposite strands makes sense. So to say that "these reads are concordant in their methylation calls", would I just shift the positions that I compare between + and - strand reads by 1 to properly align them? Only for CpG methylation, not CHH- or CHG-context methylation $\endgroup$ Commented Jul 12, 2017 at 17:31
  • $\begingroup$ Yes, you can do that for CpGs. $\endgroup$
    – Devon Ryan
    Commented Jul 12, 2017 at 17:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.