--- title: "Real data: chr19 association testing with binary, continuous, and count phenotypes" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Real data: chr19 association testing with binary, continuous, and count phenotypes} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ## Overview This 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. * For matrix-level intuition on the two splitting methods, see `vignette("split-intuition")`. * For building `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. ``` r # 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: 3313 ``` ## 1 Gene-group file and kinship The 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). ``` r 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: 3313 ``` `genes_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. ## 2 Binary phenotype (`family = binomial()`) ``` r 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.2692847 ``` ## 3 Continuous phenotype (`family = gaussian()`) Only the `family` and the simulated `y` change: ``` r 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.2485254 ``` ## 4 Count phenotype (`family = poisson()`) ``` r 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.1048616 ``` ## 5 Adjusting for covariates Covariates (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: ``` r 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.57638538 ``` Genotype 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 + ...`). ## 6 Where to get gene-group files for other species/builds `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: * [Ensembl BioMart](https://www.ensembl.org/biomart/martview) — gene/exon coordinates for hundreds of species (GRCh37/GRCh38 human builds included); the classic way to export a custom gene coordinate table. * [GENCODE](https://www.gencodegenes.org/) — comprehensive human and mouse gene annotation (GTF/GFF3), the reference set most rare-variant pipelines build their gene definitions from. * [UCSC Table Browser](https://genome.ucsc.edu/cgi-bin/hgTables) — RefSeq/ Ensembl/GENCODE gene tables for many species and genome builds, with a web UI for filtering/exporting BED files directly. * [NCBI RefSeq](https://www.ncbi.nlm.nih.gov/refseq/) — the underlying gene/transcript annotation source Ensembl/UCSC also redistribute. * [ANNOVAR](https://annovar.openbioinformatics.org/en/latest/) — if you already annotate variants with ANNOVAR (`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`). ## 7 Other ways to adapt this to your own data * **More than two populations**: nothing above changes structurally — `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). * **Full gene set / real analysis**: pass `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. * **Every `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`. * **Raw per-population results**: every call above also returns `$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_`/`p_cauchy` in `$results` — burden-only p-values, variant counts, allele frequencies, etc.