Count C versus read length

[EpiBrass]

Hi,

I've had a good look through the bioinformatics thread but as far as I can see nobody has asked this question yet. I am using bisulphite sequencing in plants (CpG, CHG & CHH contexts).

I need to compare different sequencing runs post bisulphite conversion to determine the reproducibility between runs and library preps of the same sample (technical reps). It has been suggested that one way to approach this is to compare the number of C's in a read against the read length post trimming. This can then be graphed and in theory the longer each read is the more C's will be present. If the technical reps are not biased then the trend should be the same for all reps.

Does anyone have any suggestions of how to output from .fastq or .fasta the number Cs in a read and the read length into a tab-delimited file so it can then be graphed?

Many thanks,

Justin

[dpryan]

It sounds like you're describing an "M-bias" graph, an example of which is below (this one has bias):


If you create alignments with either bismark or bison you can make those sorts of graphs easily. With bison, it's the "bison_mbias -pdf something.bam" command, which requires that you have R and the ggplot2 package installed (otherwise, drop the "-pdf" part and it'll just output a tab-separated file that you can use for graphing). Similarly, bismark will produce them but I don't recall off-hand what the command is. I think bismark uses perl GD::Graph.

There's also a recent paper in Bioinformatics on BSeqQC or something like that that can also produce M-bias graphs.

I should not that these graphs aren't looking so much at reproducibility as variance in library setup (bisulfite conversion efficiency, etc.). If you really want to look at reproducibility, it would make more sense to directly compare the methylation calls and calculate how many of them are significantly different (it should by approximately whatever alpha level you use, so ~5% of sites for p<0.05 threshold).