fragmentomics run¶

1. About¶

The fragmentomics executable is composed of several inter-related sub commands. Please see fragmentomics -h for all available options.

This part of the documentation describes options and concepts for fragmentomics run sub command in more detail. With minimal configuration, the run sub command enables you to start running fragmentomics pipeline.

Setting up the fragmentomics pipeline is fast and easy! In its most basic form, fragmentomics run only has three required arguments.

2. Synopsis¶

$ fragmentomics run [--help] [--overwrite-pipeline-template] \
      [--dry-run] [--job-name JOB_NAME] [--mode {{slurm,local}}] \
      [--sif-cache SIF_CACHE] [--singularity-cache SINGULARITY_CACHE] \
      [--silent] [--threads THREADS] [--tmp-dir TMP_DIR] \
      [--fragment-minimum FRAGMENT_MINIMUM] \
      [--fragment-maximum FRAGMENT_MAXIMUM] \
      [--left-tss-flank LEFT_TSS_FLANK] \
      [--right-tss-flank RIGHT_TSS_FLANK] \
      [--split-interval SPLIT_INTERVAL] \
      [--interval INTERVAL] \
      [--bin-size BIN_SIZE] \
      [--mapscore MAPSCORE] [--baseqscore BASEQSCORE] \
      [--cda-features CDA_FEATURES] \
      [--cda-regions CDA_REGIONS] \
      [--cda-cna-bin-size CDA_CNA_BIN_SIZE] \
      -g {{hg38,hg19}} \
      --input INPUT [INPUT ...] \
      --output OUTPUT

The synopsis for each command shows its arguments and their usage. Optional arguments are shown in square brackets.

A user must provide a list of Illumina FastQ or BAM files (globbing is supported) to analyze via --input argument, an output directory to store results via --output argument, and a reference genome build via the --genome argument.

Use you can always use the -h option for information on a specific command.

2.1 Required arguments¶

Each of the following arguments are required. Failure to provide a required argument will result in a non-zero exit-code.

--input INPUT [INPUT ...]

Input Paired-end Illumina FastQ or BAM file(s).
type: file(s)

One or more paired-end Illumina FastQ files or one or more BAM files can be provided. Only a single input type is allowed per run (i.e all inputs must either be all FastQ or BAM files). Input FastQ and BAM files cannot be mixed in a single run (i.e --input *.bam *.fastq.gz). The input type is auto-detected from the files you provide. It is worth noting that single-end FastQ files are not supported. From the command-line, each input file should be seperated by a space. Globbing is supported! This makes selecting input files easy.

