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:> readsGRanges 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 genomeFor 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?
How to interpret methylation calls from Bismark on opposite strands?
r;samtools;methylation
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.