Harmonization

GWASLab provides reference-dependent harmonization functions.

Harmonization Workflow

The harmonization process follows a sequential workflow to align and standardize summary statistics with reference data:

Load Sumstats
     │
     ▼
Basic Check (optional)
     │
     ▼
Check Reference Alignment
     │
     ├─→ NEA ≠ REF ─→ Flip Allele Stats ─┐
     │                                    │
     └─→ NEA == REF ──────────────────────┘
                         │
                         ▼
                    Assign rsID
                         │
                         ▼
                    Infer Strand
                         │
                         ├─→ Strand mismatch ─→ Flip Allele Stats ─┐
                         │                                         │
                         └─→ Strand correct ───────────────────────┘
                                             │
                                             ▼
                                    Harmonized Sumstats

Workflow Steps:

1. Basic Check (optional, basic_check=True): Fix IDs, chromosomes, positions, and alleles
2. Check Reference Alignment (check_ref()): Verify NEA matches REF in reference genome
3. Flip Allele Stats (if needed): Flip alleles and allele-specific statistics when NEA ≠ REF
4. Assign rsID (assign_rsid()): Annotate variants with rsIDs from reference
5. Infer Strand (infer_strand()): Determine strand orientation for palindromic SNPs and indels
6. Flip Allele Stats (if needed): Flip again based on strand inference results

mysumstats.harmonize(
    ref_seq="reference.fasta",
    ref_rsid_vcf="dbsnp.vcf.gz",
    ref_infer="1kg.vcf.gz",
    ref_alt_freq="AF"
)

Methods summary

Align NEA with REF in the reference genome

.check_ref(): Check if NEA is aligned with the reference sequence. After checking, the tracking status code will be changed accordingly.

mysumstats.check_ref(ref_path="ref_genome.fa")
mysumstats.flip_allele_stats()

Assign rsID according to CHR, POS, REF/ALT

.assign_rsid() : Annotated variants with rsID using a reference tsv file (1KG variants) and reference vcf file (tabix indexed, entire dbSNP).

See https://cloufield.github.io/gwaslab/AssignrsID/

mysumstats.assign_rsid(ref_rsid_tsv = gl.get_path("1kg_dbsnp151_hg19_auto"), 
                       ref_rsid_vcf = "/home/yunye/mydata/d_disk/dbsnp/GCF_000001405.25.vcf.gz",
                       threads=1)

• For TSV file, variants will be matched using SNPID (CHR:POS:NEA:EA) for quick assigning.
• For VCF file, GWASLab will first extract all variants in the reference file with matching CHR and POS. And then compare EA/NEA in sumstats with REF/ALT in reference vcf. When matching, it will annotate the variant in sumstats with the matching rsID in reference vcf.

Check palindromic SNPs or indistinguishable Indels

.infer_strand() and .infer_strand2():

• Infer the strand for palindromic SNPs (A/T, T/A, G/C, C/G) with MAF ≤ maf_threshold.
• Check the alignment status of indels with the REF allele in a reference VCF/BCF file and verify if the allele frequencies are consistent (DAF ≤ daf_tolerance).
• Update STATUS code 7th digit to reflect strand inference results (0=non-palindromic, 1=forward, 2=reverse, 3=indel forward, 4=ambiguous indel, 6=indel reverse, 7=MAF too high, 8=not found, 9=unchecked).

# Old method: per-variant VCF queries (good for small datasets)
mysumstats.infer_strand(ref_infer=gl.get_path("1kg_eas_hg19"),
                        ref_alt_freq="AF",
                        maf_threshold=0.40,
                        daf_tolerance=0.20,
                        threads=4)
mysumstats.flip_allele_stats()

# New optimized method: bulk lookup table (faster for large datasets)
mysumstats.infer_strand2(vcf_path=gl.get_path("1kg_eas_hg19"),
                         maf_threshold=0.40,
                         ref_maf_threshold=0.40,
                         daf_tolerance=0.20,
                         threads=6)
mysumstats.flip_allele_stats()

# Using cache to speed up old method
mysumstats.infer_strand(ref_infer=gl.get_path("1kg_eas_hg19"),
                        cache_options={"use_cache":True})
mysumstats.flip_allele_stats()

Check the difference in allele frequency

.check_af() : Check the allele frequency discrepancy with a reference VCF/BCF file. This function calculates the difference (DAF) between the effect allele frequency (EAF) in your sumstats and the alternative allele frequency (ALT_AF) in the reference VCF.

Purpose:

• Quality control: Identify variants with large differences in allele frequency between sumstats and reference
• Validation: Verify that EAF values in sumstats are consistent with reference population frequencies
• Detect potential issues: Large |DAF| values (> 0.2) may indicate allele mismatches, strand flips, or data quality issues

Important: Please make sure your sumstats are already harmonized, and the variants in reference VCF are also aligned. GWASLab will retrieve information only for matched variants (CHR, POS, EA-ALT, and NEA-REF).

