Pipeline overview¶

1. About¶

This section documents what the pipeline does once fragmentomics run has been invoked. Where the Commands section documents the command-line interface, this section documents the workflow behind it: how inputs are normalized, which reference files are used, and what each analysis step computes.

Every step runs inside a Singularity container defined in config/containers.json, and every step draws its threads, memory, walltime, partition and --gres request from config/cluster.json. No step depends on environment modules, so the workflow behaves identically under --mode slurm and --mode local.

2. Data flow¶

The pipeline accepts either paired-end Illumina FastQ files or ready-made alignments in BAM format, never both in the same run. The input type is auto-detected by the frontend and recorded in config.json as project.input_type. Both input types converge on a single canonical per-sample file:

bams/{sample}.sorted.bam

Every analysis step reads that file and nothing else, which is what makes the analysis half of the pipeline input-type agnostic.

inputs/{sample}.R1.fastq.gz
inputs/{sample}.R2.fastq.gz
         │
         │  align_fastq   (fastqc → fastp trim+filter → bwa-mem2 → filter
         │                 → name sort → fixmate → coordinate sort
         │                 → markdup → fastqc)
         ▼
bams/{sample}.sorted.bam

<user-provided>.bam
         │
         │  stage_bams              (filter → coordinate sort + index)
         ▼
staged_bams/{sample}.sorted.bam
         │
         │  remove_orphan_reads     (drop singletons and orphaned mates)
         ▼
staged_bams/{sample}.paired.bam
         │
         │  filter_reference_contigs  (subset to reference contigs)
         ▼
bams/{sample}.sorted.bam

From the canonical BAM, the analysis steps fan out in parallel:

                       ┌─ coverage ───────────────► merge_coverage_excel
                       ├─ frag-length-bins
                       ├─ frag-length-intervals
                       ├─ end-motifs ─────────────► mds
bams/{sample}.sorted.bam ─┼─ interval-end-motifs
                       ├─ bam_to_bed
                       ├─ delfi
                       ├─ wps ──┬─────────────────► agg-bw  (aggregate)
                       │        └─ adjust-wps ────► agg-bw  (aggregate)
                       ├─ cleavage-profile ───────► agg-bw  (aggregate)
                       ├─ bam_stats ──────────────► multiqc
                       ├─ fastqc_bam ─────────────► multiqc   (BAM input only)
                       └─ fastp_bam ──────────────► multiqc   (BAM input only)

Each analysis step is documented in the Analyses section.

The analysis outputs are not only end products. finaletoolkit_multiqc gathers them — coverage, fragment lengths, end motifs, MDS, DELFI and the aggregate TSS profiles — and summarizes them into the aggregate report as well, so a run's fragmentomics measurements can be compared across samples in one place. The two read-QC steps above exist only for BAM input; a FastQ run gets the equivalent reports from align_fastq, which runs fastqc and fastp as part of aligning. See Quality control.

3. Which steps run¶

Some steps are conditional on the reference files the selected build provides. The workflow only requests their output if the files they need are present in config/genome.json, so a partially-populated build skips those analyses instead of failing on a missing path.

The three window-based analyses are gated on reference_fa because the genomic interval BED they share is tiled from it at run time by make_intervals, at the width given by --interval, rather than read from a pre-made file. That job runs once per pipeline invocation, ahead of the analyses that consume it.

The remaining steps are unconditional: frag-length-bins, coverage, bam_to_bed (it needs no reference file beyond the BAM), and the project-level QC steps (bam_stats, merge_coverage_excel, multiqc) always run. Note that coverage needs the build's tss_interval file even though it is not gated on it, which is why that is the one reference key a build cannot leave out — a --genome config file without it is rejected before the run starts.

Both bundled builds define every reference key, so a default hg38 or hg19 run executes every step. A build of your own runs whichever subset of the table above its reference files cover.

4. Output directory layout¶

<--output>/
├── inputs/                     symlinks to the user's input files (read-only)
├── staged_bams/                sorted BAMs awaiting pairing cleanup and the
│                               contig filter (BAM input only)
├── bams/                       canonical analysis BAM + index, one per sample
├── intervals/                  windows tiled from the reference FastA at
│   └── {genome}_{size}_intervals.bed   --interval width, one BED per run
├── beds/
│   ├── {sample}.bed.gz         aligned intervals of the analysis BAM, bgzip -l 9
│   └── {sample}.bed.gz.tbi     tabix index for random access by locus
├── coverage/
│   ├── {sample}_coverage.bed
│   └── coverage_summary.xlsx   merged workbook across all samples
├── frag_length_bins/
│   ├── {sample}_frag_bin{bin_size}.tsv
│   └── {sample}_frag_bin{bin_size}.png
├── frag_length_intervals/{sample}_frag_interval.bed
├── end_motifs/{sample}_endmotif.tsv
├── interval_end_motifs/{sample}_endmotif_interval.tsv
├── mds/{sample}_mds.tsv
├── delfi/{sample}_delfi.bed
├── wps/
│   ├── {sample}_wps_out_tss.bw
│   └── {sample}_wps_out_tss_aggr.wig
├── adjust_wps/
│   ├── {sample}_wps_out_tss_adjusted.bw
│   └── {sample}_wps_out_tss_adj_aggr.wig
├── cleavage_profile/
│   ├── {sample}_cleavage_profile_tss.bw
│   └── {sample}_cleavage_profile_aggr.wig
├── qc/                         per-sample samtools, FastQC and fastp reports,
│                               contig validation reports, and finaletoolkit/
│                               (the fragmentomics summary for MultiQC)
├── multiqc/multiqc_report.html
├── config/                     resolved configuration for this run
├── workflow/                    Snakemake rules and scripts for this run
└── logfiles/                   master job log and per-job logs (slurm mode)

Sample names ({sample}) are derived from the input file basenames with the .R1/.R2 mate suffix or the .bam extension stripped. Two inputs that reduce to the same sample name are rejected up front rather than silently clobbering one another.

5. Reproducibility¶

The --output directory is self-contained. config/, workflow/ and the resolved config.json are copied into it at launch, so a completed output directory records exactly which rules, scripts, reference paths and resource requests produced its results. Re-running in the same directory resumes rather than restarts; passing --overwrite-pipeline-template refreshes the copied template from the current installation.

See also: Reference files and Quality control.