Part 3 of 5 — How can a country build its own 1000 Genomes Project? Variant calling, cohort analytics, and data organization

1
import hail as hl
2

3
# Import cohort VCF into Hail MatrixTable
4
mt = hl.import_vcf('s3://genome-project/cohort.vcf.gz')
5

6
# MatrixTable structure:
7
# - Rows: Variants (SNPs, indels, etc.)
8
# - Columns: Samples
9
# - Entries: Genotypes (0/0, 0/1, 1/1, ./.)

1
# Variant-level annotations
2
mt = hl.vep(mt)  # Functional prediction
3
mt = mt.annotate_rows(
4
    gnomad_af = gnomad.rows()[mt.row_key].AF,
5
    clinvar_sig = clinvar.rows()[mt.row_key].significance
6
)
7

8
# Sample-level annotations (sex, ancestry, phenotype)
9
mt = mt.annotate_cols(
10
    sex = sample_meta[mt.col_key].sex,
11
    ancestry = sample_meta[mt.col_key].ancestry,
12
    case_status = sample_meta[mt.col_key].phenotype
13
)

1
# Sample QC
2
sample_qc = hl.sample_qc(vds)
3
mt = mt.annotate_cols(qc = sample_qc[mt.col_key])
4

5
# Filter based on QC metrics
6
mt = mt.filter_cols(
7
    (mt.qc.n_het > 0) &  # Has heterozygous calls
8
    (mt.qc.dp_mean > 30) &  # Average depth >30x
9
    (mt.qc.r_ti_tv > 2.0)  # Ti/Tv ratio normal
10
)
11

12
# Variant QC
13
mt = mt.filter_rows(hl.agg.count_where(mt.GT.is_non_ref()) > 0)

1
# PLINK format (for GWAS)
2
hl.export_plink(mt, 'cohort.bed', 'cohort.bim', 'cohort.fam')
3

4
# BGEN format (compact, efficient)
5
hl.export_bgen(mt, 'cohort.bgen')
6

7
# VDS (Hail native, best for future Hail analyses)
8
mt.write('s3://genome-project/cohort.vds', overwrite=True)

1
# Global allele frequencies
2
mt = mt.annotate_rows(
3
    allele_freq = hl.agg.allele_stats(mt.GT).AF,
4
    ac = hl.agg.count_where(mt.GT.is_non_ref()),
5
    an = hl.agg.count() * 2
6
)
7

8
# Ancestry-stratified frequencies (if multi-population)
9
mt = mt.annotate_rows(
10
    freq_by_ancestry = hl.agg.collect_as_dict(
11
        mt.ancestry,
12
        hl.agg.allele_stats(mt.GT).AF
13
    )
14
)

1
# Hardy-Weinberg equilibrium
2
mt = mt.annotate_rows(
3
    hwe = hl.hardy_weinberg_test(
4
        hl.agg.count_where(mt.GT == hl.call(0, 0)),
5
        hl.agg.count_where(mt.GT == hl.call(0, 1)),
6
        hl.agg.count_where(mt.GT == hl.call(1, 1))
7
    )
8
)
9

10
# Principal Component Analysis (population structure)
11
pc_eigenvalues, pc_scores, pc_loadings = hl.pca(
12
    mt.GT.n_alt_alleles(),
13
    k=20,
14
    compute_loadings=True
15
)

1
# ClinVar pathogenicity
2
clinvar = hl.read_table('gs://hail-datasets/clinvar.ht')
3
mt = mt.annotate_rows(
4
    clinvar_sig = clinvar.rows()[mt.row_key].significance
5
)
6

7
# gnomAD frequencies (population reference)
8
gnomad = hl.read_table('gs://hail-datasets/gnomad.ht')
9
mt = mt.annotate_rows(
10
    gnomad_ac = gnomad.rows()[mt.row_key].AC,
11
    gnomad_af = gnomad.rows()[mt.row_key].AF
12
)
13

