Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-12T17:31:59.380Z
Answers.Text: 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
Answers.UserDisplayName: James Hawley
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-12T17:57:09.170Z
Answers.Text: Yes, you can do that for CpGs.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-11T22:22:01.657Z
Answers.LastActivityDate: 2017-07-11T22:22:01.657Z
Answers.UserDisplayName: Devon Ryan
Body: <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>
  
Tags: <r><samtools><methylation>
Title: How to interpret methylation calls from Bismark on opposite strands?
ThreadText: # 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
           <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?
  
  **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.  
  
  ---
  
  