3 Computing a GRM
We can use the SNPRelate package to compute a Genetic Relationship matrix (GRM).
library(SeqArray)
data.path <- "https://github.com/UW-GAC/analysis_pipeline/raw/master/testdata"
gdsfile <- "data/1KG_phase3_subset_chr1.gds"
if (!dir.exists("data")) dir.create("data")
if (!file.exists(gdsfile)) download.file(file.path(data.path, basename(gdsfile)), gdsfile)
gds <- seqOpen(gdsfile)
library(SNPRelate)
grm <- snpgdsGRM(gds, method="GCTA")## Genetic Relationship Matrix (GRM, GCTA):
## Calculating allele counts/frequencies ...
##
[..................................................] 0%, ETC: ---
[==================================================] 100%, completed in 0s
## Excluding 13 SNVs (monomorphic: TRUE, MAF: NaN, missing rate: NaN)
## Working space: 1,126 samples, 1,107 SNVs
## using 1 (CPU) core
## CPU capabilities: Double-Precision SSE2
## Sun Jul 22 22:14:56 2018 (internal increment: 680)
##
[..................................................] 0%, ETC: ---
[==================================================] 100%, completed in 1s
## Sun Jul 22 22:14:57 2018 Done.names(grm)## [1] "sample.id" "snp.id" "method" "grm"dim(grm$grm)## [1] 1126 1126seqClose(gds)