mysumstats.check_af(ref_infer=gl.get_path("1kg_eas_hg19"), 
                    ref_alt_freq="AF",
                    threads=2)

Terminology:

• DAF: Difference in Allele Frequency = EAF (sumstats) - ALT_AF (reference VCF)
• EAF: Effect allele frequency (from your sumstats)
• RAF: Reference ALT allele frequency (from reference VCF/BCF file)

Interpretation:

• Positive DAF: EAF in sumstats is higher than reference
• Negative DAF: EAF in sumstats is lower than reference
• Large |DAF| values (> 0.2) may indicate:
• Population differences (expected for population-specific variants)
• Allele mismatches (check EA/NEA alignment)
• Strand flips (use infer_strand() to resolve)
• Data quality issues (verify EAF calculation in sumstats)

You may want to check the allele frequency discrepancy with a reference VCF. Just specify the path and the right allele frequency for your target ancestry in INFO field.

Allele frequency correlation plot

GWASLab will simply calculate DAF = EAF (sumstats) - frequency in VCF file, and store the results in DAF column.

DAF can then be used for plotting (.plot_daf()) or filter variants.

mysumstats.plot_daf(threshold=0.12)

Flip allele-specific statistics

.flip_allele_stats() : Flip allele-specific statistics to harmonize the variants based on the tracking status code in STATUS.

mysumstats.check_ref(ref_path="ref_genome.fa")
mysumstats.flip_allele_stats()

mysumstats.infer_strand()
mysumstats.flip_allele_stats()

Sweep Mode for Large Datasets

Sweep mode is an optimized processing mode designed for large datasets that significantly improves performance when harmonizing summary statistics with reference VCF/BCF files by reducing I/O operations. Sweep mode uses bcftools for efficient VCF/BCF processing.

bcftools --version

bcftools 1.19
Using htslib 1.19
Copyright (C) 2023 Genome Research Ltd.
License Expat: The MIT/Expat license
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.

How Sweep Mode Works

When sweep_mode=True, GWASLab uses a different processing strategy that creates lookup tables instead of querying the reference file per-variant:

• Normal mode (sweep_mode=False):
• For strand inference (ref_infer): Uses _parallelize_infer_strand which queries the reference VCF/BCF file using tabix for each variant (or chunk of variants). Each query requires a separate I/O operation to the reference file.

• For rsID assignment (ref_rsid_vcf): Uses _parallelize_assign_rsid which also performs per-variant or per-chunk tabix queries against the reference VCF/BCF file.

• Sweep mode (sweep_mode=True):

• For strand inference (ref_infer): Uses _infer_strand_with_annotation which:
  1. First extracts all needed variants from the reference VCF/BCF in one pass using bcftools (bcftools view, bcftools query) to create a lookup table (TSV)
  2. Then annotates all variants in your sumstats using this lookup table (fast in-memory operations)
  3. Finally infers strand orientation using the annotated allele frequencies
• For rsID assignment (ref_rsid_vcf): Uses _assign_rsid which:
  1. First extracts a lookup table from the reference VCF/BCF in one pass using bcftools (bcftools view, bcftools query, bcftools index) via _extract_lookup_table_from_vcf_bcf
  2. Then assigns rsIDs to all variants using _assign_from_lookup (fast in-memory operations)

Key Advantages of Sweep Mode

• Reduced I/O operations: Instead of many tabix queries (one per variant/chunk), sweep mode makes a single pass through the reference file to extract all needed data
• Better performance for large datasets: The lookup table approach scales better when processing millions of variants
• Reusable lookup tables: The extracted lookup tables can be saved and reused in subsequent runs (if lookup_path is specified)
• Lower overhead: Fewer file system operations and less network I/O (if using remote files)

When to Use Sweep Mode

Use sweep mode when: - Processing large datasets (millions of variants) - Using large reference VCF/BCF files (e.g., full dbSNP, 1000 Genomes) - You want to minimize I/O operations and improve throughput - You plan to reuse lookup tables for multiple harmonization runs - Processing variants from remote/network-mounted reference files

Use normal mode when: - Processing small to medium datasets (< 1 million variants) - You need maximum parallelization across CPU cores (normal mode parallelizes tabix queries) - The reference file is small or already indexed in a way that makes per-variant queries efficient - You prefer lower memory usage (sweep mode loads lookup tables into memory)

Usage Example

# Harmonization with sweep mode for large dataset
mysumstats.harmonize(
    ref_seq="/path/to/reference.fasta",
    ref_rsid_vcf="/path/to/dbsnp.vcf.gz",
    ref_infer="/path/to/reference.vcf.gz",
    ref_alt_freq="AF",
    sweep_mode=True,  # Enable sweep mode
    threads=8
)

Assign CHR and POS according to rsID and reference data

mysumstats.rsid_to_chrpos()  

mysumstats.rsid_to_chrpos2()