BWA alignement coverage and -n option

Question

[Fabien.delapierre]

Hi, i'm currently comparing two dataset obtained from RAD sequencing (illumina). One is non-treated ( index) and the other one is treated with bisulfit. I'm using bwa to aligne the bisulfit reads against my index. I'm hoping to use samtools to perform SNP calling and identify loci with a site of bisulfite transformation ( methylated C -> T). The sequences are about 100 nt long but the alignement performed by bwa is quite strange to me.

When looking at the coverage of the reads against the index, i found that most of the alignements are only 20 nt long and show a tendency to be aligned between the 25 and 75 nt of the index loci. I can explained by the fact that the bisulfite treatement have "attacked" the edge of the DNA, making them unrecongizable but my question is "how bwa can perform such an alignement" :

I've used the command bwa aln with -n 10, thinking that it will only align bisulfit reads with a maximum of 10 differences (about 10%). How can i get then a coverage of only 20 nt (mean, some of them are close to 100) per loci (~ 80 differences) with this command ?

Thank you for your answer,

Fabien

Answer 1

[Devon Ryan]

Please don't use a standard aligner to handle bisulfite-converted reads, your results will be terrible.

If you're using BWA for other things anyway, then perhaps Brent Pederson's bwa-meth program will be the simplest thing for you to use. There are many other choices out there (such as my own package), so if there's something you don't like about bwa-meth you can search around for other options.