1. Alignment

The align command aligns preprocessed paired-end FASTQ files to a reference genome and produces cleaned, deduplicated, indexed BAM files for downstream QC and peak-level analysis.

The command automatically processes all paired FASTQ files detected in the input directory. Each FASTQ pair is aligned independently and produces one final BAM file.

For each detected paired-end FASTQ pair, the pipeline performs the following steps:

  • Reads are aligned to the reference genome using bowtie2 in local, very-sensitive-local mode, with --sam-append-comment enabled. Fragment size constraints are applied using the Bowtie2 paired-end insert size filters. Any CB/UMI values written by the upstream adapt step as SAM-compliant optional fields (CB:Z:..., UB:Z:...) in the read comment are carried directly into the BAM as read tags, with no separate extraction step required.

  • The resulting alignments are converted to BAM format, filtered by mapping quality, and cleaned by removing mitochondrial and non-standard contigs. The filtered alignments are then sorted by genomic coordinates.

  • Duplicate removal is performed after alignment. All deduplication strategies always consider genomic alignment coordinates as the primary criterion. Depending on the adapter structure used upstream, additional tag-aware grouping may be applied.

    • If reads carry a CB tag (from cell barcode information identified upstream), duplicate identification becomes CB-aware, meaning that reads with identical coordinates but originating from different cell barcodes are not treated as duplicates.

    • If reads carry a UB tag (from UMI information identified upstream), duplicate identification becomes UMI-aware, meaning that reads with identical coordinates but different UMI values are preserved as independent molecules.

    • When both CB and UMI tags are present, duplicates are identified using the combined key of genomic coordinates, CB, and UMI, ensuring that deduplication respects both cellular origin and molecular identity.

  • After duplicate removal, the coordinate-sorted BAM file is written to disk and indexed using samtools index. If indexing fails due to sorting inconsistencies, the file is resorted and indexed again to guarantee a valid output.

  • If --merge-crf is enabled, BAM files are merged by CRF target into a single per-CRF BAM after all pairs are processed. Each pair BAM (e.g. A-B.bam) contributes to both CRF A and CRF B; symmetric pairs (e.g. A-A.bam) contribute twice to CRF A, preserving the read count symmetry. Merged BAMs are coordinate-sorted and indexed. The original per-pair BAMs are retained under the pair subdirectory. Note that CRF merging is currently experimental and not recommended for routine use.

Parameters

Argument Type Default Description Example
-i, –input character — Directory containing paired-end FASTQ files (.fastq.gz) generated from the adapter identification step. -i ./adapter
-o, –output character — Output directory where final BAM files will be written (created if missing). -o ./bam
-g, –genome character — Reference genome name used to resolve the Bowtie2 index (hg38 or mm10). -g hg38
–cb flag off Indicates that reads contain cell barcode (CB) information encoded in the read name.
If enabled, deduplication becomes cell-barcode-aware, meaning reads with identical genomic coordinates but different CB values are preserved. 		–cb
–umi flag off Indicates that reads contain UMI information encoded in the read name.
If enabled, deduplication becomes UMI-aware, meaning reads with identical genomic coordinates but different UMIs are preserved as independent molecules. 		–umi
–min-len integer 10 Minimum fragment length accepted by Bowtie2 paired-end alignment. Pairs with inferred insert size below this threshold are discarded. –min-len 10
–max-len integer 800 Maximum fragment length accepted by Bowtie2 paired-end alignment (-X). Pairs with inferred insert size above this threshold are discarded. –max-len 1200
–merge-crf flag off Enable per-CRF BAM merging. If set, pair BAMs are merged into per-CRF BAMs in the output directory, with pair BAMs retained under pair subdirectory. By default, merging is skipped and only per-pair BAMs are produced. –merge-crf
-j, –threads integer auto-detect Number of CPU threads used by Bowtie2 and samtools operations. -j 16
–java-mem character auto-detect Java heap size passed to Picard (e.g., 24g, 8000m). Automatically detected if not provided. –java-mem 24g
–picard-jar character auto-detect Path to picard.jar. If missing, the script attempts to locate it automatically in the conda environment. –picard-jar /path/to/picard.jar

Output Files

By default (CRF merging disabled):

  1. Pair BAM & Index - OUT_DIR/<crf-pair>.bam and OUT_DIR/<crf-pair>.bam.bai
    • Coordinate-sorted, filtered, and deduplicated BAM file for each pair. No per-CRF merging is performed.

With --merge-crf:

  1. CRF BAM & Index - OUT_DIR/<crf>.bam and OUT_DIR/<crf>.bam.bai
    • Coordinate-sorted, filtered, and deduplicated BAM file.
  2. Pair BAM & Index - OUT_DIR/pair/<crf-pair>.bam and OUT_DIR/pair/<crf-pair>.bam.bai
    • Original pair BAM files retained under the pair subdirectory.

Example Usage

# Regular Hiplex CUT&Tag
multiEpiPrep align \
  -i ./adapter \
  -o ./bam \
  -g hg38

# UMI-containing Hiplex CUT&Tag
multiEpiPrep align \
  -i ./adapter \
  -o ./bam \
  -g hg38 \
  --umi

# CB + UMI aware deduplication
multiEpiPrep align \
  -i ./adapter \
  -o ./bam \
  -g hg38 \
  --cb \
  --umi

# Enable CRF merging (produce CRF BAMs)
multiEpiPrep align \
  -i ./adapter \
  -o ./bam \
  -g hg38 \
  --merge-crf