Skip to content

Step 2: Alignment (FASTQ to BAM)

What This Does

Aligns raw sequencing reads against the GRCh38 human reference genome and marks duplicate reads. Produces a sorted, duplicate-marked, indexed BAM file, and checks that the BAM holds every read that went in.

Why

Alignment maps each 150bp sequencing read to its position in the human genome. Required for all downstream variant calling.

Duplicates are copies of one DNA fragment (PCR or optical copies). They are not independent evidence. DeepVariant copes with them, but GATK, FreeBayes, Octopus and the structural-variant and depth steps would count each copy as another read. samtools markdup sets the duplicate flag (0x400) on all copies but one, and those tools then skip them. No read is removed.

Tools

  • minimap2: fast aligner (preferred for WGS)
  • samtools: fixmate, sort and markdup (markdup also writes the index)

Docker Images

  • MINIMAP2_IMAGE (minimap2 aligner)
  • SAMTOOLS_IMAGE (samtools fixmate, sort, markdup with its index)
  • MINIMAP2_SAMTOOLS_IMAGE (the Nextflow pipeline only: its task ALIGN_MINIMAP2 runs the same pipe in one container, so it uses one BioContainers image that holds minimap2 and samtools at the releases of the two images above)

Pinned in versions.env; Image versions lists the current tag.

Prerequisites

  • GRCh38 reference genome (${REF_FASTA}, see reference setup)
  • minimap2 index (.sr.mmi file next to the reference, generated once by the script)
  • Paired-end FASTQ files

Commands

./scripts/02-alignment.sh your_sample        # THREADS=16 to use 16 CPUs (default 8)

What the script runs, written as plain commands:

REF_FASTA=reference/GRCh38_no_alt_analysis_set.fasta   # see 00-reference-setup.md#the-reference-path-on-every-page
SAMPLE=your_sample
GENOME_DIR=/path/to/your/data
REF="${GENOME_DIR}/${REF_FASTA}"
MMI="${REF%.fasta}.sr.mmi"
OUT="${GENOME_DIR}/${SAMPLE}/aligned/${SAMPLE}_sorted.bam"

# Step 1: Create the minimap2 index with the short-read preset (one-time, 1-3 min measured).
# The preset sets the k-mer and window size of the index (k21, w11 for -x sr);
# a plain `minimap2 -d` builds k15, w10, which does not match `-x sr` mapping.
minimap2 -x sr -d "$MMI" "$REF"

# Step 2: Align, mark duplicates and sort (3 h 51 min - 5 h 21 min measured on 30-40x WGS).
# fixmate -m adds the mate tags markdup needs; it reads the pairs minimap2
# writes next to each other, then sort orders by position for markdup.
# markdup --write-index builds the index while it writes the BAM; the
# ##idx##<name>.bai suffix makes it a .bai (without it samtools writes a .csi).
minimap2 -a -x sr -t 16 \
  -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA\tLB:${SAMPLE}" \
  "$MMI" \
  ${GENOME_DIR}/${SAMPLE}/fastq/${SAMPLE}_R1.fastq.gz \
  ${GENOME_DIR}/${SAMPLE}/fastq/${SAMPLE}_R2.fastq.gz \
| samtools fixmate -@ 16 -u -m - - \
| samtools sort -u -@ 16 -m 1G - \
| samtools markdup -@ 16 --write-index - "${OUT%.bam}.tmp.bam##idx##${OUT%.bam}.tmp.bam.bai"

# Step 3: Check the BAM. test -s stops if markdup wrote no index.
samtools quickcheck "${OUT%.bam}.tmp.bam"
test -s "${OUT%.bam}.tmp.bam.bai"

# Step 4: Every read that went in is in the BAM: its primary records (not
# secondary 0x100, not supplementary 0x800) equal the FASTQ records of R1
# plus R2. The script counts the FASTQ while minimap2 reads it.
IN=$(( $(gzip -cd ${GENOME_DIR}/${SAMPLE}/fastq/${SAMPLE}_R1.fastq.gz \
         ${GENOME_DIR}/${SAMPLE}/fastq/${SAMPLE}_R2.fastq.gz | wc -l) / 4 ))
PRIMARY=$(samtools view -c -F 0x900 "${OUT%.bam}.tmp.bam")
[ "$PRIMARY" -eq "$IN" ] || { echo "the BAM lost $((IN - PRIMARY)) reads"; exit 1; }

# Step 5: Rename into place. markdup finishes the index just before the
# BAM's last block, so touch makes the index the newer file.
touch "${OUT%.bam}.tmp.bam.bai"
mv "${OUT%.bam}.tmp.bam" "$OUT" && mv "${OUT%.bam}.tmp.bam.bai" "${OUT}.bai"

