1000 Genomes FUT2 slices¶

Committed three-SNP extracts from 1000 Genomes Project panels (EAS, EUR, AFR) demonstrate ancestry-specific secretor calls and coordinate lookup when IDs are CHR:POS:REF:ALT rather than rsIDs.

Source panels live under $GWASLAB_DIR (default ~/.gwaslab). genoSecretor only fetches the three FUT2 tags, so it can also run on the full VCF in place.

Prerequisites¶

# from repository root
pip install -e .
# regenerate slices if you have the local 1KG panels
python examples/extract_1kg.py          # or: python examples/extract_1kg.py /path/to/.gwaslab
./examples/run_1kg.sh

bcftools and tabix are required only to re-extract. The committed files under examples/data/1kg/ are enough to run ./examples/run_1kg.sh. See examples/data/1kg/README.txt for panel provenance.

Pairwise haplotype r² among the three tags: Marker LD. Regenerate with ./examples/run_1kg_ld.sh.

The unphased --assume-trans rule vs phased cis/trans is in 1000 Genomes --assume-trans.

Layout¶

examples/data/1kg/
  eas_hg19.vcf.gz   # 1KG Phase 3 EAS, GRCh37, contig 19
  eur_hg19.vcf.gz
  afr_hg19.vcf.gz
  eas_hg38.vcf.gz   # 1KG 30x EAS, GRCh38, contig chr19
  eur_hg38.vcf.gz
  afr_hg38.vcf.gz

Each file holds rs1047781, rs601338, and rs602662 only.

Phase auto (EAS hg38)¶

The runner prefers, in order: GENOSECRETOR_1KG_EAS_VCF, ~/.gwaslab/EAS.ALL.split_norm_af.1kg_30x.hg38.vcf.gz (full indexed panel), then the bundled extract.

#!/usr/bin/env bash
# Run genoSecretor on 1000 Genomes 30x EAS (hg38): full local panel if present, else bundled extract.
# Uses --phase auto (detects phased GT).
set -euo pipefail
ROOT="$(cd "$(dirname "$0")" && pwd)"
ANC=EAS
# shellcheck source=_1kg_input.sh
source "$ROOT/_1kg_input.sh"

mkdir -p "$ROOT/out"
echo "Using: $VCF"
genosecretor "$VCF" \
  --build hg38 \
  --method auto \
  --phase auto \
  -o "$ROOT/out/1kg_eas"
echo "Open: $ROOT/out/1kg_eas.report.html"

./examples/run_1kg_eas.sh

EUR and AFR: ./examples/run_1kg_eur.sh, ./examples/run_1kg_afr.sh.

Allele frequencies¶

ALT allele frequencies in the source panels (se385 = T at rs1047781; se428 = A at rs601338; proxy = A at rs602662):

se385 is common in East Asians and absent in these EUR/AFR extracts. se428 (and its LD proxy) is common in EUR/AFR and nearly absent in EAS.

Expected calls (--method auto)¶

auto selects two-marker because both functional SNPs are present.

hg38 counts match hg19 for the same ancestry (build and chr19 vs 19 smoke test).

--phenotype-coding gwas recodes the 97 EAS weak-secretor calls as non-secretor (407 secretor, 0 weak-secretor, 97 non-secretor). EUR and AFR counts are unchanged. Committed HTML reports use biochemical coding.

Example report¶

Snapshot of PREFIX.report.html from the bundled hg38 extracts (--method auto --phase auto). Research demo, not clinical typing. The report gallery has the same snapshots plus the synthetic cohort.

Forced unphased (EAS)¶

--phase unphased uses the same VCF but ignores | phase. Default --assume-trans assigns the one EAS double-het to Se^w/se, so phenotype counts match the phased run (407 secretor, 97 weak-secretor). All 504 calls are phase_mode=unphased (confidence medium). Pass --no-assume-trans to leave that sample ambiguous:

#!/usr/bin/env bash
# Same 1000 Genomes 30x EAS VCF as run_1kg_eas.sh, but ignore phase (--phase unphased).
# The file is still phased (0|1); this run treats it as unphased genotypes.
set -euo pipefail
ROOT="$(cd "$(dirname "$0")" && pwd)"
ANC=EAS
# shellcheck source=_1kg_input.sh
source "$ROOT/_1kg_input.sh"

mkdir -p "$ROOT/out"
echo "Using: $VCF (forced unphased)"
genosecretor "$VCF" \
  --build hg38 \
  --method auto \
  --phase unphased \
  -o "$ROOT/out/1kg_eas_unphased"
echo "Open: $ROOT/out/1kg_eas_unphased.report.html"

./examples/run_1kg_eas_unphased.sh

Counts are a software check, not a population frequency claim. Sample-level concordance of that double-het with phased Se^w/se, and the same check in AFR/AMR/EUR, is in 1000 Genomes --assume-trans.

Marker notation¶

IDs are 19:pos:REF:ALT (hg19) or chr19:pos:REF:ALT (hg38 30x), not rsIDs. genoSecretor matches markers by coordinates and contig aliases.

In-place full panel¶

If the GWASLab 1KG VCFs are installed, genoSecretor tabix-fetches the three tags without extracting:

genosecretor ~/.gwaslab/EAS.ALL.split_norm_af.1kg_30x.hg38.vcf.gz --build hg38 -o examples/out/1kg_eas_full

The local EAS hg19 PLINK2 files (.pgen) can be used the same way. Full-panel outputs stay under examples/out/ and are not committed.

Citation¶

1000 Genomes Project Consortium. A global reference for human genetic variation. Nature. 2015;526:68–74. DOI: 10.1038/nature15393.

High-coverage 30x VCFs: Byrska-Bishop et al. High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios. Cell. 2022;185:3426–3440.e19. DOI: 10.1016/j.cell.2022.08.004.