Using SNPknock with Genetic Data

Matteo Sesia (msesia@stanford.edu)

2019-11-06

This vignette illustrates the basic usage of the SNPknock package in combination with the phasing software fastphase (Scheet and Stephens 2006) to create knockoffs of unphased genotypes or phased haplotypes (Sesia, Sabatti, and Candès 2019, Sesia et al. (2019)). Since fastphase is not available as an R package, this functionality of SNPknock requires the user to first obtain a copy of fastphase.

To learn more about the application of SNPknock to large genome-wide association studies (Sesia et al. 2019), visit: https://msesia.github.io/knockoffzoom/.

Obtaining fastphase

fastphase is a phasing and imputation tool based on the hidden Markov model described in (Scheet and Stephens 2006).

Binary executables for Linux and Mac OS are available from http://scheet.org/software.html.

Before continuing with this tutorial, download the fastphase tarball from the above link and extract the fastphase executable file into a convenient directory (e.g. “~/bin/”).

Knockoffs for unphased genotypes

Fitting the hidden Markov model on genotype data

A small synthetic dataset of 1454 unphased genotype SNPs from 100 individuals can be found in the package installation directory. We can load it with:

library(SNPknock)
X_file = system.file("extdata", "genotypes.RData", package = "SNPknock")
load(X_file)
table(X)

## X
##      0      1      2 
## 130548  90112  25540

Below, we show how to fit a hidden Markov model to this data, with the help of fastphase. Since fastphase takes as input genotype sequences in “.inp” format, we must first convert the X matrix by calling writeXtoInp. By default, this function will write onto a temporary file in the R temporary directory.

# Convert X into the suitable fastphase input format, write it into a temporary file
# and return the path to that file.
Xinp_file = writeXtoInp(X)

Assuming that we have already downloaded fastphase, we can call it to fit the hidden Markov model to X.

fp_path  = "~/bin/fastphase" # Path to the fastphase executable
# Call fastphase and return the path to the parameter estimate files
fp_outPath = runFastPhase(fp_path, Xinp_file, K = 12, numit = 15)

## ~/bin/fastphase -Pp -T1 -K12 -g -H-4 -C15 -S1 -o'/tmp/RtmpRliKJD/file2ee3295d8f00' /tmp/RtmpRliKJD/file2ee31078d299.inp

Above, the SNPknock package could not find fastphase because we did not provide the correct path (we cannot include third-party executable files within this package). However, if you install fastphase separately and provide SNPknock with the correct path, this will work.

If the previous step worked for you, you can find the parameter estimates produced by fastphase in the following files:

r_file = paste(fp_outPath, "_rhat.txt", sep="")
alpha_file = paste(fp_outPath, "_alphahat.txt", sep="")
theta_file = paste(fp_outPath, "_thetahat.txt", sep="")
char_file  = paste(fp_outPath, "_origchars", sep="")

Otherwise, for the sake of this tutorial, you can use the example parameter files provided in the package installation directory:

r_file = system.file("extdata", "genotypes_rhat.txt", package = "SNPknock")
alpha_file = system.file("extdata", "genotypes_alphahat.txt", package = "SNPknock")
theta_file = system.file("extdata", "genotypes_thetahat.txt", package = "SNPknock")
char_file  = system.file("extdata", "genotypes_origchars", package = "SNPknock")

Then, we can construct the hidden Markov model with:

hmm = loadHMM(r_file, alpha_file, theta_file, char_file)

Generating knockoff genotypes

Finally, we can use the hidden Markov model created above to generate knockoffs.

Xk = knockoffGenotypes(X, hmm$r, hmm$alpha, hmm$theta)
table(Xk)

## Xk
##      0      1      2 
## 130114  90466  25620

Knockoffs for phased haplotypes

Fitting the hidden Markov model on haplotype data

A small synthetic dataset of 1454 phased haplotype SNPs from 100 individuals can be found in the package installation directory. We can load it with:

library(SNPknock)
H_file = system.file("extdata", "haplotypes.RData", package = "SNPknock")
load(H_file)
table(H)

## H
##      0      1 
## 351208 141192

Below, we show how to fit a hidden Markov model to this data, with the help of fastphase. Since fastphase takes as input haplotype sequences in “.inp” format, we must first convert the H matrix by calling writeXtoInp. By default, this function will write onto a temporary file in the R temporary directory.

# Convert X into the suitable fastphase input format, write it into a temporary file
# and return the path to that file.
Hinp_file = writeXtoInp(H, phased = TRUE)

Assuming that we have already downloaded fastphase, we can call it to fit the hidden Markov model to X.

fp_path  = "~/bin/fastphase" # Path to the fastphase executable
# Call fastphase and return the path to the parameter estimate files
fp_outPath = runFastPhase(fp_path, Hinp_file, K = 12, numit = 15, phased = TRUE)

## ~/bin/fastphase -Pp -T1 -K12 -g -H-4 -C15 -B -S1 -o'/tmp/RtmpRliKJD/file2ee335b3c2ea' /tmp/RtmpRliKJD/file2ee347eddf25.inp

Above, the SNPknock package could not find fastphase because we did not provide the correct path (we cannot include third-party executable files within this package). However, if you install fastphase separately and provide SNPknock with the correct path, this will work.

If the previous step worked for you, you can find the parameter estimates produced by fastphase in the following files:

r_file = paste(fp_outPath, "_rhat.txt", sep="")
alpha_file = paste(fp_outPath, "_alphahat.txt", sep="")
theta_file = paste(fp_outPath, "_thetahat.txt", sep="")
char_file  = paste(fp_outPath, "_origchars", sep="")

Otherwise, for the sake of this tutorial, you can use the example parameter files provided in the package installation directory:

r_file = system.file("extdata", "haplotypes_rhat.txt", package = "SNPknock")
alpha_file = system.file("extdata", "haplotypes_alphahat.txt", package = "SNPknock")
theta_file = system.file("extdata", "haplotypes_thetahat.txt", package = "SNPknock")
char_file  = system.file("extdata", "haplotypes_origchars", package = "SNPknock")

Then, we can construct the hidden Markov model with:

hmm = loadHMM(r_file, alpha_file, theta_file, char_file)

Generating knockoff haplotypes

Finally, we can use the hidden Markov model created above to generate knockoffs.

Hk = knockoffHaplotypes(H, hmm$r, hmm$alpha, hmm$theta)
table(Hk)

## Hk
##      0      1 
## 351155 141245

See also

If you want to see some basic usage of SNPknock, see the introductory vignette.

If you want to learn about SNPknock for large genome-wide association studies (Sesia et al. 2019), see https://msesia.github.io/knockoffzoom/.

Scheet, P., and M. Stephens. 2006. “A Fast and Flexible Statistical Model for Large-Scale Population Genotype Data: Applications to Inferring Missing Genotypes and Haplotypic Phase.” Am. J. Hum. Genet. 78:629–44. https://doi.org/10.1086/502802.

Sesia, M., E. Katsevich, S. Bates, E. Candès, and C. Sabatti. 2019. “Multi-Resolution Localization of Causal Variants Across the Genome.” bioRxiv. Cold Spring Harbor Laboratory. https://doi.org/10.1101/631390.

Sesia, M., C. Sabatti, and E. J. Candès. 2019. “Gene Hunting with Hidden Markov Model Knockoffs.” Biometrika 106:1–18. https://doi.org/10.1093/biomet/asy033.