# SA5 and SA6 pseudohaploid VCFs

This folder contains pseudohaploid autosomal genotype calls for two Late
Pleistocene aurochs (*Bos primigenius*) from Šandalja II, Istria, Croatia:

- `SA5.vcf.gz` (+ `SA5.vcf.gz.tbi`) — sample HRV_SA5_11612 (11,361–11,825 cal BP)
- `SA6.vcf.gz` (+ `SA6.vcf.gz.tbi`) — sample HRV_SA6_14757 (14,163–14,847 cal BP)

## Site list (Rossi et al. 2024 aurochs panel)

Calling was restricted to positions from the published Rossi et al. (2024)
aurochs panel (Nature 635:136), kindly provided by Conor Rossi (personal
communication, 2025) as
`chr1-29-bos_QUAL20_allelefilter_noCpG_named_geno0.2_mind0.2_singletonhet_singletonhomo_repetitive_trial_transversions_filter_Rossi.vcf.gz`.
The filename encodes the upstream filtering applied by the Rossi pipeline to
the joint variant set of 21 ancient *Bos primigenius* genomes:

- restricted to chromosomes 1–29 of *Bos taurus* (autosomes)
- variant quality QUAL ≥ 20
- biallelic allele filter
- CpG sites excluded
- variants assigned IDs (`named`)
- sites with > 20% missing genotypes excluded (`geno 0.2`)
- individuals with > 20% missing data excluded (`mind 0.2`)
- singleton heterozygotes and singleton homozygotes excluded
- repetitive regions excluded
- triallelic sites excluded (`trial`)
- transversions only (transitions excluded, to mitigate ancient-DNA C→T / G→A damage)

The resulting panel comprises **7,108,856 biallelic transversion sites**.

Because the Rossi panel had been exported through PLINK v1.90, which assigns
provisional reference alleles (`PR` flag) that do not necessarily match the
true reference genome, we normalized the panel against the ARS-UCD1.2_Btau5.0.1Y
reference (ARS-UCD1.2 with the Y chromosome added) using bcftools norm
(v1.21) with `--check-ref ws`. REF/ALT alleles were swapped (and genotypes
flipped) at 247,168 positions (3.48%). The normalized panel
(`rossi_norm.vcf.gz`) was used as the input site list for ANGSD pseudohaploid
calling described below.

## How the SA5 and SA6 VCFs were generated

Raw paired-end Illumina reads were quality-checked with FastQC (v0.12.1),
adapter-trimmed and collapsed with AdapterRemoval2 (v2.3.4; `--trimns
--trimqualities --minquality 25 --minlength 18 --collapse`), and merged across
lanes per sample.

Collapsed reads were aligned to ARS-UCD1.2_Btau5.0.1Y with BWA (v0.7.18) using
the `aln` algorithm (`-l 1024 -n 0.01 -o 2`) followed by `samse`. Alignments
were converted to BAM and unmapped reads removed (samtools v1.18, `-F 4`);
BAM files were sorted and indexed.

Read groups were assigned and PCR duplicates removed with Picard (v3.3.0;
`MarkDuplicates REMOVE_DUPLICATES=true`). Local realignment around indels was
performed with GATK (v3.8; `RealignerTargetCreator` + `IndelRealigner`). Reads
with mapping quality below 25 were removed (`samtools view -q 25`). Post-mortem
damage was assessed with mapDamage (v2.2.3) and coverage with Qualimap (v2.2.1).

Pseudohaploid autosomal genotypes were called with ANGSD (v0.940) on the
final MAPQ ≥ 25 BAM, restricted to autosomes (chromosomes 1–29) and to the
7,108,856-site normalized Rossi panel, using `-doHaploCall 1 -doCounts 1 -GL 1
-minMapQ 30 -minQ 25 -C 50 -uniqueOnly 1 -remove_bads 1 -trim 6`. At each
covered site a single read was randomly sampled, and the sampled base was
returned as the pseudohaploid call. The site list is **not LD-pruned and not
masked for mappability or GC content**; any such downstream filtering is applied
in subsequent analyses, not at the per-sample calling stage represented by these
VCFs.

ANGSD haploid calls were converted to VCF (v4.2) by the following rule: if the
sampled base equalled the REF allele the genotype was recorded as `0/0`; if it
equalled the ALT allele the genotype was recorded as `1/1`; positions where the
sampled base matched neither (third-allele observation, 1.3–1.6% of calls) were
discarded. VCFs were bgzipped and tabix-indexed. Final yields: 129,654 emitted
records for SA5 and 51,131 for SA6.