14
# CADD scores (variant effect prediction)
15
cadd = hl.read_table('path/to/cadd.ht')
16
mt = mt.annotate_rows(
17
    cadd_score = cadd.rows()[mt.row_key].score
18
)

1
# Flag likely pathogenic variants
2
mt = mt.annotate_rows(
3
    is_likely_pathogenic = (
4
        (mt.clinvar_sig == 'pathogenic') |
5
        ((mt.cadd_score > 30) & (mt.gnomad_af < 0.01))
6
    ),
7
    is_rare = mt.gnomad_af < 0.01,
8
    is_common = mt.gnomad_af > 0.05
9
)

1
# For each rare disease sample, find matching variants in cohort
2
# Example: Exome analysis with population context
3

4
rare_disease_mt = hl.read_matrix_table('rare_disease_samples.vds')
5

6
# Annotate with cohort allele frequencies
7
rare_disease_mt = rare_disease_mt.annotate_rows(
8
    cohort_af = mt.rows()[rare_disease_mt.row_key].allele_freq,
9
    cohort_ac = mt.rows()[rare_disease_mt.row_key].ac
10
)
11

12
# Prioritize rare variants (absent or very rare in cohort)
13
rare_disease_mt = rare_disease_mt.filter_rows(
14
    (rare_disease_mt.cohort_ac < 5) |  # <5 people in cohort
15
    hl.is_missing(rare_disease_mt.cohort_af)  # Not in cohort at all
16
)

1
# Researcher has 1k individuals, 500 with Type 2 Diabetes, 500 without
2
# Goal: Find variants associated with diabetes in your population
3

4
from hail.methods import linear_regression_rows
5

6
# Load PLINK format
7
mt = hl.read_table('cohort.plink')
8

9
# Fit linear regression: genotype ~ diabetes status
10
mt = mt.annotate_rows(
11
    beta = linear_regression_rows(
12
        y=mt.phenotype.diabetes,
13
        x=mt.GT.n_alt_alleles()
14
    )
15
)
16

17
# Find significant associations
18
gwas_hits = mt.filter_rows(mt.p_value < 5e-8)
19

20
# Result: 50-200 genome-wide significant variants unique to your population

1
# Researcher: "I think inactivating mutations in Gene X cause Disease Y"
2
# Goal: Test if individuals with mutations in Gene X have higher disease rates
3

4
# Filter to rare variants in Gene X
5
gene_x_variants = mt.filter_rows(
6
    (mt.gene == "GENE_X") &
7
    (mt.gnomad_af < 0.01)  # Rare in world population
8
)
9

10
# Count carriers in cases vs. controls
11
case_carriers = mt.filter_cols(mt.phenotype.is_case).filter_rows(
12
    mt.GT.is_non_ref()
13
).count_rows() / mt.count_cols()
14

15
# Perform burden test (logistic regression)
16
burden_test = hl.logistic_regression_rows(...)
17

18
# Result: "Individuals with rare variants in Gene X are 5x more likely to have Disease Y"

1
# Researcher: "How does allele frequency of SNP rs123 differ by ancestry?"
2

3
# PCA on your 1k cohort (calculated in 2.3)
4
mt = mt.annotate_cols(
5
    ancestry = predicted_ancestry[mt.col_key]  # PC-based ancestry labels
6
)
7

8
# Allele frequencies stratified by ancestry
9
mt = mt.annotate_rows(
10
    af_by_ancestry = hl.agg.collect_as_dict(
11
        mt.ancestry,
12
        hl.agg.allele_stats(mt.GT).AF
13
    )
14
)
15

16
# Result: "SNP rs123 is 10% in European ancestry, 3% in African ancestry"
17
# Clinical implication: Frequency interpretation depends on patient ancestry

1
# Clinician: "I have a patient with a rare variant; how common is it?"
2
# Query your public allele frequency database (built from 2.4)
3

4
# Web interface: Enter variant rs123456
5
# System returns: "This variant found in 2 of 1000 samples (0.2%)"
6
#                 "gnomAD says 0.05%, so it's enriched in your population"
7
#                 "ClinVar classification: Likely pathogenic"
8
#                 "VEP prediction: This affects splicing"
9

