PyMaSC

Python implementation to calc mappability-sensitive cross-correlation for fragment length estimation and quality control for ChIP-Seq.

🐾 Visit PyMaSC web site for more information and to get human genome mappability tracks.

📃 If you use this software in your work, please cite the following paper.

Anzawa, Hayato, Hitoshi Yamagata, and Kengo Kinoshita. "Theoretical characterisation of strand cross-correlation in ChIP-seq." BMC bioinformatics 21.1 (2020): 417. https://doi.org/10.1186/s12859-020-03729-6

1. Introduction
2. Install
3. Usage
   • pymasc command
   • pymasc-precalc command
   • pymasc-plot command
4. Computation details
5. References

Introduction

To estimate mean fragment length for single-ended sequencing data, cross-correlation between positive- and negative-strand reads is commonly used. One of the problems with this approach is phantom peak, which is the occasionally observed peak corresponding to the read length. In the ChIP-Seq guidelines by ENCODE consortia, cross-correlation at fragment length and read length are used for quality control metrics. Additionally, library length itself is one of the important parameters for analysises. However, estimating correct flagment length is not a easy task because of phantom peak occurrence. P Ramachandran et al. proposed MaSC, mappability-sensitive cross-correlation to remove the bias caused from ununiformity of mappability throughout the genome. This method provides cross-correlation landscape without phantom peak and much accurate mean fragment length estimation. PyMaSC is a tool implemented by python and cython to visualize (mappability-sensitive) cross-correlation and estimate ChIP-Seq quality metrics and mean fragment length with MaSC algorithm.

Install

Python version 3.8 or higher is recommended (Python 3.8-3.14 officially supported). C compiler needs to build C sources (recommend GCC).

Install using pip

Install PyMaSC from PyPI

$ pip install pymasc

If cython is installed, PyMaSC will be built with Cython native compiled sources.

Usage

After installation, PyMaSC provides pymasc, pymasc-precalc and pymasc-plot command. Note that pymasc-precalc is not essential to calculate mappability-sensitive cross-correlation.

pymasc command

pymasc [-h]
       [-v {DEBUG,INFO,WARNING,ERROR,CRITICAL}]
       [--disable-progress]
       [--color {TRUE,FALSE}]
       [--version]
       [-p PROCESS]
       [--successive]
       [--skip-ncc]
       [--skip-plots]
       [-r READ_LENGTH]
       [--readlen-estimator {MEAN,MEDIAN,MODE,MIN,MAX}]
       [-l LIBRARY_LENGTH]
       [-m REGION_FILE]
       [--mappability-stats MAPPABILITY_STATS]
       [-q MAPQ]
       [-i CHROM [CHROM ...]]
       [-e CHROM [CHROM ...]]
       [-d MAX_SHIFT]
       [--chi2-pval CHI2_PVAL]
       [-w SMOOTH_WINDOW]
       [--mask-size MASK_SIZE]
       [--bg-avr-width BG_AVR_WIDTH]
       [-n NAME [NAME ...]]
       [-o OUTDIR]
       reads [reads ...]

Usage example

Calculate only naïve cross-correlation for ENCFF000VPI.bam with max shift size = 1000. Output files are ENCFF000VPI_stats.tab, ENCFF000VPI_nreads.tab, ENCFF000VPI_cc.tab and ENCFF000VPI.pdf.

$ pymasc -d 1000 ENCFF000VPI.bam

Calculate both naïve and mappability-sensitive cross-correlation by 4 worker processes. Output ENCFF000VPI_mscc.tab Additionally.

$ pymasc -p 4 -d 1000 -m wgEncodeCrgMapabilityAlign36mer.bigWig ENCFF000VPI.bam

Main input file

SAM and BAM file format are acceptable.

• Input alignment file must be sorted.
• Additionally, for parallel processing, input file must be BAM format and indexed.
• If multiple files specified, PyMaSC processes each file with common parameters.
• Unmapped or duplicated reads will be discarded.
• If input file contains paired-end reads, the last (second) segment read will be discarded.

Output files

