Real data: chr19 association testing with binary, continuous, and count phenotypes

This vignette was run on data that cannot be shared, so its code is shown with the output it produced, but it is not re-run when the package is built. vignette("lantern") runs the same steps on data included in the package.

Overview

We analyse chromosome 19 in 3,313 African American participants of the Jackson Heart Study (JHS), whose local ancestry was inferred as African (AFR) or European (EUR) with RFMix. The phenotypes are simulated, so we do not expect any associations; the aim is to show the workflow on data of realistic size, and how to analyse binary, continuous and count traits. Only the family argument changes between them.

1 Split and write the genotypes

ancestry_split() reads the phased genotypes and the local ancestry and returns one dosage matrix per ancestry. write_ancestry_gds() saves each to a GDS file for GMMAT.

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  = FALSE
)
gds_paths <- write_ancestry_gds(split_real, tempfile(), verbose = FALSE)
dim(split_real$AFR)   # variants x samples
#> [1] 44271  3313

2 Gene groups and kinship

The gene-group file lists, for each gene, the variants to test together (six columns: gene, chromosome, position, reference allele, alternate allele and weight). The full file for chromosome 19 has 1,442 genes. To keep the examples quick we first use the first five genes:

genes_full <- read.table(file.path(RAW, "genes_chr19_smmat.tsv"), 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], ]
table(genes_small$gene)
#> 
#>    A1BG   ABCA7 ABHD17A   ABHD8   ACER1 
#>      41     205      18      30      12

ancestry_smmat() accepts the gene groups either as a data frame, as here, or as a path to a tab-separated file without a header.

The kinship matrix holds the genetic relatedness between every pair of participants, which JHS needs because it includes relatives. We keep the participants present in both the genotypes and the kinship matrix:

kinship <- readRDS(file.path(RAW, "LANTERN_kinship_rds.rds"))
common_ids <- intersect(split_real$sample_ids, rownames(kinship))
length(common_ids)
#> [1] 3313

3 Binary phenotype (family = binomial())

We simulate a binary trait with 30% cases and test the five genes. The trait has no genetic effect, so the p-values should look like random draws between 0 and 1. MAF.range = c(0, 0.05) restricts the test to variants with a minor allele frequency below 5%.

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)
res_binary$results
#>      gene     p_AFR     p_EUR  p_cauchy
#> 1    A1BG 0.1910291 0.8511402 0.5807953
#> 2   ABCA7 0.2370768 0.2414580 0.2392513
#> 3 ABHD17A 0.8004630 0.2707219 0.5784290
#> 4   ABHD8 0.5268295 0.3263750 0.4186764
#> 5   ACER1 0.1019934 0.4238325 0.1752785

4 Continuous phenotype (family = gaussian())

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)
res_cont$results
#>      gene     p_AFR     p_EUR  p_cauchy
#> 1    A1BG 0.3506940 0.4525619 0.3989808
#> 2   ABCA7 0.7445442 0.5025977 0.6443147
#> 3 ABHD17A 0.3875901 0.5085255 0.4461201
#> 4   ABHD8 0.5353446 0.7577901 0.6674988
#> 5   ACER1 0.0807992 0.7691255 0.2014519

5 Count phenotype (family = 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)
res_count$results
#>      gene      p_AFR      p_EUR   p_cauchy
#> 1    A1BG 0.93625150 0.01910404 0.05382056
#> 2   ABCA7 0.44142308 0.62752028 0.53760581
#> 3 ABHD17A 0.36705820 0.62199107 0.49354335
#> 4   ABHD8 0.12410494 0.52148951 0.22338299
#> 5   ACER1 0.02775486 0.92972375 0.08877860

6 Adjusting for covariates

Covariates are added to the formula, and the phenotype table must contain a column for each. Here we add a simulated age and sex to the continuous trait. Genotype principal components are added the same way, for example y ~ age + sex + PC1 + PC2.

set.seed(2)
pheno_cov <- data.frame(id  = common_ids,
                        y   = rnorm(length(common_ids)),
                        age = round(rnorm(length(common_ids), mean = 55, sd = 12)),
                        sex = factor(sample(c("F", "M"), length(common_ids), replace = TRUE)))

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)
res_cov$results
#>      gene      p_AFR     p_EUR   p_cauchy
#> 1    A1BG 0.15641519 0.9971322 0.99416662
#> 2   ABCA7 0.52732801 0.9354582 0.87777687
#> 3 ABHD17A 0.03468871 0.5286584 0.06923531
#> 4   ABHD8 0.76451429 0.6051571 0.69848074
#> 5   ACER1 0.09922975 0.1134745 0.10589289

7 Testing every gene, with ancestry weights

A real analysis tests all genes. Passing ancestry_counts and variant_info from ancestry_split() weights each gene’s Cauchy combination by how many participants of each ancestry carry its variants, so the African-ancestry test, with more data in this cohort, counts for more. This takes a few minutes.

