Practical 7 - Rare Variant Analysis

Last updated: 2024-06-14

Checks: 6 1

Knit directory: SISG2024_Association_Mapping/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20230530)

The command set.seed(20230530) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute	relative
/Users/joelle.mbatchou/SISG/2024/SISG2024_Association_Mapping/data/	data

Repository version: c788d93

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version c788d93. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    analysis/.DS_Store
    Ignored:    data/run_regenie.r
    Ignored:    data/sim_rels_geno.bed
    Ignored:    exe/
    Ignored:    gwas_plink.log
    Ignored:    gwas_regenie.log
    Ignored:    lectures/
    Ignored:    mk_website.R
    Ignored:    notes.txt
    Ignored:    step2_gwas_regenie.log
    Ignored:    tmp/

Untracked files:
    Untracked:  .Rhistory
    Untracked:  .notes.txt.swp
    Untracked:  analysis/SISGM15_prac4Solution.Rmd
    Untracked:  analysis/SISGM15_prac5Solution.Rmd
    Untracked:  analysis/SISGM15_prac6Solution.Rmd
    Untracked:  analysis/SISGM15_prac9Solution.Rmd
    Untracked:  analysis/Session02_practical_Key_cache/
    Untracked:  sim.config

Unstaged changes:
    Modified:   .gitignore
    Modified:   analysis/Session08_practical_Key.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/Session07_practical.Rmd) and HTML (docs/Session07_practical.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
html	e8bf041	Joelle Mbatchou	2024-06-13	Build site.
Rmd	a497e68	Joelle Mbatchou	2024-06-13	update practical
html	a497e68	Joelle Mbatchou	2024-06-13	update practical
html	3499240	Joelle Mbatchou	2024-06-05	update session pages
html	de4ef79	Joelle Mbatchou	2024-06-05	add docs dir

Before you begin:

Make sure that R is installed on your computer
For this lab, we will use the following R libraries:

library(data.table)
library(dplyr)
library(BEDMatrix)
library(SKAT)
library(ACAT)
library(ggplot2)

Introduction

We will look into a dataset collected on a quantitative phenotype which was first analyzed through GWAS and a signal was detected in a region on chromosome 1. Let’s determine whether the signal is present when we focus on rare variation at the locus. In our analyses, we will define rare variants as those with \(MAF \leq 5\%\).

The file “rv_pheno.txt”” contains the phenotype measurements for a set of individuals and the file “rv_geno_chr1.bed” is a binary file in PLINK BED format with accompanying BIM and FAM files which contains the genotype data.

Data preparation

Let’s first load the files into the R session. We first need to define the path to the directory containing the phenotype and genotype files (change the path to the files location on your machine).

files_dir <- "/SISGM19/data/"

Also specify the paths to the PLINK2 binary:

plink2_binary <- "/SISGM19/bin/plink2"

We can now read the phenotype file:

pheno_file <- fread(sprintf("%s/rv_pheno.txt", files_dir), header = TRUE) 
head(pheno_file, 3)

    FID  IID       Pheno
1: 5257 5257  0.73085382
2: 4686 4686  0.38374848
3: 5818 5818 -0.03473697

Exercises

Using PLINK, extract rare variants in a new PLINK BED file. We use options --max-maf to select rare variants and --maj-ref force so that the minor allele is the effect allele.

# first fill in the thresholds to use for each filter
filter_maf = 

cmd <- sprintf('%s --bfile "%s/rv_geno_chr1" --max-maf %g --maj-ref force --make-bed --out chr1_region_rv', plink2_binary, files_dir, filter_maf)
system(cmd)

This generates a new set of PLINK BED genotype files containing only the variants that passed the MAF filter chr1_region_rv.{bed,bim,fam}.

Import the data files in R.

Read the genotype data using R function BEDMatrix(). We use option simple_names = TRUE to easily filter by sample IDs later.

G <- BEDMatrix("chr1_region_rv", simple_names = TRUE)

Keep only samples who are present both in the genotype as well as phenotype data and who don’t have missing values for the phenotype.

# identify ID of samples with non-missing phenotypes
ids.keep <- with(pheno_file, IID[ !is.na(Pheno) ])
# subset the genotype & phenotype data based on IDs
G <- G[match(ids.keep, rownames(G)), ]
pheno_file <- pheno_file[match(ids.keep, pheno_file$IID), ]

Examine the genotype data:

Compute the minor allele frequency (MAF) for each SNP and plot histogram. We use the colMeans() function to G and specify the argument na.rm = TRUE in case missing genotypes are present.

# Recall the MAF formula: maf(g) = sum(g) / (2*N) = mean(g) / 2
maf <- colMeans(G, na.rm=TRUE)/2
# we can use the 'hist' function in R to plot histograms
hist(maf)

Run the single variant association tests in PLINK (only for the extracted variants).

cmd <- sprintf("%s --bfile chr1_region_rv --pheno %s/rv_pheno.txt --pheno-name Pheno --glm allow-no-covars  --out test_plink", plink2_binary, files_dir)
system(cmd)
sv_pvals <- fread("test_plink.Pheno.glm.linear")
str(sv_pvals) # what variables are used to store p-values and effect sizes?

What would be your significance threshold after applying Bonferroni correction for the multiple tests (assume the nominal significance level is 0.05)? Is anything significant after this correction?

# determine the appropriate significance threshold
bonf_threshold <- 
sv_pvals[ sv_pvals$P < bonf_threshold, ]

Make a volcano plot (i.e. plot log10 p-values against the effect sizes). Which of the Burden/SKAT/ACAT tests do you expect will give us most power?

# volcano plot hint: to make a scatterplot of variable var2 against var1 stored in the data frame, you can use:
with(sv_pvals, plot(x=var1, y=var2)) # change 'var1' and 'var2' to the right names in sv_pvals

We will first compare three collapsing/burden approaches:

CAST (Binary collapsing approach): for each individual, count if they have a rare allele at any of the sites
MZ Test/GRANVIL (Count based collapsing): for each individual, count the total number of sites where a rare allele is present
Weighted burden test: for each individual, take a weighted count of the rare alleles across sites (for the weights, use weights <- dbeta(MAF, 1, 25))

For each approach, first need to generate the burden scores.

# CAST : count number of rare alleles for each person and determine if it is > 0
sum_alleles_per_sample <- rowSums(G)
burden.cast <- as.numeric( sum_alleles_per_sample > 0 )

# MZ : count number of sites with rare alleles for each person
burden.mz <- rowSums(G > 0)

# Weighted burden : weighted sum of genotype counts across sites
weights <- dbeta(maf, 1, 25) 
burden.weighted <- G %*% weights

Run a test for association between the phenotype and each burden score using the lm() R function, e.g.

# e.g. for CAST
summary(lm(pheno_file$Pheno ~ burden.cast))

Are there notable differences across the test results?

Now use SKAT to test for an association.

# first fit the null model
skat.null <- SKAT_Null_Model( pheno_file$Pheno ~ 1 , out_type = "C")
# Run SKAT association test (returns a list - p-value is in `$p.value`)
skat_sumstats <- SKAT(G, skat.null )
str(skat_sumstats)
skat_sumstats$p.value

Run the omnibus SKAT, but consider setting \(\rho\) (i.e.r.corr) to 0 and then 1.

Compare the results to using the CAST,MZ/GRANVIL and Weighted burden collapsing approaches in Question 5 as well as SKAT in Question 6. What tests do these \(\rho\) values correspond to?

# Example for rho = 0
rho <- 0
skat_sumstats_rho <- SKAT(G, skat.null, r.corr = rho )
skat_sumstats_rho$p.value

Now consider the omnibus version of SKAT, but use the “optimal.adj” approach which searches across a range of rho values.

# Run SKATO association test using grid of rho values
skat_sumstats_rho_grid <- SKAT(G, skat.null, method="optimal.adj")
skat_sumstats_rho_grid$p.value

Run ACATV on the single variant p-values. The basic command would look like

# `weights` vector is from Question 5
acat.weights <- weights * weights * maf * (1 - maf)
ACAT(sv_pvals$P, weights = acat.weights)

Run ACATO combining the SKAT and BURDEN p-values (from Question 7) with the ACATV p-value (from Question 9).

# Fill in the p-values
SKAT_pvalue <- 
Burden_pvalue <- 
ACATV_pvalue <-
# compute ACATO
ACAT( c(SKAT_pvalue, Burden_pvalue, ACATV_pvalue) )

Extra

You can explore the impact of different genetic architectures by analyzing 3 simulated phenotypes in “data/rv_pheno_extended.txt”. More specifically, you can run the same analyses done above for each of these traits and compare the performance of the tests relative to how variant effects were simulated:

“Pheno_sparse”: a single causal variant (i.e. sparse genetic architecture)
“Pheno_burden_skat”: 60% of the variant are causal, and they all have positive effects on the trait
“Pheno_skat”: 20% of the variants are causal and the effect direction (+/-) is randomly assigned

sessionInfo()

R version 4.3.0 (2023-04-21)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS 14.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_3.4.2     ACAT_0.91         SKAT_2.2.5        RSpectra_0.16-1  
[5] SPAtest_3.1.2     Matrix_1.5-4      BEDMatrix_2.0.3   dplyr_1.1.2      
[9] data.table_1.14.8

loaded via a namespace (and not attached):
 [1] sass_0.4.6       utf8_1.2.3       generics_0.1.3   stringi_1.7.12  
 [5] lattice_0.21-8   digest_0.6.31    magrittr_2.0.3   evaluate_0.21   
 [9] grid_4.3.0       fastmap_1.1.1    rprojroot_2.0.3  workflowr_1.7.0 
[13] jsonlite_1.8.5   whisker_0.4.1    promises_1.2.0.1 fansi_1.0.4     
[17] scales_1.2.1     jquerylib_0.1.4  cli_3.6.1        rlang_1.1.1     
[21] munsell_0.5.0    withr_2.5.0      cachem_1.0.8     yaml_2.3.7      
[25] tools_4.3.0      colorspace_2.1-0 httpuv_1.6.11    crochet_2.3.0   
[29] vctrs_0.6.2      R6_2.5.1         lifecycle_1.0.3  git2r_0.32.0    
[33] stringr_1.5.0    fs_1.6.2         pkgconfig_2.0.3  pillar_1.9.0    
[37] bslib_0.5.0      later_1.3.1      gtable_0.3.3     glue_1.6.2      
[41] Rcpp_1.0.10      xfun_0.39        tibble_3.2.1     tidyselect_1.2.0
[45] highr_0.10       rstudioapi_0.14  knitr_1.43       htmltools_0.5.5 
[49] rmarkdown_2.22   compiler_4.3.0