Additionaly, PyMaSC generates a JSON file as cache for mappability analyses. See the --mappability-stats option for details.

General options

Processing settings

Input alignment file settings

Input mappability file settings

Input file filtering arguments

Analysis Parameters

Output options

pymasc-precalc command

pymasc-precalc [-h]
               [-p PROCESS]
               [-v {DEBUG,INFO,WARNING,ERROR,CRITICAL}]
               [--disable-progress]
               [--color {TRUE,FALSE}]
               [--version]
               [-m REGION_FILE]
               [--mappability-stats MAPPABILITY_STATS]
               [-d MAX_SHIFT]
               [-r MAX_READLEN]

Usage example

Calculate total length of doubly mappable regions. wgEncodeCrgMapabilityAlign36mer_mappability.json will be write.

$ pymasc -p 4 -r 50 -d 1000 -m wgEncodeCrgMapabilityAlign36mer.bigWig

Options

Almost same as pymasc command. Note that actual max shift size is,

• 0 to read_length (if max_shift < read_len * 2)
• 0 to max_shift - read_len + 1 (if max_shift => read_len * 2)

pymasc-plot command

pymasc-plot [-h]
            [-v {DEBUG,INFO,WARNING,ERROR,CRITICAL}]
            [--disable-progress]
            [--color {TRUE,FALSE}]
            [--version]
            [--stats STATS]
            [--cc CC]
            [--masc MASC]
            [--nreads NREADS]
            [-s SIZES]
            [-m MAPPABILITY_STATS]
            [-i CHROM [CHROM ...]]
            [-e CHROM [CHROM ...]]
            [--chi2-pval CHI2_PVAL]
            [-w SMOOTH_WINDOW]
            [--mask-size MASK_SIZE]
            [--bg-avr-width BG_AVR_WIDTH]
            [-l LIBRARY_LENGTH]
            [-n NAME]
            [-o OUTDIR]
            [-f [{all,stats,cc,mscc} [{all,stats,cc,mscc} ...]]]
            [statfile]

Usage example

(Re)plot figures from pymasc outputs output/ENCFF000VPI_* with the specified smoothing window size and library length. Note that you need to specify a tab delimited file or a SAM/BAM format file to obtain chromosome lengths, and/or, specify a mappability stats (JSON) file which was generated for a mappability BigWig file by PyMaSC to obtain mappable region lengths.

$ pymasc-plot -w 50 -l 250 -s ENCFF000VPI.bam output/ENCFF000VPI

Input argument

Specify a prefix to pymasc output files. For example, set output/ENCFF000VPI to plot figures from output/ENCFF000VPI_stats.tab, output/ENCFF000VPI_nreads.tab and output/ENCFF000VPI_cc.tab (and/or output/ENCFF000VPI_mscc.tab). *_stats.tab and either or both of *_cc.tab and *_mscc.tab must be exist. To specify these files individually, use --stats, --nreads, --cc and --masc options.

Computation details

BitArray and successive implementation

PyMaSC provides two algorithms for calculation.

BitArray

BitArray approach is based on allocated binary arrays with length of references and computation time mostly depends on the size of reference genome and the maximum shift size. Typically (hg19, chr1), PyMaSC consumes about 250MB RAM per worker.

Successive implementation

Successive implementation is based on buffer with length of maximum shift size and calculate overwrapped bases successively while reading reads and mappable regions. Computation time mostly depends on number of reads and tracks in input files and BitArray implementation is faster in most cases. Successive approach surpasses in processing small input data, quite low memory efficiency and robustness for shift size.

References

• PyMaSC paper
  • Anzawa, Hayato, Hitoshi Yamagata, and Kengo Kinoshita. "Theoretical characterisation of strand cross-correlation in ChIP-seq." BMC bioinformatics 21.1 (2020): 417.
• Original paper on the MaSC algorithm
  • Ramachandran, Parameswaran, et al. "MaSC: mappability-sensitive cross-correlation for estimating mean fragment length of single-end short-read sequencing data." Bioinformatics 29.4 (2013): 444-450.