BS-seq problems – the solution

In the previous post I described a problem with BS-seq data. After much messing about I found the solution in the analysis tool methylKit. Rather than using the final output of the Bismark aligner tool (the ‘coverage’ or ‘cytosineReport’ files) you can use methylKit to process the BAM file Bismark makes in order to generate a file that can be re-read into methylKit for analysis. This function, processBismarkAln, includes the option to insist that each read base covering a particular position be of a specified PHRED quality. As the issue described in the previous post is characterised by poor quality bases, setting this to 30+ excludes the problematic bases from the analysis without any need to chop away good data using standard QC tools.