paired-end mode calling

[avilella]

Hi,

There is a new sample prep that generates read pairs where the CpG
modification needs to be read for both reads in a pair. E.g.

ff = forward/forward read pairs
rr = reverse/reverse read pairs

RefC:Base1Base2 RefG:Base1Base2

C:..  G:A.  ff  =>  ModA
C:,t  G:,,  rr  =>  ModA
C:.T  G:A.  ff  =>  C
C:t,  G:,a  rr  =>  C
C:..  G:..  ff  =>  ModB
C:,,  G:,,  rr  =>  ModB

I couldn't find any way to obtain this currently in MethylDackel,
maybe there is...

Since it seems that MethylDackel is being adopted more and more in the
community, I thought I'd propose to have this new sample prep
implemented in it (I wouldn't want to hack this together in the
"other" methylation caller, with the name starting with a b and ends
with a k, and be stuck with not being able to align the reads with
bwa-meth).

So far I scripted this in a very slow way:

• Take coordinate sorted bam and sort it by read name, so pairs of
  reads are consecutive.

• Create a mini-sam entry for each pair of reads, run samtools
  mpileup, intersect it with the list of CpGs, and record what type is
  each CpG for this read pair in a tabular file.

• Collate the tabular file for summary statistics, per-CpG
  percentages, etc.

On a scale of 1 to 10, how much different would this be from the
current codebase and how feasible to implement?

[dpryan79]

Is this method published anywhere? I'd be interested in reading exactly how all of the library prep is working.

The ModA and ModB case should be easy to implement, but I'm not sure I understand enough of what's going on in the C case to say how hard that would be to implement.

[dpryan79]

Interesting. Feel free to email me directly so I can keep everything confidential.