cond.test <- function(meta.stat, stat, reference, nsamples, model, ref.info, options){
res <- NULL
for(i in 1:nrow(model)){
msg <- paste0("Performing conditional test (", i, '/', nrow(model), '): ', date())
message(msg)
cond <- unlist(strsplit(model$cond[i], ','), use.names = FALSE)
test <- unlist(strsplit(model$test[i], ','), use.names = FALSE)
snps <- c(cond, test)
ref <- ref.info[snps, ]
# Reference.ID refer to plink files that contain genotypes of snps
ref.id <- unique(ref$Reference.ID)
ref.geno <- NULL
for(j in ref.id){
if(is.null(ref.geno)){
ref.geno <- read.bed(reference$bed[j], reference$bim[j], reference$fam[j], sel.snps = ref$SNP[ref$Reference.ID == j])
}else{
ref.geno <- cbind(ref.geno, read.bed(reference$bed[j], reference$bim[j], reference$fam[j], sel.snps = ref$SNP[ref$Reference.ID == j]))
}
}
names(snps) <- ref$SNP
colnames(ref.geno) <- snps[colnames(ref.geno)]
maf <- apply(ref.geno, 2, function(u){m <- mean(u, na.rm = TRUE)/2; min(m, 1-m)})
rm.snps <- colnames(ref.geno)[which(maf < options$maf)]
if(length(rm.snps) > 0){
cond <- setdiff(cond, rm.snps)
test <- setdiff(test, rm.snps)
if(length(cond) == 0 || length(test) == 0){
next
}
model$cond[i] <- paste(cond, sep = '', collapse = ',')
model$test[i] <- paste(test, sep = '', collapse = ',')
snps <- c(cond, test)
ref.geno <- ref.geno[, snps, drop = FALSE]
}
for(s in test){
res <- rbind(res, recover.stat(meta.stat, stat, nsamples, cond, s, ref.geno, ref.info, options))
}
#print(dim(res))
}
res
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.