Strategy for noisy ATAC-seq: Fragment size selection for pioneer factor DA analysis

Hello everyone.

I have standard bulk ATAC-seq experiment with KD treatment, the exp design is 4 replicates per one condition. I aligned, and filtered the bams from mito reads + dup + non-standard chroms + black listed regions and did a TN5 shift (9 bp) and my primary QC looks similar to the post: ATAC-SEQ FRAGMENT SIZE I called peaks using macs2 and have ~20 k peaks per replicate and ~ 90k peaks when calling peaks on the pooled bam from replicates. My tracks are noisy and it is hard to distinguish peaks from background, however, because there is nucleosomal pattern of the noise, I fractioned my bam into different insert sizes to get rid of di- and tri- etc nuclosomal fragments. My TF of which I am investigating is a pioneer factor (Ascl1) that is known to bind mono-neocleosomes in addition to NFR, so my thoughts were to keep such inserts as they might relevant and not to lose much seq depth. Would it be biologically sound to focus only on the fraction from ~ 0-180 bp to mitigate the noise effects and do a DA on that subset? and if I lose considerable part of coverage ~30-40% of filtered total 40M fragments per replicate, would that be losing critical power in downstream analysis (if I start see disperse genome-wide coverage where some regions of background lack reads) with tools like deseq2 on counts over consensus peak list or metaplots/profile aggregates. NPC cell-line track is a control from matching cells from literature deposited datasets

p.s. I tagged deeptools as many analysis steps using profiles and heatmaps for genomic metaplots are to be done by deeptools.

deeptools deseq2 atac-seq bam

• 1,051 views • link

written 7 months ago by Abeer H •

1 answer

Your data are indeed noisy. Let me say that I do not agree with the answer in the linked thread on insert sizes stating that shown plots indicate good quality. You should at least see a banding pattern beyond the first two major peaks, something like this below. Below is two random examples from our own data, which both I would call "good", based on QC, and how things looked towards callable peaks, FRiPs and on the IGV.

enter image description here

Unlike the ones in the link here you see several peaks at expected locations correcponding to nucleosome-free, 1n, 2n- and so on histones.

As for your question, I don't think that size selection denoises data. If the library is not good then this is continuous, meaning subsetting to certain sizes will then still be suboptimal quality for this chunk, as I see no reason why certain insert sizes are more likely to be signal than noise.

Just do the analysis and see what you get. If inconclusive, or trends but not significant, consider repeating the experiment. There is no in silico magic to magically rescue a poor library.

As for the peaks, pooling bam files (for me) always vastly exaggerated the number of callable peaks as compared to what I see in the IGV. I had some success with Genrich which can assess reproducibility of peaks when calling in replicate mode. Alternatively, use call per sample, and only keep peaks called in at least 75% of samples per group, then merge. Lots of ways, but generally I would do some sort of assessment that filters for reproducibility across replicates, simply to remove nonsense calls which macs can produce a lot of and which then noises the DA analysis.

• 0 views • link

written 7 months ago by ATpoint

Hello ATpoint , thank you for your reply.

I will definitely give Genrich a try to see if replicate-aware caller (or the peaks in 3/4 reps strategy) will give a more confident list of peaks as having a consensus peak list is a critical must for DA.

So regarding the bam fractions, eventhough the primary QC from my samples shows nuclosomal pattern (image attached from ataqv), I have good portion of reads distributed in mono- di- tri nuclosomal regions, that's why the peaks are broad and I am losing coverage on those cutting events. Please note there are wavey descending signal (in the fraction containing insert sizes >180 bp) around the original peaks, this fraction is contributing significantly to the noise around peaks, affecting peak calling (absent/present, peak borders). If you compare tracks (counting after the refseq track), track 1 (all fractions included) with track 5 (NFR+mono), with track 4 (inserts > 180 bp) which is merely extended genomic diffused signal. This reduced my TSS enrichment metrics, and those are inserts useful for nuclosomal analysis but not for DA. As I understand they increase the complexity of my libraries on the expense of losing coverage on the basic events that matter for detecting changes. The treatment changes are very limited due to short treatment, so I don't want those changes to be lost due to noise in very close regions to real peak borders where TF binding events occur.

So when I separate them depending on insert size, the NFR fraction peaks seem similar to the clean track (on the top of Refseq track) from literature data, though ~50% less signal/coverage. I know that many analysis proceed with NFR fraction for more downstream analysis (inserts below ~100/140 bp) but in my case the TF can bind mono-nuclosomal regions as well. So would DA on that fraction (inserts of <180 bp / NFR+mono) be biologically defendable/sound?

enter image description here

• 0 views • link

written 7 months ago by Abeer H •