res_full <- ancestry_smmat(gds_paths, pheno_cov, y ~ age + sex, kinship, genes_full,
                           ancestry_counts = split_real$ancestry_counts,
                           variant_info    = split_real$variant_info,
                           MAF.range = c(0, 0.05), verbose = FALSE)
nrow(res_full$results)
#> [1] 1442
head(res_full$results[order(res_full$results$p_cauchy), ])
#>                gene        p_AFR      p_EUR w_AFR w_EUR     p_cauchy
#> 541           JSRP1 8.323902e-05 0.98326189  2283   121 8.767381e-05
#> 301 ENSG00000266953 1.060329e-03 0.29933312  2262   132 1.122045e-03
#> 789          PDCD2L 1.294066e-03 0.29521556  2262   132 1.369337e-03
#> 911            RDH8 2.275999e-03 0.10761493  2268   122 2.395803e-03
#> 29              AMH 3.948893e-03 0.02433174  2283   121 4.122772e-03
#> 393         FAM187B 4.087365e-03 0.24260357  2270   120 4.300353e-03

The trait is simulated without genetic effects, so the smallest of these p-values are what we expect by chance among 1,442 genes.

8 Other SMMAT tests and the raw SMMAT output

By default SMMAT runs its efficient hybrid test (tests = "E"), and ancestry_smmat() combines its p-values (column E.pval). Any other GMMAT::SMMAT() argument can be passed through. For example, to combine the SKAT-O p-values instead, request the test with tests and name the p-value column with p_col:

res_o <- ancestry_smmat(gds_paths, pheno_cov, y ~ age + sex, kinship, genes_small,
                        tests = "O", p_col = "O.pval",
                        MAF.range = c(0, 0.05), verbose = FALSE)
res_o$results
#>      gene      p_AFR     p_EUR  p_cauchy
#> 1    A1BG 0.10482040 1.0000000 1.0000000
#> 2   ABCA7 0.55127668 0.8140770 0.7219490
#> 3 ABHD17A 0.09469666 0.4024312 0.1622402
#> 4   ABHD8 0.86583872 0.6980776 0.8102436
#> 5   ACER1 0.31388522 0.1143865 0.1724220

The unsummarised output of each ancestry’s SMMAT run is kept in smmat_results, for example the number of variants tested in each gene and the burden test statistics:

head(res_o$smmat_results$AFR)
#>     group n.variants miss.min miss.mean miss.max     freq.min    freq.mean
#> 1    A1BG         41        0         0        0 0.000000e+00 0.0004614505
#> 2   ABCA7        205        0         0        0 0.000000e+00 0.0006208557
#> 3 ABHD17A         18        0         0        0 0.000000e+00 0.0002688709
#> 4   ABHD8         30        0         0        0 0.000000e+00 0.0004144215
#> 5   ACER1         12        0         0        0 6.036825e-05 0.0005768488
#>      freq.max    B.score     B.var     B.pval     S.pval     O.pval     O.minp
#> 1 0.003760942  161.15499  70793.09 0.54472257 0.05548930 0.10482040 0.05548930
#> 2 0.008212858  661.37964 484767.25 0.34215607 0.57782108 0.55127668 0.34215607
#> 3 0.001266088 -228.15606  15523.39 0.06706888 0.08400736 0.09469666 0.05159878
#> 4 0.002712375  -87.38494  45539.21 0.68217974 0.70218112 0.86583872 0.68217974
#> 5 0.003697555 -185.41816  21148.99 0.20231204 0.38163155 0.31388522 0.20231204
#>   O.minp.rho
#> 1       0.00
#> 2       1.00
#> 3       0.25
#> 4       1.00
#> 5       1.00

9 More than two populations

This cohort only has African and European local ancestry. With an MSP file that declares more populations, ancestry_split() returns one matrix per population and write_ancestry_gds() one GDS file per population, and ancestry_smmat() is called in the same way. For example, with RFMix output for African, European and Native American ancestry (the file names are placeholders):

split3 <- ancestry_split("cohort.vcf.gz", "cohort.msp.tsv", chrom = "chr19")
gds3   <- write_ancestry_gds(split3, "gds")   # gds/AFR.gds, gds/EUR.gds, gds/NAT.gds
res3   <- ancestry_smmat(gds3, pheno, y ~ age + sex, kinship, genes,
                         ancestry_counts = split3$ancestry_counts,
                         variant_info    = split3$variant_info)

vignette("toy-vcf-msp-example") builds a four-population example.

10 Gene-group files for your own data

A gene-group file is specific to your variants: it pairs each variant in your VCF with the gene whose boundaries contain it. The gene coordinates can be downloaded, for example from Ensembl BioMart, GENCODE or the UCSC Table Browser, and intersected with your variant positions, for instance with bedtools intersect or in R. The file used above was made this way. A weight of 1 for every variant gives an unweighted test; SMMAT can also weight variants by allele frequency (see MAF.weights.beta in ?GMMAT::SMMAT).