Hi everyone, first of all Iam a biologist so sorry for the (potentially) very basic question. 

We performed ATAC-seq (treatment vs. control in duplicates) and processed the data through the Encode pipeline. The QC (especially the TSS enrichment) shows for both control samples a excellent result whereas both treatment samples are borderline ( but not so bad to trash it). Also the visual inspection of the peaks looks okay (very similar peak pattern between control and treatment). However we noticed that the peaks in the treatment samples are overall smaller compare to the control (bigwig TPM normalized) so we concluded that we have a potential enrichment bias. Since we do not knock down any global factor it is most likely a technical artefact. We continued to identify differential accessible regions using the csaw package ("trended" normalization). The resulting regions are highly enriched for binding motives of biological meaningful transcription factors (identified with HOMER) and also close to relevant genes. 

So far so good... however now I have problems to visualise those regions. Since the TPM normalization of the bigwig files does not account for the global enrichment bias (correct ?) the plotted heatmaps of the differential accessible regions do not reflect the result from csaw (control has always higher signal compare to treatment). Please correct me when Iam wrong but in this case we could perform instead of the TPM normalization a quantile normalization right ? Is there a easy way to do this either with bam or bigwig files ? Otherwise I could get a count matrix from deeptools and feed it into the CQN package. However the cqn function requires a covariate and to be honest Iam not sure what this is in my data (sequencing depths ?). After the normalization I could convert the resulting file back into a bigwig file and us it for plotting. 

Thanks for any suggestions.  

What you describe seems to be a difference in signal-to-noise ratio which is not uncommon. <br>
This is where more elaborate normalization methods such as TMM from `edgeR` or  RLE from `DESeq2` come into play. <br>

See the following links on why these methods are superior. The videos talk about RNA-seq but simply consider genes as "regions" in any assay such as ATAC-seq or ChIP-seq: <br>

 - [StatQuest video DESeq2][1] 
 - [StatQuest video edgeR][2]

Also, this is where more simple methods such as TPM/RPKM/FPKM which only correct for total library read count <br> (=sequencing depth) typically make limited sense as they fail to correct for the systematic differences in library composition due to the differences in either the peak landscape or signal/noise ratio. I cannot comment on Quantile Normalization as I have no experience with it but I had good results with both TMM and RLE. Therefore my recommendation to routinely use something like `edgeR` to calculate scaling factors which you then use to divide the column 4 of your bedGraph/bigwig by.

**BedGraph/Bigwig creation**:

Among other tools this could be done by `bamCoverage` from deeptools or `genomecov` from bedtools. Both tools start from BAM files and offer plenty of options for customization. 

**Scaling Factors:**

Here is a code suggestion in `R` to calculate scaling factors to scale the bedGraph/bigwig files. It assumes you have already created a count matrix with rows = regions and columns = samples, e.g. based on a list of peaks called from the underlying experiment. This could be done with tools such as `macs2` or `Genrich` followed by `featureCounts`. It is assumed that the matrix of **raw** counts is now loaded in `R` as object `raw.counts`:

    ## edgeR:: calcNormFactors
    NormFactor <- calcNormFactors(object = raw.counts, method = "TMM")

    ## if you prefer to use the DESeq2 strategy use method="RLE" instead
    
    ## raw library size:
    LibSize <- colSums(raw.counts)
    
    ## calculate size factors:
    SizeFactors <- NormFactor * LibSize / 1000000
     
    ## Reciprocal, please read section below:   
    SizeFactors.Reciprocal <- 1/SizeFactors

`SizeFactors` will now contain a factor for every sample which can be used to divide the 4th colun of a bedGraph/bigwig by.
Both the aforementioned tools (`bamCoverage` and `genomecov`) have options though to directly scale the output by a factor (`--scaleFactor` or `--scale` respectively).
**!! Note though** that these options will multiply rather than divide the 4th column with the factor, therefore you would need to provide the reciprocal as mentioned in the code chunk above.

----------

**Illustration:**

Here is a picture I made a while ago addressing exactly this problem. It shows two biological replicates with different signal-to-noise ratios (replicate two has worse SN therefore overall lower peak counts as many reads fall into non-peak regions) for a region which I know is not differential. Using CPM (so per-million scaling based on only the read counts) shown in the top panel fails to correct for the systematic difference while TMM on the counts in peak regions (bottom panel) manages to correct for it. The middle panel is a bin-based normalization with edgeR suggested in the `csaw` package which might in certain situations be beneficial if library composition is very different between samples (check csaw vignette for details if you want), but lets not further discuss this here as it does not add to the topic right now.

![enter image description here][3]

**Bigwigs in R:**

If you want to visualize bigwigs in R one might want to check this handy tool that was recently posted here which seems to be easy to use yet powerful and visually appealing:

https://www.biostars.org/p/475853/


  [1]: https://www.youtube.com/watch?v=UFB993xufUU
  [2]: https://www.youtube.com/watch?v=Wdt6jdi-NQo
  [3]: https://i.imgur.com/6Yzvboo.png