Help with Bismark methylation extractor

[Dipro]

Hi,

Could I please get some advice on the following issue. I am working on targted capture genome-wide DNA methylation assay. The aligner I use is Bismark/bowtie2. After obtaining output from bismark, I isolated targeted regions in the bismark BAM output file by using bedtools 'intersect' between the bismark BAM output file and a targeted region bed file. I then ran the bismark methylation extractor. I obtained this error:

The IDs of Read 1 (HWI-D00356:33:C3D9RACXX:1:1101:1482:2234_1:N:0:TGACCA) and Read 2 (HWI-D00356:33:C3D9RACXX:1:1101:1551:2226_1:N:0:TGACCA) are not the same. This might be a result of sorting the paired-end SAM/BAM files by chromosomal position which is not compatible with correct methylation extraction. Please use an unsorted file instead

I have not done any sorting prior to methylation extractor so I am unsure about the cause of this issue. If I ommit bedtools 'intersect' and used bismark methylation extractor to process the bismark BAM output, the methylation extractor works well. It therefore suggests that something is going on with the bedtools 'intersect' which I am not aware of. Could someone please advice?

Thanks very much.

[fkrueger]

Hi Dipro,

If you run the methyation extractor in paired-end mode it expects read 1 and read 2 of a pair to follow each other line by line in the BAM file. If you disturb this R1/R2 order you would team up two reads that have nothing to do with each other, and thus features such as the --no_overlap detection would not work. Also, this will interfere with the strand origin of the reads. Because of this we included a check at the start that would not allow the methylation extraction of such files since it almost never is what you are intending to do. What appears to have happened with bedtools intersect is that it accepted some hits where only R1 or R2 of a pair overlapped with a desired region, but its partner read didn't, thus getting the read pairs out of sync.

As solutions to this you could:

a) run the methylation extractor first and simply use all data for downstream analysis

b) run the methylation extractor first and then filter on the position of the methylation calls afterwards (even though this might not work out of the box with bedtools intersect if it requires BAM files)

c) write a quick script to go through the filtered BAM file that reads in two lines at a time and sorts R1 and R2 into a paired-end file if they follow each other, and sort singleton reads into another file. Then you should be able to process both files in paired-end (-p) and single-end mode (-s), respectively.

I personally would probably stick with a) since it is the most straight forward method and you don't exclude data that is actually there. Let me know if I can be of any help with whatever you choose to do.

[Dipro]

Hi fkrueger,

Thanks very much for your reply. I'll work on a) as adviced.

Much appreciated.