nanoporetech
clarification on apparent discrepancy between snp.vcf and bedmethyl outputs
Hello, after running the epi2me human-variation workflow with --mod and --phased, in the snp.vcf file I see this location where the individual is homozygous alternate (C>T). The relevant vcf line is
chr1 10931 . C T 14.89 PASS ,,, GT:GQ:DP:AD:AF 1/1:14:11:1,10:0.9091
Visual inspection of the cram file in IGV confirms the presence of 10 reads with a T and 1 with a C at this location.
In the bedmethyl merged file (from mods.1, mods.2 and the ungrouped), at this location I see
chr1 10930 10931 h 4 . 10930 10931 255,0,0 4 0.00 0 4 0 0 2 0 0 chr1 10930 10931 m 4 . 10930 10931 255,0,0 4 0.00 0 4 0 0 2 0 0
Thus, N_canonical=4 and N_diff=0. These two outputs don't seem to match up, maybe I'm not interpreting correctly. Can somebody clarify? Thank you
Hello @e-manduchi,
It looks like these bedMethyl records have . as the strand, meaning that Modkit pileup was run with --combine-strands or --preset traditional. What that means is you're summing the counts from the (+)-strand C and the (-)-strand C. You should be seeing N_diff >= 10, however. Could you share the IGV shot with me and/or a BAM with this small region?
Thank you for your response. Indeed I was expecting an N_diff of at least 10 in the bedmethyl file. Here is the IGV screenshot for the region containing the C base of interest, at chr1:10931)
In the above screenshot, reads are arranged by alignment quality. How would you like me to share the bam file with you @ArtRand?
@e-manduchi If you're willing to send over the BAM for this small window that would be helpful in debugging, thank you!
email it to art.rand[at]nanoporetech.com please.
P.S. This issue was also posted on https://github.com/epi2me-labs/wf-human-variation/issues/283.