Example: --input .tests/*.bam
Example: --input .tests/*.R1.fastq.gz .tests/*.R2.fastq.gz

How each input type is processed. The pipeline converges both input types on a single coordinate-sorted, indexed BAM per sample (bams/{sample}.sorted.bam) that feeds every downstream fragmentomics analysis: - FastQ input runs the full pipeline. Paired-end reads are quality filtered, aligned to the selected --genome build with bwa-mem2, and reduced to properly-paired, primary, mapped reads with duplicates marked (but kept) before the analysis steps run. See FastQ alignment. - BAM input runs the full pipeline minus alignment. Provided alignments are staged (quality filtered, coordinate-sorted and indexed), stripped of singletons and of the orphaned mates that per-read filtering leaves behind, checked against the selected genome build's sequence dictionary, and subset to the contigs the two share. See BAM normalization and Reference contig filter.

--output OUTPUT

Path to an output directory.
type: path

This location is where the pipeline will create all of its output files, also known as the pipeline's working directory. If the provided output directory does not exist, it will be created automatically.

Example: --output pipeline_output

--genome {hg38,hg19,GENOME_CONFIG.json}

Reference genome build.
type: string
default: hg38

Selects the set of reference files the pipeline uses to characterize cfDNA fragmentation features. Choosing a genome build determines which chromosome sizes, 2bit reference sequence, genomic interval and TSS files, and any blacklist or gap files are used throughout the analysis. For FastQ input, it additionally selects the bwa-mem2 reference index that reads are aligned against. For BAM input, it selects the sequence dictionary the input alignments are validated against. This argument is required for both FastQ and BAM input (BAM inputs still need the build to select the analysis references, and are assumed to already be aligned to it).

Two builds are bundled with the pipeline and named by alias: hg38 or hg19.

A reference set of your own can be used instead, by giving the path to a JSON file describing it rather than an alias. The file holds a single build, in the shape of one entry of the pipeline's config/genome.json:

{
    "chrom_sizes":   "/refs/mm10.chrom.sizes",
    "ref2bit":       "/refs/mm10.2bit",
    "reference_fa":  "/refs/mm10.fa",
    "dict":          "/refs/mm10.dict",
    "tss":           "/refs/tss.mm10_sorted.bed",
    "tss_interval":  "/refs/tss.mm10_interval_sorted.bed"
}

The build is named after the config file, so mm10.json labels its output mm10; add a "name" key to choose the label yourself. Only tss_interval is required — every other key gates just the analyses that read it, so a build that leaves one out skips those rather than failing the run. Reference paths become container bind points automatically. Note that cfDNAanalyzer supports no build but hg19 or hg38, so a build of your own has to declare which of the two its coordinates match, with a "cda_genome" key, before --cda-features can select anything.

Every file a build can provide, which analysis uses it, and the other accepted spellings of a config file, are documented in Reference files.

Example: --genome hg38
Example: --genome /refs/mm10.json

2.2 Analysis options¶

Each of the following arguments are optional, and do not need to be provided. What each analysis does with these values is documented per analysis in the Analyses section.

--fragment-minimum FRAGMENT_MINIMUM

Minimum fragment length.
type: int
default: 50

Minimum fragment length, in base pairs. Fragments shorter than this length are excluded from the fragmentation analyses: coverage, fragment-length bins and intervals, end motifs and interval end motifs, and cleavage profiles.

Two analyses do not use this value: wps uses a fixed 120–180 bp window (the definition of L-WPS), and delfi partitions short from long internally.

Example: --min 50

--fragment-maximum FRAGMENT_MAXIMUM

Maximum fragment length.
type: int
default: 500

Maximum fragment length, in base pairs. Fragments longer than this length are excluded from the same analyses listed under --fragment-minimum, and with the same two exceptions. Raise this if you need the dinucleosome and trinucleosome range represented in the fragment-length distribution.

Example: --max 500

--left-tss-flank LEFT_TSS_FLANK

Left TSS flank size.
type: int
default: 2000

Size, in base pairs, of the region to include to the left of each transcription start site (TSS) when computing cleavage profiles. The TSS reference file holds point positions, so a region has to be built around each one before a profile can be computed. Left and right flanks are set separately, so the window need not be symmetric.

Example: --left-tss-flank 2000

--right-tss-flank RIGHT_TSS_FLANK

Right TSS flank size.
type: int
default: 2000

Size, in base pairs, of the region to include to the right of each transcription start site (TSS) when computing cleavage profiles. Widening this extends the profile further into the gene body, where the downstream nucleosome array is most regular. Note that cleavage_profile is the most expensive step in the pipeline, and its cost scales with the total flank width.

Example: --right-tss-flank 2000

--split-interval SPLIT_INTERVAL

WPS interval size.
type: int
default: 5000

Interval size, in base pairs, used when computing and adjusting window protection scores around transcription start sites. Each TSS is expanded to an interval of this size, centered on the site, and WPS is computed across every base of it. The same value is passed to adjust-wps so the two stages agree.

Setting this to 0 also disables the frag-length-intervals analysis.

Example: --split-interval 5000

-i INTERVAL, --interval INTERVAL

Genomic window size.
type: SI-prefixed base unit
default: 1mb

Width of the fixed-size windows the reference genome is tiled into. These windows are the coordinate space that frag-length-intervals, interval-end-motifs and delfi summarize over: fragment-length statistics and end-motif frequencies are reported per window, and DELFI bins over them.

The interval BED is built on the fly from the genome build's reference FastA at the start of each run, rather than read from a pre-made file. The window size is therefore a property of the run and not of the genome build — it is recorded in config.json and in the output filename, so two runs at different sizes never overwrite or silently reuse each other's windows.

Accepts a base count with an optional SI prefix, in either case. 2MB, 2mb, 2m, 2000kb and 2000000 are all two megabases; a bare number is read as a plain base count. The size is normalized before it names the file, so every spelling of a size resolves to the same BED:

--interval 2MB   ->  intervals/hg38_2mb_intervals.bed
--interval 5kb   ->  intervals/hg38_5kb_intervals.bed
--interval 500   ->  intervals/hg38_500bp_intervals.bed

Larger windows mean fewer of them, so each one accumulates more fragments and its per-window statistics are less noisy, at the cost of spatial resolution. The default of 1mb is what DELFI's published analyses bin at.

Example: --interval 2MB

--bin-size BIN_SIZE

Fragment-length bin size.
type: int
default: 1

Bin size, in base pairs, for the fragment-length distribution histogram. At the default of 1, every distinct fragment length gets its own row, which is what preserves the ~10 bp periodicity in the distribution. This value is also embedded in the fragment-length bin output filenames, so runs at different bin sizes do not overwrite each other.

Example: --bin-size 1

--mapscore MAPSCORE, --m MAPSCORE

Minimum mapping quality.
type: int
default: 20

Minimum mapping quality (MAPQ) a read must have to be kept: a read survives when its MAPQ >= MAPSCORE. The filter is applied with samtools on both input paths — on the FastQ path immediately after alignment, and on the BAM path while the input alignments are staged. Pass 0 to keep every alignment regardless of mapping quality.

The same value is passed to every FinaleToolkit analysis as its -q threshold, so one setting governs both what is written into the analysis BAM and what the analyses read out of it. See shared conventions.

Example: --mapscore 20

--baseqscore BASEQSCORE, --b BASEQSCORE

Minimum mean base quality.
type: int
default: 20

Minimum mean base quality (Phred) a read must have to be kept: a read survives when its mean base quality is >= BASEQSCORE. On the FastQ path this is applied by fastp before the reads are aligned, so a failing pair is never aligned at all; on the BAM path it is applied to the staged alignments by samtools using the filter expression avg(qual) >= BASEQSCORE. Pass 0 to keep every read regardless of base quality.

Both paths therefore mean the same thing by this option: the mean quality over a read's bases, not a per-base cutoff.

On the FastQ path the mean is measured after adapter trimming, which fastp performs in the same pass, so low-quality adapter read-through at the 3′ end cannot push an otherwise good pair below the threshold. Adapter trimming is not configurable by this option and is always on for FastQ input.

Example: --baseqscore 20

--cda-features CDA_FEATURES

Comma-separated list of cfDNAanalyzer features to extract.
type: str
default: none

cfDNAanalyzer is an additional, and fairly expensive, feature extraction suite the pipeline can run over the same analysis BAMs. Only the features named here are extracted, and none are extracted by default — an empty selection leaves cfDNAanalyzer out of the run entirely, so no rules are defined and its container is never pulled.

Names are case-insensitive. all selects every feature, none selects none. Each feature is written to cfdnaanalyzer/features/ as one or more CSV matrices of samples by measurements.

Feature extraction is all that runs

cfDNAanalyzer's feature processing, feature selection and machine learning modules are never invoked. No classifier is fit, no cross-validation is performed, and no sample is ever assigned a class or a probability. The pipeline passes --noDA and then asserts that no inference output directory was created, failing the job if one was. What you get are feature matrices; any inference over them is left entirely to you.

Genome-wide

Name	Measures
CNA	Copy number alterations, called with the bundled ichorCNA at --cda-cna-bin-size
EM	Fragment end motif frequencies and motif diversity score (MDS)
FP	Short/long fragmentation profile in 100 kb windows

Region-specific — measured over --cda-regions

Name	Measures
NOF	Nucleosome occupancy and fuzziness
NP	Nucleosome profile (mean coverage, central coverage, amplitude) at the bundled Griffin transcription factor site lists
WPS	Windowed protection score, long and short
OCF	Orientation-aware cfDNA fragmentation
EMR	End motif frequencies and MDS, aggregated over and reported per region
FPR	Fragmentation profile per region

Transcription start sites

Name	Measures
PFE	Promoter fragmentation entropy
TSSC	Average coverage around each TSS, using --left-tss-flank and --right-tss-flank

WPS run time scales with the number of regions

WPS is scored one region at a time, so its run time grows with the size of --cda-regions. The default region set is the genome build's TSS intervals — roughly 62,000 regions — which is fine for the other region-specific features but makes WPS extremely slow. Pair it with a focused region BED.

PFE expects a deep targeted panel

cfDNAanalyzer drops any sample with less than 500x median depth over the promoters it scores, so a typical whole-genome sample will be absent from the PFE matrix. This is not an error; the matrix is written with only its header.

Most features are paired-end only

EM, FP, NP, OCF, EMR, FPR and PFE are paired-end only, and cfDNAanalyzer fails rather than skipping them if handed single-end reads. On single-end data select only CNA, NOF, WPS or TSSC.

A selected feature whose reference files the chosen --genome does not supply is dropped from the selection rather than failing the run. Both bundled builds supply everything, so in practice every feature is available. cfDNAanalyzer supports hg19 and hg38 only, so a --genome config file of your own must set "cda_genome" to whichever of the two its coordinates match; without it every feature here is dropped, with a warning saying so.

Example: --cda-features CNA,EM,OCF
Example: --cda-features all

--cda-regions CDA_REGIONS

BED file of regions the region-specific cfDNAanalyzer features are measured over.
type: file
default: the genome build's TSS intervals

Applies to NOF, NP, WPS, OCF, EMR and FPR. Only the first three columns are used, and the regions are sorted before use, so any BED flavor will do; comment, track and browser lines are stripped.

Without this option the genome build's TSS interval file is used, which is every Ensembl TSS padded by 2 kb — roughly 62,000 regions. Ignored when --cda-features selects no region-specific feature.

Example: --cda-regions resources/promoters.bed

--cda-cna-bin-size CDA_CNA_BIN_SIZE

Bin size, in kilobases, for the cfDNAanalyzer CNA feature.
type: int
default: 1000

Restricted to 10, 50, 500 and 1000. cfDNAanalyzer reads pre-computed GC and mappability tracks for the requested bin size out of the ichorCNA it bundles, and those exist only at those four sizes, so any other value is rejected by the frontend rather than failing mid-run.

Unrelated to --bin-size, which bins the fragment length histogram. Ignored when --cda-features does not select CNA.

Example: --cda-cna-bin-size 500

2.3 Orchestration options¶

Each of the following arguments are optional, and do not need to be provided.

--dry-run

Dry run the pipeline.
type: boolean flag

Displays what steps in the pipeline remain or will be run. Does not execute anything!

Example: --dry-run

--silent

Silence standard output.
type: boolean flag

Reduces the amount of information directed to standard output when submitting master job to the job scheduler. Only the job id of the master job is returned.

Example: --silent

--mode {slurm,local}

Execution Method.
type: string
default: slurm

Execution Method. Defines the mode or method of execution. Vaild mode options include: slurm or local.

slurm
The slurm execution method will submit jobs to the SLURM workload manager. It is recommended running fragmentomics in this mode as execution will be significantly faster in a distributed environment. This is the default mode of execution.

local
Local executions will run serially on compute instance. This is useful for testing, debugging, or when a users does not have access to a high performance computing environment. If this option is not provided, it will default to a slurm execution mode.

Example: --mode slurm

--job-name JOB_NAME

Set the name of the pipeline's master job.
type: string
default: pipeline_fragmentomics

When submitting the pipeline to a job scheduler, like SLURM, this option always you to set the name of the pipeline's master job. By default, the name of the pipeline's master job is set to pipeline_fragmentomics.

Example: --job-name pl_id-42

--singularity-cache SINGULARITY_CACHE

Overrides the $SINGULARITY_CACHEDIR environment variable.
type: path
default: /path/to/output/directory/.singularity

Singularity will cache image layers pulled from remote registries. This ultimately speeds up the process of pull an image from DockerHub if an image layer already exists in the singularity cache directory. By default, the cache is set to the value provided to the --output argument. Please note that this cache cannot be shared across users. Singularity strictly enforces you own the cache directory and will return a non-zero exit code if you do not own the cache directory! See the --sif-cache option to create a shareable resource.

Example: --singularity-cache /data/$USER/.singularity

--sif-cache SIF_CACHE

Path where a local cache of SIFs are stored.
type: path

Uses a local cache of SIFs on the filesystem. This SIF cache can be shared across users if permissions are set correctly. If a SIF does not exist in the SIF cache, the image will be pulled from Dockerhub and a warning message will be displayed. The fragmentomics cache subcommand can be used to create a local SIF cache. Please see fragmentomics cache for more information. This command is extremely useful for avoiding DockerHub pull rate limits. It also remove any potential errors that could occur due to network issues or DockerHub being temporarily unavailable. We recommend running fragmentomics with this option when ever possible.

Example: --sif-cache /data/OpenOmics/SIFs

--threads THREADS

Max number of threads for local processes.
type: int
default: 2

Max number of threads for local process. This option is more applicable when running the pipeline with --mode local. It is recommended setting this vaule to the maximum number of CPUs available on the host machine.

Example: --threads 12

--tmp-dir TMP_DIR

Path for writing temporary files.
type: path
default: /lscratch/$SLURM_JOBID

Path on the file system for writing temporary output files. By default, the temporary directory is set to '/lscratch/$SLURM_JOBID' for backwards compatibility with the NIH's Biowulf cluster; however, if you are running the pipeline on another cluster, this option will need to be specified. Ideally, this path should point to a dedicated location on the filesystem for writing tmp files. On many systems, this location is set to somewhere in /scratch. If you need to inject a variable into this string that should NOT be expanded, please quote this options value in single quotes.

Example: --tmp-dir /data/scratch/$USER/

--overwrite-pipeline-template

Overwrite pipeline template in output directory.
type: boolean flag

Overwrite pipeline template files in an existing output directory. When this option is provided, the pipeline replaces the existing workflow/, resources/, and config/ directories in the output directory with fresh copies from the current pipeline installation. This is mainly useful for developers who want to re-run or test an existing output directory after updating the pipeline template or workflow.

Example: --overwrite-pipeline-template

2.4 Miscellaneous options¶

Each of the following arguments are optional, and do not need to be provided.

-h, --help

Display Help.
type: boolean flag

Shows command's synopsis, help message, and an example command

Example: --help

3. Example¶

3.1 Running the pipeline with input BAM files¶

# Step 1.) Grab an interactive node,
# do not run on head node!
srun -N 1 -n 1 --time=1:00:00 --mem=8gb  --cpus-per-task=2 --pty bash
module purge
module load singularity snakemake

# Step 2A.) Dry-run the pipeline with
# Input BAM files
./fragmentomics run --input .tests/*.bam \
    --sif-cache /data/OpenOmics/SIFs \
    --output /data/$USER/output \
    --mode slurm \
    --dry-run

# Step 2B.) Run the fragmentomics pipeline
# The slurm mode will submit jobs to 
# the cluster. It is recommended running 
# the pipeline in this mode.
./fragmentomics run --input .tests/*.bam \
    --sif-cache /data/OpenOmics/SIFs \
    --output /data/$USER/output \
    --mode slurm

3.2 Running the pipeline with input FastQ files¶

# Step 1.) Grab an interactive node,
# do not run on head node!
srun -N 1 -n 1 --time=1:00:00 --mem=8gb  --cpus-per-task=2 --pty bash
module purge
module load singularity snakemake

# Step 2A.) Dry-run the pipeline with
# Input paired-end FastQ files, please
# note single-end data is not supported!
./fragmentomics run \
    --input .tests/*.fastq.gz \
    --output /data/$USER/output \
    --genome hg38 \
    --sif-cache /data/OpenOmics/SIFs \
    --mode slurm \
    --dry-run

# Step 2B.) Run the fragmentomics pipeline
# The slurm mode will submit jobs to 
# the cluster. It is recommended running 
# the pipeline in this mode.
./fragmentomics run \
    --input .tests/*.fastq.gz \
    --output /data/$USER/output \
    --genome hg38 \
    --sif-cache /data/OpenOmics/SIFs \
    --mode slurm

4. Results¶

The --output directory holds one subdirectory per analysis, plus project-level QC. The full layout is documented in Pipeline overview §4; the files to look at first are:

multiqc/multiqc_report.html

Aggregate QC across all samples. Alignment rates, duplicate rates, read quality, insert-size distributions and per-contig read counts, plus a Fragmentomics section summarizing the analysis results themselves — fragment lengths, end motifs, MDS, DELFI and the TSS profiles — for every sample side by side. Written for either input type. See Quality control.

coverage/coverage_summary.xlsx

Merged coverage workbook. A per-sample summary sheet and the full interval × sample coverage matrix.

frag_length_bins/{sample}_frag_bin1.png

Fragment-length histogram. The fastest per-sample sanity check — look for the ~167 bp mononucleosome peak.

qc/{sample}.contig_validation.txt

Contig validation record (BAM input only). Which contigs were kept, which were dropped, and how many reads were affected. See Reference contig filter.

Each analysis output is documented on its own page under Analyses.