vignettes/real-data-chr19.Rmd
real-data-chr19.RmdThis vignette runs the complete Step 1 → 2 → 3 pipeline on real
genotype data from the Jackson Heart Study (JHS) cohort — chromosome 19,
3 313 samples, local ancestry inferred as African (AFR) / European (EUR)
— and demonstrates ancestry_smmat() with all three response
types GMMAT supports: binary (logistic), continuous (linear), and count
(Poisson). Only the family argument (and the simulated
phenotype) changes between them; everything else about the call is
identical, which is the point — swap those two inputs and the same code
handles your own binary, continuous, or count phenotype.
vignette("split-intuition").ancestry_split() inputs from your own VCF
+ MSP files (including K > 2 populations), see
vignette("toy-vcf-msp-example") — this JHS cohort only has
AFR/EUR calls, so Steps 1–2 below are shown for two populations, but
nothing changes for K > 2 beyond the MSP header.
# Step 1: split. Step 2: write one GDS per population.
split_real <- ancestry_split(
vcf_path = file.path(RAW, "LANTERN.phased.bcf"),
msp_path = file.path(RAW, "LANTERN_chr19.msp.tsv.gz"),
mode = "dosage",
chrom = "chr19",
verbose = TRUE
)
#> === LANTERN Ancestry Split (Step 1, mode = dosage) ===
#> Step 1: Parsing RFMix MSP file...
#> Parsing MSP file: ../../raw/LANTERN_chr19.msp.tsv.gz
#> Found 3416 samples, 2 populations
#> Pop codes: AFR = 0, EUR = 1
#> MSP samples: 3416 Tracts: 653
#>
#> Step 2: Querying VCF sample IDs...
#> VCF samples: 3314
#>
#> Step 3: Intersecting samples...
#> Common: 3313
#>
#> Step 4: Parsing VCF genotypes...
#> 44275 variant lines read
#>
#> Step 5: Mapping variants to ancestry tracts...
#> Filtered to chrom chr19: 44275 variants
#> Variants with no tract: 2 (dropped)
#>
#> Step 6: Broadcasting ancestry to variants...
#>
#> Step 7: Splitting genotypes by ancestry (dosage)...
#>
#> Step 8: Filtering monomorphic variants...
#> Removed 2 monomorphic variants
#>
#> === Step 1 Complete: 44271 variants, 3313 samples, 2 populations ===
gds_paths <- write_ancestry_gds(split_real, tempfile(), verbose = TRUE)
#> === LANTERN Write Ancestry GDS (Step 2) ===
#> Writing GDS for AFR (1/2)...
#> Writing GDS for EUR (2/2)...
#>
#> === Step 2 Complete: 2 GDS file(s) written to /tmp/Rtmp7oMjf3/file699222c4d9016 ===
ids <- split_real$sample_ids
cat("Samples:", length(ids), "\n")
#> Samples: 3313The full chr19 gene-group file covers 1 442 genes; to keep this
vignette fast we use a 5-gene subset here (the pattern is identical for
the full file — just point gene_group_file at it and expect
several minutes per SMMAT() call instead of seconds).
genes_full <- read.table(file.path(RAW, "genes_chr19_smmat.tsv"), header = FALSE,
sep = "\t",
col.names = c("gene", "chrom", "pos", "ref", "alt", "weight"))
genes_small <- genes_full[genes_full$gene %in% unique(genes_full$gene)[1:5], ]
kinship <- readRDS(file.path(RAW, "LANTERN_kinship_rds.rds"))
common_ids <- intersect(ids, rownames(kinship))
cat("Genes:", length(unique(genes_small$gene)), " Common samples:", length(common_ids), "\n")
#> Genes: 5 Common samples: 3313genes_small is passed straight to
ancestry_smmat()’s gene_group_file argument
below as a data.frame — no
tempfile()/write.table()/unlink()
needed; ancestry_smmat() manages the temporary file GMMAT
requires internally and cleans it up automatically. A plain file path
still works too, if that’s more convenient for your own pipeline.
family = binomial())
set.seed(1)
pheno_binary <- data.frame(id = common_ids, y = rbinom(length(common_ids), 1, 0.3))
res_binary <- ancestry_smmat(gds_paths, pheno_binary, y ~ 1, kinship, genes_small,
family = binomial(link = "logit"),
MAF.range = c(0, 0.05), verbose = FALSE)
print(res_binary$results)
#> gene p_AFR p_EUR p_cauchy
#> 1 A1BG 0.1492616 0.7938859 0.3997315
#> 2 ABCA7 0.2031324 0.1900808 0.1964182
#> 3 ABHD17A 0.8941616 0.2588411 0.7459638
#> 4 ABHD8 0.5078140 0.3048892 0.3958519
#> 5 ACER1 0.2196545 0.3375616 0.2692847family = gaussian())
Only the family and the simulated y
change:
pheno_cont <- data.frame(id = common_ids, y = rnorm(length(common_ids)))
res_cont <- ancestry_smmat(gds_paths, pheno_cont, y ~ 1, kinship, genes_small,
family = gaussian(link = "identity"),
MAF.range = c(0, 0.05), verbose = FALSE)
print(res_cont$results)
#> gene p_AFR p_EUR p_cauchy
#> 1 A1BG 0.32271246 0.5564923 0.4305589
#> 2 ABCA7 0.71670271 0.3862533 0.5684131
#> 3 ABHD17A 0.31455720 0.6888582 0.5024665
#> 4 ABHD8 0.56071838 0.8476944 0.7593479
#> 5 ACER1 0.09634068 0.7767378 0.2485254family = poisson())
pheno_count <- data.frame(id = common_ids, y = rpois(length(common_ids), lambda = 3))
res_count <- ancestry_smmat(gds_paths, pheno_count, y ~ 1, kinship, genes_small,
family = poisson(link = "log"),
MAF.range = c(0, 0.05), verbose = FALSE)
print(res_count$results)
#> gene p_AFR p_EUR p_cauchy
#> 1 A1BG 0.89654667 0.02082793 0.0513362
#> 2 ABCA7 0.45190129 0.66226067 0.5638687
#> 3 ABHD17A 0.29662855 0.75212153 0.5428024
#> 4 ABHD8 0.11031509 0.65554127 0.2321887
#> 5 ACER1 0.04307048 0.81238745 0.1048616Covariates (age, sex, genotype PCs, batch, …) are just additional
terms in formula, as long as the matching columns exist in
pheno — nothing else about the call changes. Here we add a
simulated age and sex to the
continuous-phenotype example from Section 3:
set.seed(2)
age <- round(rnorm(length(common_ids), mean = 55, sd = 12))
sex <- factor(sample(c("F", "M"), length(common_ids), replace = TRUE))
pheno_cov <- data.frame(id = common_ids, y = rnorm(length(common_ids)),
age = age, sex = sex)
res_cov <- ancestry_smmat(gds_paths, pheno_cov, y ~ age + sex, kinship, genes_small,
family = gaussian(link = "identity"),
MAF.range = c(0, 0.05), verbose = FALSE)
print(res_cov$results)
#> gene p_AFR p_EUR p_cauchy
#> 1 A1BG 0.78973635 0.63646754 0.72827528
#> 2 ABCA7 0.61704448 0.70291389 0.66320132
#> 3 ABHD17A 0.05755877 0.08958273 0.07014252
#> 4 ABHD8 0.61776792 0.42714304 0.52459912
#> 5 ACER1 0.10424888 0.90977046 0.57638538Genotype PCs from your own PCA (e.g. PC1, …,
PC10) work the same way: add them as columns to
pheno and to formula
(y ~ age + sex + PC1 + PC2 + ...).
ancestry_smmat()’s gene_group_file (a
data.frame or a path — both resolve to the same six columns, no header:
gene chrom pos ref alt weight) is fundamentally
your own called variants intersected with gene
boundaries — it’s specific to your VCF, not something you
download ready-made for an arbitrary variant set.
genes_full, used above, was built exactly this way: gene
coordinates for chr19 intersected against the JHS variant positions.
What you can download are the gene/exon coordinate
annotations to intersect against your own variants (e.g. with
bedtools intersect or a GTF-parsing script), for human or
any other species:
refGene/ensGene databases), its per-variant
gene assignments can be reshaped directly into the
gene_group_file format.Once you have gene boundaries in hand, build the file per chromosome:
for each variant in your VCF, look up which gene’s boundary it falls in
(chromosome + position), write one row per (gene, variant) pair, and set
weight (commonly 1 for an unweighted burden
test, or a MAF-based weight — see MAF.weights.beta in
?ancestry_smmat).
gds_paths just needs one named entry per
population (see vignette("toy-vcf-msp-example"), Part B,
for building the underlying
ancestry_split()/write_ancestry_gds() inputs
at K > 2).genes_full instead of the 5-gene subset as
gene_group_file, and pass
ancestry_counts = split_real$ancestry_counts /
variant_info = split_real$variant_info (from Step 1) to
ancestry_smmat() to weight the Cauchy combination by each
gene’s local ancestry proportions instead of equal weights.GMMAT::SMMAT() argument is
available: MAF.range, miss.cutoff,
method, and ncores are named parameters here
with the same meaning as in a direct SMMAT() call; anything
else SMMAT() accepts — tests,
rho, MAF.weights.beta,
missing.method, use.minor.allele,
auto.flip, Garbage.Collection,
group.file.sep, meta.file.prefix — forwards
through ... unchanged. See ?ancestry_smmat and
?GMMAT::SMMAT.$smmat_results, a named list of the untouched
GMMAT::SMMAT() data.frames
(e.g. res_binary$smmat_results$AFR) if you need columns
beyond the summarized p_<POP>/p_cauchy
in $results — burden-only p-values, variant counts, allele
frequencies, etc.