This is the first post in a series on processing an ancient DNA sample for use with ADMIXTOOLS. Here I go from paired-end FASTQ files to a filtered, duplicate-removed BAM aligned to hs37d5. The workflow is based on the run I used for ERR14088885, from the Başur Höyük study PRJEB83032.

Ancient DNA needs a different alignment strategy from ordinary modern whole-genome data. The molecules are short, the ends may carry post-mortem damage, and paired reads often overlap because the DNA insert is shorter than the sequencing cycles. For this sample I therefore clean poly-G tails, trim adapters, merge overlapping mates, align the merged molecules with bwa aln, remove low-confidence alignments, and deduplicate using both observed ends of each molecule.

The commands below run on any Debian or Ubuntu Linux machine.

Hardware and software

The example run used 12 threads and 31 GiB of RAM. A machine with fewer cores will also work, but alignment will take longer. The samtools sort -m value is memory per thread: 1500M with 12 threads permits sorting to use about 18 GiB, leaving room for the operating system and the other tools.

If you’re running this on Windows Subsystem for Linux, work exclusively within the Linux filesystem (~/ or /home/username/), not in Windows directories such as /mnt/c/. Accessing the Windows filesystem from WSL incurs substantial I/O overhead and can make alignment extremely slow. Keep all downloads, reference genomes, and BAM files in your Linux home directory.

Install the required packages:

sudo apt-get update
sudo apt-get install -y \
  adapterremoval bwa curl default-jre-headless fastp samtools wget

The recorded run used fastp 0.24.0, AdapterRemoval 2.3.4, bwa 0.7.18, samtools 1.21, and DeDup 0.12.9.

Set up the working directories

Run the pipeline from a project directory with enough free space for the reference, FASTQs, and intermediate files:

set -Eeuo pipefail

PROJECT_DIR=$PWD
SAMPLE=ERR14088885
DATA_DIR=${PROJECT_DIR}/data/${SAMPLE}
REFERENCE_DIR=${PROJECT_DIR}/hs37d5
OUTPUT_DIR=${PROJECT_DIR}/results/${SAMPLE}_paired_merged
POLYG_DIR=${OUTPUT_DIR}/fastp_polyg
TRIM_DIR=${OUTPUT_DIR}/adapterremoval
TMP_DIR=${OUTPUT_DIR}/tmp
TOOL_DIR=${PROJECT_DIR}/tools

THREADS=12
SORT_MEMORY=1500M
JAVA_MEMORY=20G

mkdir -p \
  "$DATA_DIR" "$REFERENCE_DIR" "$OUTPUT_DIR" \
  "$POLYG_DIR" "$TRIM_DIR" "$TMP_DIR" "$TOOL_DIR"

Adjust THREADS, SORT_MEMORY, and JAVA_MEMORY for your machine.

Download and index hs37d5

I use hs37d5 because the downstream target dataset is AADR, which uses GRCh37 coordinates. The decoy sequences in hs37d5 also give reads from repetitive or non-reference regions somewhere more appropriate to align than the primary chromosomes.

Download the compressed reference and verify it before indexing:

cd "$REFERENCE_DIR"

wget -c \
  https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz

gzip -t hs37d5.fa.gz
echo 'd10eebe06c0dbbcb04253e3294d63efc  hs37d5.fa.gz' | md5sum -c -

# Keep the downloaded .gz file and write an uncompressed FASTA for the tools.
gzip -dc hs37d5.fa.gz > hs37d5.fa

Build the BWA and samtools indexes. This is a one-time step for the reference:

bwa index hs37d5.fa
samtools faidx hs37d5.fa

test -s hs37d5.fa.fai
for suffix in amb ann bwt pac sa; do
  test -s "hs37d5.fa.${suffix}" || {
    echo "Missing BWA index: hs37d5.fa.${suffix}" >&2
    exit 1
  }
done

Set the reference path for the remaining commands:

REFERENCE=${REFERENCE_DIR}/hs37d5.fa

Download and check the paired FASTQs

ENA provides ERR14088885 as two compressed FASTQ files containing the paired-end reads: R1 and R2.

cd "$DATA_DIR"

wget -c \
  ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR140/085/ERR14088885/ERR14088885_1.fastq.gz \
  ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR140/085/ERR14088885/ERR14088885_2.fastq.gz

R1=${DATA_DIR}/${SAMPLE}_1.fastq.gz
R2=${DATA_DIR}/${SAMPLE}_2.fastq.gz

gzip -t "$R1"
gzip -t "$R2"

Remove poly-G tails

