TSS heatmap inconsistency between ChIPseeker and nf-core ATAC-seq pipeline

[hyjforesight]

Hello ATAC-seq pipeline,
Previously, we did the ATAC-seq analysis manually by trimming adapters (Adapterremoval), alignment (Bowtie2), QC (remove blacklist, duplicates, unmapped reads, mito genes…), peak calling (MACS2), and peak annotation (ChIPseeker).

Luckily, we found this amazing and easy-to-use ATAC-seq pipeline on nf-core. We tried this pipeline and it works super great.

However, we are confused with 2 issues:

1. The TSS heatmaps are inconsistent between our ChIPseeker results and nf-core ATAC-seq pipeline. In our hands, we can see heatmaps with signals, but the pipeline shows no signal at all. Which heatmap is right?
   

2. We found that the pipeline only calls broad peak, but ENCODE ATAC-seq pipeline calls narrow peak. Why does nf-core pipeline only call broad peak?
   

Thank you!
Best,
Yuanjian Huang

[bjlang]

With respect to your 2nd issue:
The pipeline defaults to calling broad peaks. To activate narrowPeak mode, use --narrow_peak (see https://nf-co.re/atacseq/1.2.2/parameters#peak-calling-options)

[lnblum]

I've experienced the same issue. I think it has to do with the genes.tss.bed file. I ran deeptools manually and swapped it out for refTSS (attached below - refTSS_v3.1_human_coordinate.hg38.bed) and got a result that makes sense. I have not been able to diagnose the problem with genes.tss.bed but maybe deeptools is looking at the wrong columns in the bed file?
50K-15-S-40-CEB_R1.plotHeatmap.pdf
50K-15-S-40-CEB_R1.mLb.clN.plotHeatmap.pdf
refTSS_v3.1_human_coordinate.hg38.bed.txt

[robsyme]

There are three issues, I think:

1. We are supplying the gene bed rather than a TSS bed, which is the primary cause of the poor results.
2. We are not using the --missingDataAsZero parameter to computeMatrix scale-regions
3. Even if we switched the bed file to the TSS bed emitted in PREPARE_GENOME.out.tss_bed, that would not work as the bin width in computeMatrix scale-regions is 10bp and so the 1bp wide regions in PREPARE_GENOME.out.tss_bed would all be filtered out.

I'd suggest that we use the PREPARE_GENOME.out.tss_bed file, but expand all regions to at least 10 bp. This would involve either some awk magic in the DEEPTOOLS_COMPUTEMATRIX process, or an upstream bedtools process that then feeds into DEEPTOOLS_COMPUTEMATRIX.

If we expand the PREPARE_GENOME.out.tss_bed so that it meets the computeMatrix bin width minimum, we get sensible results that look like the heatmap attached.

[drpatelh]

Fixed in #228

Be great if you can confirm by trying out the code on the dev branch. We are still fixing some stuff there but the plots should be generated as expected.

nextflow pull nf-core/atacseq -r dev
nextflow run nf-core/atacseq <YOUR_OPTIONS> -r dev

[lnblum]

Can confirm. I ran the dev version on some data last week and the heatmaps were generated as expected.

[hyjforesight]

Hello @drpatelh
Thank you for fixing the issue!

Could you please let me know how to adjust the scale bars to the same scale between different samples? Now they're 1.6, 1.72, 1.2 and 1.2.

Thank you in advance!