# Output: ~80-120 GB BAM + ~9 MB BAI index

The script writes the BAM and the index under temporary names, both in the one pass of samtools markdup, and renames them only after samtools quickcheck and the read count pass, so a run that is killed never puts a partial BAM under ${SAMPLE}_sorted.bam. A BAM and index from an earlier complete run stay in place, and later steps still read them. run-all.sh starts the pipeline from that BAM when it and its .bai exist (and validate-setup.sh checks its header against the reference); with no BAM it starts from the FASTQ, and the pipeline aligns them. An index built by an older version (reference/GRCh38.mmi, default preset) is no longer used and can be deleted.

Read count

minimap2 and bwa-mem2 write every read they are given once, mapped or not, and samtools fixmate, sort and markdup drop none. So the finished BAM holds exactly as many primary records as R1 and R2 together hold FASTQ records. The unit is reads, not pairs: both reads of a pair count, which is also what fastp's summary.after_filtering.total_reads counts.

Steps 02 and 02a count the records of the two FASTQ files they align, with one gzip per file running beside the aligner, and count the BAM's primary records with samtools view -c -F 0x900 once markdup is done. That second count reads the whole BAM once more. When the two differ, the step exits 1, keeps no new BAM, and leaves an earlier BAM in place. The message says how many reads are missing and how many records each FASTQ had. The Nextflow task ALIGN_MINIMAP2 applies the same rule and fails, without a retry.

Both counts go to qc/<sample>_read_count.tsv (step 02a: qc/<sample>_bwamem2_read_count.tsv), a key/value table with input_reads, input_source, primary_records, missing_reads and read_count_check (ok or mismatch). The Nextflow pipeline publishes the same table under <outdir>/<sample>/qc/. The text report (generate-report.sh) and the HTML report of step 24 show the two counts in their quality-control block, and the missing reads in red. The pipeline's own HTML_REPORT task does not get the table yet, so the Nextflow HTML report leaves the counts out.

One way to lose reads without an error: when R1 and R2 hold different numbers of records, minimap2 and bwa-mem2 stop at the end of the shorter file with only a warning. The e2e case align-read-count-gate-1 cuts R2 of a 20,000-pair slice to 12,000 records: minimap2 writes 24,000 records of the 32,000 reads, and step 02 and ALIGN_MINIMAP2 stop.

For an input you know is odd, READ_COUNT_CHECK=warn keeps the BAM with a warning, and the table says mismatch. In the Nextflow pipeline the same switch is ext.read_count_check = 'warn' on ALIGN_MINIMAP2, in a config file given with -c (Troubleshooting).

validate-setup.sh runs the same comparison for an existing BAM when fastp's report sits beside it (fastq_trimmed/<sample>_fastp.json): the BAM's primary records against fastp's total_reads, an [OK] or a [FAIL] line ([WARN] under READ_COUNT_CHECK=warn). That covers a BAM aligned before this check. Counting reads the whole BAM, a few minutes on a 30x one. A BAM aligned from the untrimmed reads (FASTQ_SUBDIR=fastq) is not comparable with fastp's count: when qc/<sample>_read_count.tsv, not older than the BAM, says its input_source is fastq/..., validate-setup.sh skips the comparison with an [INFO] line, because step 02 already compared that BAM with its own input. A BAM from the raw reads with no such table, aligned before this check, still gets the [FAIL]: realign it with step 02, or go on with READ_COUNT_CHECK=warn.

Resource Requirements

  • CPU: 16+ cores recommended (THREADS)
  • RAM: minimap2 loads the whole index into memory. With the sr index it peaked at 10 GB on the test reference (1.8 Gb, 57% of GRCh38), against 5.5 GB with the old default index. On GRCh38 the whole task, minimap2 and samtools sort with 1 GB per thread, peaked at 18.3-18.4 GB in two measured 30-40x runs (runtime); the container's cap is 32 GB.
  • Disk and time: see Hardware and storage requirements (the BAM is about 80-120 GB). Sort spills go to aligned/${SAMPLE}.sort_tmp/ in the sample directory and are removed at the end.

Notes

  • Use -x sr for Illumina short reads (short-read preset), for the index and for mapping
  • Alternative: bwa-mem2 is equally valid but minimap2 is faster (see scripts/02a-alignment-bwamem2.sh). It writes aligned_bwamem2/ with the same fixmate, sort and markdup steps, piped, with no SAM file on disk. Its one-time index build needs about 90 GB of RAM for GRCh38 (28 GB per Gbp of reference, see reference setup); the script stops with that message when the build is killed for lack of memory.
  • The BAM index (.bai) must always accompany the BAM file