<p>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).</p>

<p>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.</p>

<p>As an aside, you might find <a href="https://github.com/dpryan79/MethylDackel" rel="nofollow noreferrer">MethylDackel</a> 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.</p>

CC BY-SA 3.0
2017-07-12T17:31:59.380Z
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)](https://www.bioinformatics.babraham.ac.uk/projects/bismark/Bismark_User_Guide.pdf). 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
James Hawley
CC BY-SA 3.0
2017-07-12T17:57:09.170Z
Yes, you can do that for CpGs.
Devon Ryan
CC BY-SA 3.0
2017-07-11T22:22:01.657Z
2017-07-11T22:22:01.657Z
Devon Ryan
<p>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 <code>R</code> with <code>Rsamtools</code> and <code>GenomicRanges</code>, I get something like this:</p>

<pre><code>&gt; reads
GRanges object with 2 ranges and 5 metadata columns:
      seqnames         ranges strand |                     seq                                   XM          XR          XG        NM
         &lt;Rle&gt;      &lt;IRanges&gt;  &lt;Rle&gt; |          &lt;DNAStringSet&gt;                          &lt;character&gt; &lt;character&gt; &lt;character&gt; &lt;integer&gt;
  [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
</code></pre>

<p>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?</p>

<p>For example, the <code>z</code> in the methylation string of the - strand is one position ahead of the <code>z</code> 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?</p>

<r><samtools><methylation>
How to interpret methylation calls from Bismark on opposite strands?
# How to interpret methylation calls from Bismark on opposite strands?

(Asked by: username_0 on 2017-07-11)

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
                  |                                       
  [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?

**Tags:** r, samtools, methylation

---

## 1 answers

### Answer by username_1 on 2017-07-11. Score: 4 | Accepted answer

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](https://github.com/dpryan79/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.

> **Comments:**
> **username_0**: 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)](https://www.bioinformatics.babraham.ac.uk/projects/bismark/Bismark_User_Guide.pdf). 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  
> **username_1**: Yes, you can do that for CpGs.  

---

