Aligning ChIP-Seq reads to repeats for downstream peak analysis

Question

This is a brief question regarding the above. I have previously used bowtie to map reads from paired-end ChIP-Seq sequencing, then used the positions for peak-calling. However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

For peak-calling purposes, however, I need to know reasonably confidently where precisely those repeats are located, so that MACS2 treats the coverage pileup sensibly. So my question is this: what's the best way of mapping these (short) paired-end reads from a ChIP-Seq experiment so I can do this?

EDIT: I want to be able to run a peak-caller after aligning the data, so I want a good approximation to the coverage at repeating elements after mapping for instance - I don't want to throw away all reads which multimap, as may mean that some reads could lose their mates, and so I won't be able to calculate coverage or call peaks using the data.

I want an aligner to keep both ends of the paired-end pair within about 180bp of each other (as that is the expected length of the DNA fragment that is sequenced). The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced. With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I have tried bowtie (and bowtie2) for this, which is where I understood that the dataset maps to repeats.

Answer 1

This is very good question, but as far as I know there really isn't a good answer. I will attempt to provide some comments and tips.

However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

If this dataset is also PE and the alignment is outputting either uniquely mapped reads or a single (random) position for multimappers, it is pretty much as good as it gets. A lot of repeat regions are not exact copies, so using PE there is a good change that one read of the pair mapping to a single position is helping with mapping of the other read.

I need to know reasonably confidently where precisely those repeats are located

It is probably a matter of phrasing, but what one needs is to know the precise genomic location of the reads. The location of the repeat elements can be obtained from the Repeat Masker track of UCSC, but you probably know this already.

With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I think bowtie already does this. Did you map the read pairs together or separately?

Note that:

• There will always be a number of multimapping reads (fragments) that will be discarded. The information of how many are being discarded for your samples would be helpful.
• Mapping can be improved using both reads and tunning bowtie settings (see bellow)

Finding a valid paired-end alignment where both mates align to repetitive regions of the reference can be very time-consuming. By default, Bowtie avoids much of this cost by imposing a limit on the number of "tries" it makes to match an alignment for one mate with a nearby alignment for the other. The default limit is 100. This causes bowtie to miss some valid paired-end alignments where both mates lie in repetitive regions, but the user may use the --pairtries or -y options to increase Bowtie's sensitivity as desired. source

The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced.

As indicated above, the aligner will take care of that sort of. It will gave the mapping information of each read pair, which is then used by MACS.

Other tools

I tested ShortStack which uses the coverage of uniquely mapped reads to place multimapping reads for ChIP-seq, and the results were mixed. That said, it was SE data and the code has been through some changes, so it be worth a try.

When requested, the authors of RepEnrich provided me with the scripts for an utility called RepConsensus which maps reads to a consensus reference of repeat regions. The caveats are:

• mappings are better for elements whose sequences have not diverged a lot yet
• it will tell the position of the protein in the repeat sequence, but no information about the surrounding genomic regions because there isn't any.

---

Whatever you do, do not use all alignments, that is double counting of reads, and think very carefully before using custom methods. A blog post detailing some of the potential issues can be found here (I wrote it).