QC Based on Nucleosome Signal?

[salwanbutrus]

Hello,

Thank you for developing this invaluable tool! We noticed that one of the R packages for analyzing ATAC-seq data known as Signac computes a NucleosomeSignal score https://stuartlab.org/signac/reference/nucleosomesignal for each cell, filtering cells with a high number of mono/di-nucleosome fragments relative to nucleosome-free fragments. Is there an analogous feature in pycisTopic?

A related question: can you explain the macs2 arguments --shift 73 --extsize 146 in the context of ATAC-seq data? Are fragments >146bp not being considered in the peak calling?

Thanks for your time in advance,
Salwan

[SeppeDeWinter]

Hi @salwanbutrus

I transferred your question from the issues tab to the discussions one.
At a sample level pycisTopic has functions to generate plots on the distribution of fragments at different lengths, however we don't quantify the fraction of mono- and di-nucleosome fragments per cell nor do we use it to filter cell barcodes.

Related to your second question, please see the relevant documentation on the MACS github page:

--extsize

While --nomodel is set, MACS3 uses this parameter to extend reads in 5'->3' direction to fix-sized fragments. For example, if the size of the binding region for your transcription factor is 200 bp, and you want to bypass the model building by MACS3, this parameter can be set as 200. This option is only valid when --nomodel is set or when MACS3 fails to build model and --fix-bimodal is on.

--shift

(...) You can set an arbitrary shift in bp here to adjust the alignment positions of reads in the whole library. Please use discretion while setting it other than the default value (0). When --nomodel is set, MACS3 will use this value to move cutting ends (5') then apply --extsize from 5' to 3' direction to extend them to fragments. When this value is negative, the cutting ends (5') will be moved toward 3'->5' direction, otherwise 5'->3' direction. Recommended to keep it as default 0 for ChIP-Seq datasets, or -1 * half of EXTSIZE together with --extsize option for detecting enriched cutting loci such as certain DNAseI-Seq datasets. Note, you can't set values other than 0 if the format is BAMPE or BEDPE for paired-end data. The default is 0.

Here are some examples for combining --shift and --extsize:

To find enriched cutting sites such as some DNAse-Seq datasets. In this case, all 5' ends of sequenced reads should be extended in both directions to smooth the pileup signals. If the wanted smoothing window is 200bps, then use --nomodel --shift   -100 --extsize 200.

For certain nucleosome-seq data, we need to pile up the centers of nucleosomes using a half-nucleosome size for wavelet analysis (e.g. NPS algorithm). Since the DNA wrapped on nucleosome is about 147bps, this option can be used: --nomodel --shift 37 --extsize   73.

In short these parameters are only used to move and resize fragments before calling peaks, fragments > 146 bp are still considered in the analysis.

I hope this answers your questions?

All the best,

Seppe

[ghuls]

@salwanbutrus nucleosome_signal statistics are now calculated when running pycistopic qc when you install the latest git version.