The reads were produced on a NovaSeq instrument. NovaSeq uses two-color chemistry, where a lack of signal is read as G; once a short molecule has been sequenced through, this can produce artificial poly-G tails. A scan of the first 100,000 pairs found terminal runs of at least ten Gs in 69.7% of R1 and 66.3% of R2.

I use fastp only for poly-G trimming. Adapter trimming and quality handling are left to AdapterRemoval in the next step:

POLYG_R1=${POLYG_DIR}/${SAMPLE}_1.polyg.fastq.gz
POLYG_R2=${POLYG_DIR}/${SAMPLE}_2.polyg.fastq.gz

fastp \
  --in1 "$R1" \
  --in2 "$R2" \
  --out1 "$POLYG_R1" \
  --out2 "$POLYG_R2" \
  --thread "$THREADS" \
  --trim_poly_g \
  --poly_g_min_len 10 \
  --disable_adapter_trimming \
  --disable_quality_filtering \
  --disable_length_filtering \
  --dont_eval_duplication \
  --json "${POLYG_DIR}/${SAMPLE}.fastp.json" \
  --html "${POLYG_DIR}/${SAMPLE}.fastp.html"

test -s "$POLYG_R1"
test -s "$POLYG_R2"

The input contained 32,987,452 read pairs. After poly-G cleanup, 32,949,986 pairs were written for the next stage. Mean read lengths changed from 150 bases for both mates to 122 bases for R1 and 131 bases for R2.

Trim adapters and merge overlapping mates

Short aDNA inserts are often sequenced from both ends with extensive overlap. AdapterRemoval can trim the adapter sequence and combine the two observations into one consensus read. This prevents an overlapping pair from counting the same molecule twice.

AR_PREFIX=${TRIM_DIR}/${SAMPLE}
COLLAPSED_FASTQ=${AR_PREFIX}.collapsed.gz

AdapterRemoval \
  --file1 "$POLYG_R1" \
  --file2 "$POLYG_R2" \
  --basename "$AR_PREFIX" \
  --threads "$THREADS" \
  --gzip \
  --trimns \
  --trimqualities \
  --minquality 20 \
  --minlength 30 \
  --collapse \
  --collapse-deterministic \
  --preserve5p

test -s "$COLLAPSED_FASTQ"

These choices have specific jobs:

  • --minlength 30 discards molecules shorter than 30 bases after trimming.
  • --collapse merges mates with at least 11 bases of overlap, the AdapterRemoval default recorded in this run.
  • --collapse-deterministic writes N with quality 0 when overlapping bases disagree and have equal quality, instead of resolving the tie randomly.
  • --preserve5p protects the physical molecule ends used for aDNA duplicate detection. With AdapterRemoval 2.3.4, it also means collapsed reads are not quality-trimmed at either end.

The .settings report is important: read it before deciding whether a merged-only pipeline is suitable for another library.

less "${AR_PREFIX}.settings"

For ERR14088885, AdapterRemoval processed 32,949,986 pairs and produced 31,355,473 full-length collapsed reads, a collapse rate of about 95.2%. The retained reads averaged 67 bases. Because this workflow sends only ${SAMPLE}.collapsed.gz to BWA, non-overlapping pairs and singletons do not enter the final BAM. That is a deliberate, conservative choice for this library, not a general rule for every aDNA dataset. If the collapse rate is low, align the non-collapsed pairs and singletons as separate branches rather than silently discarding most of the data.

Align the collapsed molecules with BWA

The original version of this tutorial used bwa mem -k 16 -L 16,16 directly on the raw paired FASTQs. I replaced that step with bwa aln, using the settings from the completed run:

SAI=${OUTPUT_DIR}/${SAMPLE}.collapsed.sai

bwa aln \
  -t "$THREADS" \
  -l 1024 \
  -n 0.01 \
  -o 2 \
  "$REFERENCE" \
  "$COLLAPSED_FASTQ" \
  > "$SAI"

test -s "$SAI"

The parameters are commonly used for short, damaged ancient DNA:

  • -l 1024 makes the seed longer than the reads, effectively disabling seeding. Damage near a read end can then contribute to the full-read alignment instead of causing the seed to fail.
  • -n 0.01 is BWA’s floating-point edit-distance setting. BWA chooses the maximum edit distance by read length, using a 1% missing-alignment threshold under its error model.
  • -o 2 allows up to two gap openings instead of the default one.

On 12 threads, bwa aln processed the 31,355,473 collapsed reads in 2 hours 43 minutes. Runtime will vary with CPU speed, storage, read length, and sample size.

Convert the SAI output to SAM, add a read group, retain mapped reads with MAPQ at least 30, and coordinate-sort the result:

