forked from Blue-Lab/assoc_pipeline
-
Notifications
You must be signed in to change notification settings - Fork 0
/
king_grm.R
executable file
·42 lines (33 loc) · 1.21 KB
/
king_grm.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#! /usr/bin/env Rscript
library(argparser)
library(magrittr)
argp <- arg_parser("Run KING-robust") %>%
add_argument("gds_file", help = "GDS file") %>%
add_argument("--out_prefix", help = "Prefix for output files",
default = "") %>%
add_argument("--variant_id", help = "File with vector of variant IDs") %>%
add_argument("--sample_id", help = "File with vector of sample IDs") %>%
add_argument("--num_core", help = "num.thread argument for snpgdsIBDKING (if NA, detect the number of cores automatically)")
argv <- parse_args(argp)
sessionInfo()
print(argv)
library(SNPRelate)
library(SeqArray)
if (!is.na(argv$variant_id)) {
variant_id <- readRDS(argv$variant_id)
} else {
variant_id <- NULL
}
if (!is.na(argv$sample_id)) {
sample_id <- readRDS(argv$sample_id)
} else {
sample_id <- NULL
}
gds <- seqOpen(argv$gds_file)
king <- snpgdsIBDKING(gds, snp.id = variant_id, sample.id = sample_id,
type = "KING-robust",
num.thread = as.numeric(argv$num_core))
rownames(king$kinship) <- king$sample.id
colnames(king$kinship) <- king$sample.id
saveRDS(king, paste0(argv$out_prefix, "king_out.rds"))
saveRDS(king$kinship, paste0(argv$out_prefix, "king_grm.rds"))