Quality control metrics for ATAC-seq data including fragment size distribution, TSS enrichment, FRiP, and library complexity. Use when assessing ATAC-seq library quality before or after peak calling to identify problematic samples.
Before using code patterns, verify installed versions match. If versions differ:
Python: pip show <package> then help(module.function) to check signatures
R: packageVersion('<pkg>') then ?function_name to verify parameters
CLI: <tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
ATAC-seq Quality Control
"Check the quality of my ATAC-seq library" → Evaluate fragment size distribution (nucleosome periodicity), TSS enrichment, FRiP, and library complexity to assess chromatin accessibility experiment quality.
Goal: Assess ATAC-seq library quality by visualizing the characteristic nucleosome periodicity in fragment sizes.
Approach: Extract insert sizes from the BAM file using Picard or samtools, producing a distribution that should show NFR (<100 bp) and mono-nucleosome (~200 bp) peaks.
Goal: Measure library complexity to detect over-amplification or low-diversity libraries.
Approach: Calculate NRF (unique/total reads), PBC1 (1-read locations / all locations), and PBC2 (1-read / 2-read locations) using Picard or custom counting.
# Using Picard EstimateLibraryComplexity
java -jar picard.jar EstimateLibraryComplexity \
I=sample.bam \
O=complexity.txt
# Or calculate from BAM
# NRF = unique reads / total reads
# PBC1 = locations with exactly 1 read / locations with >= 1 read
# PBC2 = locations with exactly 1 read / locations with exactly 2 reads
import pysam
def calculate_complexity(bam_file):
'''Calculate library complexity metrics.'''
bam = pysam.AlignmentFile(bam_file, 'rb')
positions = {}
total = 0
for read in bam.fetch():
if read.is_unmapped or read.is_secondary:
continue
total += 1
pos = (read.reference_name, read.reference_start)
positions[pos] = positions.get(pos, 0) + 1
distinct = len(positions)
m1 = sum(1 for v in positions.values() if v == 1)
m2 = sum(1 for v in positions.values() if v == 2)
nrf = distinct / total if total > 0 else 0
pbc1 = m1 / distinct if distinct > 0 else 0
pbc2 = m1 / m2 if m2 > 0 else 0
return {'NRF': nrf, 'PBC1': pbc1, 'PBC2': pbc2}
Goal: Generate a single QC summary combining all major ATAC-seq quality metrics.
Approach: Run samtools and bedtools commands to collect total reads, mapping rate, mitochondrial fraction, FRiP, and peak count, then write a consolidated report.
import subprocess
import pandas as pd
def atac_qc_report(bam_file, peaks_file, output_prefix):
'''Generate comprehensive ATAC-seq QC report.'''
metrics = {}
# Total reads
result = subprocess.check_output(f'samtools view -c -F 4 {bam_file}', shell=True)
metrics['total_reads'] = int(result.strip())
# Mapped reads
result = subprocess.check_output(f'samtools view -c -F 4 -F 256 {bam_file}', shell=True)
metrics['mapped_reads'] = int(result.strip())
# Mitochondrial reads
result = subprocess.check_output(f'samtools view -c {bam_file} chrM', shell=True)
metrics['mt_reads'] = int(result.strip())
metrics['mt_fraction'] = metrics['mt_reads'] / metrics['total_reads']
# Reads in peaks (FRiP)
result = subprocess.check_output(
f'bedtools intersect -a {bam_file} -b {peaks_file} -u | samtools view -c', shell=True)
metrics['reads_in_peaks'] = int(result.strip())
metrics['frip'] = metrics['reads_in_peaks'] / metrics['total_reads']
# Peak count
result = subprocess.check_output(f'wc -l < {peaks_file}', shell=True)
metrics['peak_count'] = int(result.strip())
# Write report
with open(f'{output_prefix}_qc.txt', 'w') as f:
for k, v in metrics.items():
if isinstance(v, float):
f.write(f'{k}: {v:.4f}\n')
else:
f.write(f'{k}: {v}\n')
return metrics