FILTERED_BAM=${OUTPUT_DIR}/${SAMPLE}.mapq30.sorted.bam
READ_GROUP="@RG\\tID:${SAMPLE}\\tSM:${SAMPLE}\\tLB:${SAMPLE}\\tPL:ILLUMINA"

bwa samse -r "$READ_GROUP" \
    "$REFERENCE" "$SAI" "$COLLAPSED_FASTQ" \
  | samtools view --bam --with-header -F 4 -q 30 - \
  | samtools sort \
      -@ "$THREADS" \
      -m "$SORT_MEMORY" \
      -T "${TMP_DIR}/${SAMPLE}.sort" \
      -o "$FILTERED_BAM" -

test -s "$FILTERED_BAM"
samtools quickcheck -v "$FILTERED_BAM"

Here, -F 4 removes unmapped records and -q 30 removes alignments below MAPQ 30. Filtering before duplicate removal keeps DeDup focused on alignments that will be used downstream.

This sample produced 219,005 mapped reads at MAPQ 30 or higher. That is only about 0.7% of the collapsed input, so it should not be mistaken for a typical mapping rate. A low rate like this deserves follow-up QC: confirm the sample identity and reference, inspect read composition and adapter reports, and quantify endogenous human DNA before drawing biological conclusions.

Remove PCR duplicates

I use DeDup because it was designed for short ancient DNA and can compare both ends of a merged molecule. Download the pinned release once:

DEDUP_VERSION=0.12.9
DEDUP_JAR=${TOOL_DIR}/DeDup-${DEDUP_VERSION}.jar

if [[ ! -s "$DEDUP_JAR" ]]; then
  curl --fail --location --retry 3 \
    --output "${DEDUP_JAR}.partial" \
    "https://github.com/apeltzer/DeDup/releases/download/${DEDUP_VERSION}/DeDup-${DEDUP_VERSION}.jar"
  mv "${DEDUP_JAR}.partial" "$DEDUP_JAR"
fi

java -jar "$DEDUP_JAR" -h >/dev/null

Run it in merged-read mode:

java "-Xmx${JAVA_MEMORY}" -jar "$DEDUP_JAR" \
  -m \
  -i "$FILTERED_BAM" \
  -o "$OUTPUT_DIR"

DEDUP_BAM=${OUTPUT_DIR}/${SAMPLE}.mapq30.sorted_rmdup.bam
test -s "$DEDUP_BAM"

The -m flag tells DeDup that every input record is a merged molecule, so read-name prefixes are not required. Do not apply it to an ordinary single-end library: the observed 3′ read end is not necessarily the original molecule end in that case.

For this run, DeDup examined 219,005 mapped reads, removed 100,188 duplicates, and retained 118,817 reads. The reported duplication rate was 0.46.

Write and verify the final BAM

Sort DeDup’s output once more, index it, and write two standard samtools reports:

FINAL_BAM=${OUTPUT_DIR}/${SAMPLE}.final.bam

samtools sort \
  -@ "$THREADS" \
  -m "$SORT_MEMORY" \
  -T "${TMP_DIR}/${SAMPLE}.final-sort" \
  -o "$FINAL_BAM" \
  "$DEDUP_BAM"

samtools index -@ "$THREADS" "$FINAL_BAM"
samtools flagstat -@ "$THREADS" "$FINAL_BAM" \
  > "${OUTPUT_DIR}/${SAMPLE}.flagstat.txt"
samtools stats -@ "$THREADS" "$FINAL_BAM" \
  > "${OUTPUT_DIR}/${SAMPLE}.stats.txt"

samtools quickcheck -v "$FINAL_BAM"
test -s "${FINAL_BAM}.bai"
cat "${OUTPUT_DIR}/${SAMPLE}.flagstat.txt"

The files to carry forward are:

results/ERR14088885_paired_merged/ERR14088885.final.bam
results/ERR14088885_paired_merged/ERR14088885.final.bam.bai

Result

The final BAM contains 118,817 mapped reads after MAPQ filtering and duplicate removal, including 237 reads on chromosome Y. PMDtools 0.60 gave a mean PMD score of -0.745 across the BAM, with 1,027 reads (0.86%) reaching PMD >= 3; none of the Y reads reached that threshold. The single reads supporting M694/CTS5611 and CTS7400/PF6469 scored -0.797 and -1.403, respectively, under the default double-stranded model.

For practical purposes, this BAM is unusable for reliable autosomal or Y-DNA analysis. Conservative PMD, end, and substitution filtering leaves only 22 v66 SNPs, all at depth one and without independent read support. YBYRA found an R1b-M269-like best-scoring path, but its score remained below the reporting threshold, so no formal Y-DNA haplogroup could be assigned. The original study also excluded this sample from downstream genetic analysis.