-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathCOITR.r
55 lines (45 loc) · 1.93 KB
/
COITR.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
43
44
45
46
47
48
49
50
51
52
53
54
55
library(optparse)
library(RJSONIO)
library(stringr)
library(e1071)
library(SDMTools)
option_list <- list(
make_option("--ind1", action = "store", default = NULL, type = 'character', help = "json-file of the first parent"),
make_option("--ind2", action = "store", default = NULL, type = 'character', help = "json-file of the second parent"),
make_option("--phen", action = "store", default = "Phenotypes.json", type = 'character', help = "probabilities"),
make_option("--out", action = "store", default = "output_all/", type = 'character', help = "folder for results")
)
args <- parse_args(OptionParser(option_list = option_list))
source_local <- function(fname){
argv <- commandArgs(trailingOnly = FALSE)
base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
source(paste(base_dir, fname, sep = "/"))
}
source_local('functions.r')
#get probabilities
phen <- fromJSON(args$phen)
#get parents' genotypes
ind1 <- fromJSON(args$ind1)
ind2 <- fromJSON(args$ind2)
#result for child
result_ch <- vector(mode = "list", length = 3)
names(result_ch) <- c("eyes","hair","skin")
result_ch <- lapply(names(result_ch), function(x){
#list of snps
snps <- names(phen[[x]])[substring(names(phen[[x]]), 1, 2) %in% "rs" | substring(names(phen[[x]]), 1, 3) %in% "gen"]
mom <- rep(NA, length(snps))
names(mom) <- snps
mom[names(ind1[[x]]$Snps)] <- ind1[[x]]$Snps
mom = sapply(mom,drop.null)
mom[which(mom == "NA")] <- NA
dad <- rep(NA, length(snps))
names(dad) <- snps
dad[names(ind2[[x]]$Snps)] <- ind2[[x]]$Snps
dad = sapply(dad,drop.null)
dad[which(dad == "NA")] <- NA
res_mid <- child_pred(mother = mom, father = dad, snps = snps, trait = x, phenes = phen)
result_ch[[x]] <- stat_child(comb_prob = res_mid$comb_prob, comb_prob_2 = res_mid$comb_prob_2, phenes = phen, trait = x)
})
names(result_ch) <- c("eyes","hair","skin")
exportJson <- toJSON(result_ch)
write(exportJson, file = paste0(args$out,"/","child.json"))