10
# Result: Clinician has population-specific evidence for clinical interpretation

1
# Researcher: "Do people with high polygenic risk scores develop heart disease earlier?"
2
# Requires: 5+ years of follow-up data
3

4
# Build PRS from GWAS hits
5
prs = hl.sum_prs(mt_gwas_hits)
6

7
# Predict time-to-disease (Kaplan-Meier curves)
8
survivors_by_prs = mt.annotate_cols(
9
    prs = prs,
10
    years_to_disease = ehr_data[mt.col_key].years_to_heart_disease
11
)
12

13
# Result: "High PRS group: 30% disease by age 60; Low PRS group: 5%"
14
#         "This enables risk stratification in your healthcare system"

1
# Annual QC check: Re-call a subset of samples with latest pipeline
2
mt_revalidation = hl.import_vcf('revalidated_samples.vcf.gz')
3

4
# Compare allele frequencies to historical records
5
mt_historical = mt.filter_rows(mt_revalidation.row_key)
6

7
# Should be nearly identical (small differences OK, large drift is problem)
8
af_concordance = hl.corr(
9
    mt_historical.allele_freq,
10
    mt_revalidation.allele_freq
11
)
12

13
# Flag if >0.01 drift (suggests pipeline changes or batch effects)
14
if af_concordance < 0.99:
15
    alert("Variant frequency drift detected; investigate pipeline changes")

1
s3://genome-project/archive/gvcfs/
2
├── 2026/01/sample_001.gvcf.gz
3
├── 2026/01/sample_002.gvcf.gz
4
...
5
├── 2026/12/sample_12000.gvcf.gz
6

7
Purpose: Immutable source data (5-10 year retention)
8
Access: Restricted to data stewards
9
Storage class: Cold storage (infrequent access)

1
s3://genome-project/processed/
2
├── cohort_2026_q1.vcf.gz       (Jan-Mar samples)
3
├── cohort_2026_q1.vcf.gz.tbi
4
├── cohort_2026_q2.vcf.gz       (Apr-Jun samples)
5
...
6

7
Purpose: QC baseline, joint-called source
8
Access: Project researchers with approval
9
Retention: 5+ years
10
Storage class: Standard (frequent access)

1
s3://genome-project/analytics/
2
├── cohort_combined.plink.*     (PLINK bed/bim/fam)
3
├── cohort_combined.bgen        (BGEN format)
4
├── cohort_combined.vds/        (Hail VDS)
5
├── population_stats/           (AF, AC, HWE tables)
6
└── public_release_v3.0/        (Anonymized for sharing)
7

8
Purpose: Analysis-ready, indexed, compressed
9
Access: Researchers, public (stratified by tier)
10
Retention: Indefinite for current
11
Storage class: Standard

1
samples.tsv
2
col1: sample_id
3
col2: institution
4
col3: sex
5
col4: ancestry
6
col5: case_status
7
col6: sequencing_date
8
col7: pipeline_version
9
col8: qc_pass (true/false)

1
- Pipeline version (nf-germline-short-read-variant-calling v1.2.3)
2
- Reference genome (GRCh38)
3
- Joint genotyping tool (GLnexus v1.5.2)
4
- Annotation version (VEP version 106)
5
- Date processed (2026-03-15)

1
v1.0 (2026-06): 500 samples, 30M variants
2
v1.1 (2026-09): 750 samples, 35M variants
3
v2.0 (2026-12): 1,000 samples, 40M variants (complete 1k cohort)
4
v2.1 (2027-06): 1,000 samples, 40M variants + updated annotations
5
v3.0 (2028-06): 1,000 samples, 40M variants + multi-omic layers (if extended work pursued)
6

7
Each release:
8
- Frozen cohort VCF (immutable)
9
- Allele frequency tables (stratified by ancestry if multi-population)
10
- QC metrics summary (depth, call rates, Hardy-Weinberg stats)
11
- Known issues & caveats
12
- Citation guidance for researchers using the data