10
$\begingroup$

We have heard in the group that it is important to keep track of and to filter artifact regions when analysing data from functional genomics experiments, especially ChIP-seq.

Here, we have seen pipelines that remove the ENCODE tracks i) before cross-correlations QCs, ii) after cross-correlations QC but before peak calling and iii) after peak calling.

We have noticed that removal of the tracks does not affect significantly cross-correlation and peak-independent QCs. However, we are not sure whether peak calling should be done on the filtered tracks or not?

$\endgroup$

2 Answers 2

10
$\begingroup$

Aside: Cross-correlation is largely meaningless, regardless of what some of the ENCODE folks might argue. When we process our DEEP samples we don't even look at that value.

Regardless, if you're using SPP/phantomPeakQual for cross-correlation then note that it already removes the highest peaks from your dataset before computing the cross-correlation (in fact, it can remove most of the actual peaks too, which makes one further wonder what it's actually telling you). I don't know that this is actually documented anywhere, it's something I noticed when going through the code while pondering whether to implement it in deepTools. But at least it's ignoring these regions already :)

In general, it's most convenient to just remove peaks overlapping blacklisted regions. In an ideal world you'd filter out the blacklisted reads before peak calling, but (1) this is really inconvenient (more time and disk required) and (2) I've never seen an appreciable gain in peak calling performance. In theory at least you should be losing sensitivity right around blacklisted regions if you don't remove reads overlapping blacklisted regions, but you have to ask yourself whether you want to trust such peaks anyway. For other QC steps, at least with deepTools we provide a parameter with every tool to specify a BED file of blacklisted regions to skip.

As an aside, there are many fewer blacklisted regions in more recent genome builds (GRCh38 and GRCm38 at least), so this is less of an issue in general with them.

$\endgroup$
0
$\begingroup$

I filtered my bams out of curiosity. I had around 500 peaks before filtering. When i called on the filtered bams, i lost 100 (from the blacklist) and gained 38. The peak call was using macs2. Looking at the gained peaks, they are legit and a binding partner of my TF had peaks there. Note however i ran out of ram for bam > 8GB using bedtools intersect to filter. And if you do it, sort the bam and make sure the chromosome names are ordered using the same tool between blacklist bed and your bams.

$\endgroup$
1
  • $\begingroup$ It is much more efficient to use samtools view. It accepts a BED file and only keeps reads overlapping this. So you take a BED file of the entire genome, make the complement of this with your blacklist and what is left is the regions to keep. Provide that to samtools view, this is fast with almost no memory footprint. $\endgroup$
    – ATpoint
    May 22, 2023 at 11:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.