Multi-mappers about methylpy HOT 3 CLOSED

The term "multi-mapped reads" is a bit confusing. It can refer to two not mutually exclusive cases:

1. reads mapped to multiple locations in the C-to-T (forward reference) genome or G-to-A (reverse reference) genome.
2. reads mapped to both C-to-T genome and G-to-A genome.

If you are interested in case 2, you can set --keep-temp-files to True to keep the results of read alignment to C-to-T genome and G-to-A genome. The files are * _forward_strand_hits.sam and * _reverse_strand_hits.sam. Please make sure there are enough hard drive space since they will usually occupy ~500g or more.

If you care about reads falling in case 1 or both case 1 and 2, in addition to setting --keep-temp-files to True, you will need to set --aligner-options` such that the aligner (e.g. bowtie2) will report reads mapped to multiple locations.

from methylpy.

You are right, the term was confusing. I am interested in case 1, when reads map to multiple locations in the genome.

I tried the --aligner-options with bowtie2, but I do not think that positively worked, because of 2 reasons:

1. I think after the bam alignment files are made, the methylpy program filters through the reads to only analyze unique-mapped reads. This is done by using samtools after the alignment files are made. So, correct me if I am wrong please, changing the aligner option might not have the desired effect.

2. When I remove one of the regions in my genome (delete the sequence), that I know causes the loss of many reads due to multi-loci mapping, as it exactly matches another region of the genome, I actually get more reads in the bam file, compared to when I use the "--all" option in the --aligner-options. This demonstrates that "--all" option is not working as expected.

That is why I was wondering, if there is any way to not exclude the multi-loci mapping reads in the methylpy pipeline, not just at the aligner stage but during methylation calling too.

from methylpy.

You are right that methylpy filter out multi-map reads. So, the multi-map reads (either case 1 or case 2) won't be available in the final bam file. However, with --keep-temp-files set to True and the bowtie2 options that allow reporting multimap reads, you should be able to find multimap reads in * _forward_strand_hits.sam and * _reverse_strand_hits.sam files.

from methylpy.