Nothing
setGeneric(name = "enrichSNP", def = function (List, Table, EnrichmentRatio, Z, PValue, Resampling) {standardGeneric("enrichSNP")})
setGeneric(name = "chromosome", def = function (Data, LD, eSNP, xSNP) {standardGeneric("chromosome")})
setGeneric(name = "enrichment", def = function (Loss, Call, eSNP, xSNP, Chromosomes) {standardGeneric("enrichment")})
setGeneric(name = "doLDblock", def = function (object, mc.cores = 1) {standardGeneric("doLDblock")})
setGeneric(name = "excludeSNP", def = function (object, excludeFile, mc.cores = 1) {standardGeneric("excludeSNP")})
setGeneric(name = "computeER", def = function (object, sigThresh = 0.05, mc.cores = 1) {standardGeneric("computeER")})
setGeneric(name = "reSample", def = function (object, nSample = 100, empiricPvalue = TRUE, MAFpool = c(0.05, 0.10, 0.2, 0.3, 0.4, 0.5), mc.cores = 1, onlyGenome = TRUE, ...) {standardGeneric("reSample")})
setGeneric(name = "compareEnrichment", def = function (object.x, object.y, pattern = "Chrom", nSample = 100, empiricPvalue = TRUE, mc.cores = 1, onlyGenome = TRUE) {standardGeneric("compareEnrichment")})
setGeneric(name = "is.enrichment", def = function (object) {standardGeneric("is.enrichment")})
setGeneric(name = "is.chromosome", def = function (object) {standardGeneric("is.chromosome")})
setGeneric(name = "is.EnrichSNP", def = function (object) {standardGeneric("is.EnrichSNP")})
setGeneric(name = "reset", def = function (object, i) {standardGeneric("reset")})
setGeneric(name = "getEnrichSNP", def = function (object, type = "eSNP") {standardGeneric("getEnrichSNP")})
setMethod(f = "reSample", signature = "ANY", definition = function (object, nSample, sigThresh, MAFpool, mc.cores) {
if (!(is.enrichment(object) & is.chromosome(object))) {
stop('[Method:reSample] not available for "', class(object), '" object.', call. = FALSE)
} else {}
})
setMethod(f = "excludeSNP", signature = "ANY", definition = function (object, excludeFile, mc.cores = 1) {
if (!is.enrichment(object)) {
stop('[Method:excludeSNP] not available for "', class(object), '" object.', call. = FALSE)
} else {}
})
setMethod(f = "reset", signature = "ANY", definition = function (object, i) {
if (!(is.enrichment(object) & is.chromosome(object))) {
stop('[Method:reset] not available for "', class(object), '" object.', call. = FALSE)
} else {}
})
setMethod(f = "getEnrichSNP", signature = "ANY", definition = function (object, type = "eSNP") {
if (!(is.enrichment(object))) {
stop('[Method:getEnrichSNP] not available for "', class(object), '" object.', call. = FALSE)
} else {}
})
.verbose <- function (expr) {return(invisible(capture.output(expr)))}
GC <- function (verbose = getOption("verbose"), reset = FALSE) {
while (!identical(gc(verbose, reset)[, 4], gc(verbose, reset)[, 4])) {}
return(gc(verbose, reset))
}
.checkFilePath <- function (path) {
# END <- unlist(regmatches(path, regexec(".$", path)))
# START <- unlist(regmatches(path, regexec("^.", path)))
# if (START != "/") {
# path <- paste0("/", path)
# } else {}
# if (END != "/") {
# path <- paste0(path, "/")
# } else {}
path <- gsub("/*$", "/", path)
return(path)
}
maxCores <- function (mc.cores = 1) {
if (Sys.info()[["sysname"]] == "Linux") {
nbCores <- detectCores()
mc.cores.old <- mc.cores
if (file.exists("/proc/meminfo")) {
memInfo <- readLines("/proc/meminfo")
sysMemFree <- memInfo[grep('^MemFree:', memInfo)]
sysMemCached <- memInfo[grep('^Cached:', memInfo)]
sysMemAvailable <- 0.95*(as.numeric(gsub("[^0-9]*([0-9]*)", "\\1", sysMemFree)) + as.numeric(gsub("[^0-9]*([0-9]*)", "\\1", sysMemCached)))
sysProc <- as.numeric(unlist(strsplit(system(paste("ps v", Sys.getpid()), intern = TRUE)[2], " +"))[8])
mc.cores <- max(min(as.numeric(mc.cores), floor(sysMemAvailable/sysProc)), 1)
if (mc.cores > nbCores) {
mc.cores <- nbCores
} else {}
if (mc.cores != mc.cores.old) {
warning(paste0('To avoid memory overload "mc.cores" was decreased to "', mc.cores, '".'), call. = FALSE)
} else {}
} else {
mc.cores <- ifelse(mc.cores.old>nbCores, nbCores, mc.cores.old)
}
} else {
mc.cores <- 1
}
return(mc.cores)
}
mclapply2 <- function (X, FUN, ..., mc.preschedule = TRUE, mc.set.seed = TRUE, mc.silent = FALSE, mc.cores = getOption("mc.cores", 2L), mc.cleanup = TRUE, mc.allow.recursive = TRUE) {
if (Sys.info()[["sysname"]] != "Linux") {
mc.cores <- 1
} else {
mc.cores <- min(detectCores(), mc.cores)
}
return(mclapply(X = X, FUN = FUN, ...,
mc.preschedule = mc.preschedule, mc.set.seed = mc.set.seed, mc.silent = mc.silent,
mc.cores = maxCores(mc.cores), mc.cleanup = mc.cleanup, mc.allow.recursive = mc.allow.recursive))
}
.writeSignal <- function (pattern, snpInfoDir, signalFile) {
tmpDir <- gsub("\\\\", "/", tempdir())
if (length(unlist(strsplit(readLines(signalFile, n = 1), split = "\t")))>1) {
signal <- read.delim(file = signalFile, header = TRUE, sep = "\t", stringsAsFactors = FALSE,
colClasses = c("character", "numeric"), na.string = c("NA", ""),
check.names = FALSE, strip.white = TRUE, col.names = c("SNP", "PVALUE"))
} else {
if (length(unlist(strsplit(readLines(signalFile, n = 1), split = " ")))>1) {
signal <- read.delim(file = signalFile, header = TRUE, sep = " ", stringsAsFactors = FALSE,
colClasses = c("character", "numeric"), na.string = c("NA", ""),
check.names = FALSE, strip.white = TRUE, col.names = c("SNP", "PVALUE"))
} else {
stop('[Enrichment:initFiles] only " " and "\t" are allowed as columns separator in Signal file.', call. = FALSE)
}
}
chrom.bim <- read.delim(file = paste0(snpInfoDir, pattern, ".bim"), header = FALSE, stringsAsFactors = FALSE,
colClasses = c("numeric", "character", "NULL", "NULL", "NULL", "NULL"), na.string = c("NA", ""),
check.names = FALSE, strip.white = TRUE, col.names = c("CHR", "SNP", "", "", "", ""))
signal <- merge(signal, chrom.bim, by = "SNP")[, c(3, 1, 2)]
eval(parse(text = paste0('write.table(signal, file = "', tmpDir, '/snpEnrichment/', pattern,
'.signal", quote = FALSE, row.names = FALSE, col.names = FALSE, sep = "\t")')))
return(invisible())
}
.readSignal <- function (pattern) {
tmpDir <- gsub("\\\\", "/", tempdir())
signal <- read.delim(file = paste0(tmpDir, "/snpEnrichment/", pattern, ".signal"), header = TRUE, stringsAsFactors = FALSE,
colClasses = c("NULL", "character", "numeric"), na.string = c("NA", ""),
check.names = FALSE, strip.white = TRUE, col.names = c("", "SNP", "PVALUE"))
return(signal)
}
.writeFreq <- function (pattern, snpInfoDir) {
tmpDir <- gsub("\\\\", "/", tempdir())
IN <- paste0(snpInfoDir, pattern)
OUT <- paste0(tmpDir, "/snpEnrichment/", pattern)
plinkData <- read.plink(bed = paste0(IN, ".bed"), bim = paste0(IN, ".bim"), fam = paste0(IN, ".fam"), select.snps = .readSignal(pattern)[, "SNP"])
plinkFreq <- col.summary(plinkData$genotypes)
plinkFreq <- cbind(snp.name = rownames(plinkFreq), MAF = plinkFreq[, "MAF"])
plinkRes <- merge(plinkData$map, plinkFreq, by = "snp.name")
plinkRes <- plinkRes[, c("chromosome", "snp.name", "position", "MAF")]
plinkRes[, "MAF"] <- as.numeric(as.character(plinkRes[, "MAF"]))
colnames(plinkRes) <- c("CHR", "SNP", "POS", "MAF")
write.table(plinkRes, paste0(OUT, ".all"), row.names = FALSE, sep = "\t")
return(invisible())
}
.checkSnpInfoDir <- function (snpInfoDir) {
snpInfoDir <- .checkFilePath(snpInfoDir)
if (length(list.files(snpInfoDir, pattern = "*.bim"))!=22) {
stop(paste0("[Enrichment:initFiles] Only ", length(list.files(snpInfoDir, pattern = "*.bim")), " 'bim' files found when 22 is needed."), call. = FALSE)
} else {}
if (length(list.files(snpInfoDir, pattern = "*.bed"))!=22) {
stop(paste0("[Enrichment:initFiles] Only ", length(list.files(snpInfoDir, pattern = "*.bed")), " 'bed' files found when 22 is needed."), call. = FALSE)
} else {}
if (length(list.files(snpInfoDir, pattern = "*.fam"))!=22) {
stop(paste0("[Enrichment:initFiles] Only ", length(list.files(snpInfoDir, pattern = "*.fam")), " 'fam' files found when 22 is needed."), call. = FALSE)
} else {}
return(invisible())
}
.checkSignalFile <- function (signalFile) {
if (!file.exists(signalFile)) {
stop(paste0("[Enrichment:initFiles] ", signalFile, " doesn't exist."), call. = FALSE)
} else {}
return(invisible())
}
.checkSnpListDir <- function (snpListDir, pattern) {
snpListDir <- .checkFilePath(snpListDir)
if (length(list.files(snpListDir, pattern = pattern))==0) {
stop(paste0("[Enrichment:readEnrichment] No snp list file found when at least one is needed."), call. = FALSE)
} else {}
return(invisible())
}
initFiles <- function (pattern = "Chrom", snpInfoDir, signalFile, mc.cores = 1) {
if (missing(snpInfoDir) | missing(signalFile)) {
stop('[Enrichment:initFiles] argument(s) missing.', call. = FALSE)
} else {}
snpInfoDir <- .checkFilePath(snpInfoDir)
FILES <- list.files(snpInfoDir, pattern = ".bim")
.checkSnpInfoDir(snpInfoDir)
.checkSignalFile(signalFile)
tmpDir <- gsub("\\\\", "/", tempdir())
dir.create(paste0(tmpDir, "/snpEnrichment/"), showWarnings = FALSE)
cat("All files are ready for chromosome:\n ")
resParallel <- mclapply2(X = seq(22), mc.cores = min(22, mc.cores), FUN = function (iChr) {
# newPattern <- unlist(strsplit(grep(paste0(pattern, iChr, "[^0-9]*.bim"), FILES, value = TRUE), ".bim"))[1]
newPattern <- gsub(".bim", "", grep(paste0(pattern, iChr, "[^0-9]"), FILES, value = TRUE))
err1 <- try(.writeSignal(pattern = newPattern, snpInfoDir = snpInfoDir, signalFile = signalFile), silent = TRUE)
err2 <- try(.writeFreq(pattern = newPattern, snpInfoDir = snpInfoDir), silent = TRUE)
cat(iChr, " ", sep = "")
if (class(err1)=="try-error" | class(err2)=="try-error") {
return(invisible("ERROR"))
} else {
return(invisible())
}
})
if (any(unlist(resParallel)=="ERROR")) {
stop("[Enrichment:initFiles] initialize files failed.", call. = FALSE)
} else {}
cat("\n\n")
return(invisible())
}
writeLD <- function (pattern = "Chrom", snpInfoDir, signalFile, ldDir = NULL, ldThresh = 0.8, depth = 1000, mc.cores = 1) {
if (missing(pattern) | missing(snpInfoDir) | missing(signalFile) | missing(ldThresh)) {
stop('[Enrichment:writeLD] argument(s) missing.', call. = FALSE)
} else {}
tmpDir <- gsub("\\\\", "/", tempdir())
dir.create(paste0(tmpDir, "/snpEnrichment/"), showWarnings = FALSE)
snpInfoDir <- .checkFilePath(snpInfoDir)
.checkSnpInfoDir(snpInfoDir)
.checkSignalFile(signalFile)
if (missing(ldDir) | is.null(ldDir)) {
ldDir <- paste0(tmpDir, "/snpEnrichment/")
} else {
ldDir <- .checkFilePath(ldDir)
}
FILES <- list.files(snpInfoDir, pattern = ".bim")
cat("Compute LD for chromosome:\n ")
resParallel <- mclapply2(X = seq(22), mc.cores = min(22, mc.cores), FUN = function (iChr) {
# newPattern <- unlist(strsplit(grep(paste0(pattern, iChr, "[^0-9]*.bim"), FILES, value = TRUE), ".bim"))[1]
newPattern <- gsub(".bim", "", grep(paste0(pattern, iChr, "[^0-9]"), FILES, value = TRUE))
# isThereSignals <- grep(".signal", list.files(paste0(tmpDir, "/snpEnrichment/"), full.names = TRUE), value = TRUE)
isThereSignals <- list.files(paste0(tmpDir, "/snpEnrichment/"), full.names = TRUE, pattern = ".signal")
if (length(isThereSignals) != 22) {
.writeSignal(pattern = newPattern, snpInfoDir, signalFile)
} else {}
IN <- paste0(snpInfoDir, newPattern)
OUT <- paste0(tmpDir, "/snpEnrichment/", newPattern)
plinkData <- read.plink(bed = paste0(IN, ".bed"), bim = paste0(IN, ".bim"), fam = paste0(IN, ".fam"), select.snps = .readSignal(newPattern)[, "SNP"])
ldData <- ld(x = plinkData$genotypes, depth = min(ncol(plinkData$genotypes)-1, depth), stats = "R.squared")
if (any(isNA <- is.na(ldData))) {
ldData <- replace(ldData, grep(TRUE, isNA), 0)
} else {}
ldData <- apply(ldData, 1, function (li) { which(li>ldThresh) })
resLD <- data.frame(matrix(unlist(strsplit(names(unlist(ldData)), "\\.")), ncol = 2, byrow = TRUE, dimnames = list(NULL, c("SNP_A", "SNP_B"))), stringsAsFactors = FALSE)
resLD <- resLD[which(resLD[, 1]!=resLD[, 2]), ]
### Check LD distance
# snpNames <- grep("[sS][nN][pP]", colnames(plinkData$map), value = TRUE)
# posNames <- grep("[pP][oO][sS].*", colnames(plinkData$map), value = TRUE)
# tmpPlink <- plinkData$map[which(plinkData$map[, snpNames]%in%unique(unlist(resLD[, c("SNP_A", "SNP_B")]))), c(snpNames, posNames)]
# resLD[, "Dist"] <- apply(resLD[, c("SNP_A", "SNP_B")], 1, function (li) {
# abs(tmpPlink[tmpPlink[, snpNames]%in%li["SNP_A"], posNames] - tmpPlink[tmpPlink[, snpNames]%in%li["SNP_B"], posNames])
# })
# table(resLD[, "Dist"]<1000000)
write.table(resLD[, c("SNP_A", "SNP_B")], file = paste0(ldDir, newPattern, ".ld"), sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
cat(iChr, " ", sep = "")
})
cat("\n\n")
return(invisible())
}
.readSNP <- function (pattern, snpListDir) {
# snpListFile <- grep(pattern, sort(list.files(snpListDir, full.names = TRUE)), value = TRUE)[1]
snpListFile <- list.files(gsub("/$", "", snpListDir), pattern = paste0(pattern, "[^0-9]"), full.names = TRUE)
if (is.na(snpListFile)) {
snpList <- data.frame()
} else {
snpList <- read.delim(file = snpListFile, header = FALSE, stringsAsFactors = FALSE,
na.string = c("NA", ""), check.names = FALSE, strip.white = TRUE)
switch(EXPR = as.character(ncol(snpList)),
"1" = {colnames(snpList) <- "SNP"},
"2" = {colnames(snpList) <- c("CHR", "SNP")},
"3" = {colnames(snpList) <- c("CHR", "SNP", "TRANSCRIPT")},
colnames(snpList)[1:3] <- c("CHR", "SNP", "TRANSCRIPT")
)
if (ncol(snpList)>=3) {
snpList <- snpList[, c("SNP", "TRANSCRIPT")]
} else {}
if (nrow(snpList)!=0) {
snpList[, "eSNP"] <- 1
} else {}
}
return(snpList)
}
.readTranscript <- function (transcriptFile) {
if (all(class(try(close(file(transcriptFile)), silent = TRUE))!="try-error")) {
transcript <- read.delim(file = transcriptFile, header = TRUE, stringsAsFactors = FALSE, na.string = c("NA", ""), check.names = FALSE, strip.white = TRUE)
transcript <- transcript[, 1:4]
colnames(transcript) <- c("TRANSCRIPT", "CHR", "START", "END")
} else {
if (class(transcriptFile) %in% c("matrix", "data.frame")) {
transcript <- transcriptFile[, 1:4]
colnames(transcript) <- c("TRANSCRIPT", "CHR", "START", "END")
} else {
stop('[Enrichment:readEnrichment] "transcriptFile" required "matrix", "data.frame" or "txt file".', call. = FALSE)
}
}
return(transcript)
}
.readFreq <- function (pattern, snpInfoDir) {
tmpDir <- gsub("\\\\", "/", tempdir())
freq <- read.delim(file = paste0(tmpDir, "/snpEnrichment/", pattern, ".all"), header = TRUE, stringsAsFactors = FALSE,
colClasses = c("numeric", "character", "integer", "numeric"), na.string = c("NA", ""),
check.names = FALSE, strip.white = TRUE, col.names = c("CHR", "SNP", "POS", "MAF"))
return(freq)
}
.readLD <- function (pattern, snpInfoDir, ldDir) {
tmpDir <- gsub("\\\\", "/", tempdir())
if (missing(ldDir) | is.null(ldDir)) {
IN <- paste0(tmpDir, "/snpEnrichment/", pattern, ".ld")
} else {
IN <- paste0(.checkFilePath(ldDir), pattern, ".ld")
}
nbCol <- as.character(ncol(read.table(text = readLines(con = IN, n = 1), stringsAsFactors = FALSE)))
switch(EXPR = nbCol,
"2" = {ldData <- read.delim(file = IN, header = TRUE, stringsAsFactors = FALSE, colClasses = c("character", "character"),
na.string = c("NA", ""), check.names = FALSE, strip.white = TRUE, sep = "\t")
},
"7" = {ldData <- read.delim(file = IN, header = TRUE, stringsAsFactors = FALSE, colClasses = c("NULL", "NULL", "character", "NULL", "NULL", "character", "NULL"),
na.string = c("NA", ""), check.names = FALSE, strip.white = TRUE, sep = "")
},
stop(paste0('[Enrichment:readEnrichment] "', pattern, '.ld" structure must be a matrix file with columns:\n c("SNP_A", "SNP_B") or c("CHR_A", "BP_A", "SNP_A", "CHR_B", "BP_B", "SNP_B", "R2")".'), call. = FALSE)
)
return(ldData)
}
.checkTranscript <- function (data, transcriptFile, distThresh) {
if (any(transcriptFile != FALSE)) {
transcript <- .readTranscript(transcriptFile = transcriptFile)
transcriptCHR <- na.exclude(transcript[transcript[, "CHR"] == unique(data[, "CHR"]), c("START", "END")])
cisFunc <- function (line, distThresh, dataTranscript) {
position <- as.numeric(line[4])
CIS <- any(position > (dataTranscript[, "START"] - distThresh) & (position < (dataTranscript[, "END"] + distThresh)))
return(c(line, CIS = as.numeric(CIS)))
}
temp <- unique(as.data.frame(t(apply(data, 1, cisFunc, distThresh = distThresh*10^3, dataTranscript = transcriptCHR)), stringsAsFactors = FALSE))
temp <- temp[temp[, "CIS"] == 1, -grep("CIS", colnames(temp))]
} else {
temp <- unique(data)
}
return(temp)
}
.readFiles <- function (pattern, snpInfoDir, snpListDir, distThresh) {
fullPattern <- gsub(".bim", "", grep(paste0(pattern, "[^0-9]"), list.files(snpInfoDir, pattern = ".bim"), value = TRUE))
signalPattern <- gsub(".signal", "", grep(paste0(pattern, "[^0-9]"), list.files(paste0(gsub("\\\\", "/", tempdir()), "/snpEnrichment/"), pattern = ".signal"), value = TRUE))
eSNP <- .readSNP(pattern = fullPattern, snpListDir = snpListDir)
signal <- .readSignal(pattern = signalPattern)
plinkData <- .readFreq(pattern = fullPattern, snpInfoDir = snpInfoDir)
signalPlink <- merge(signal, plinkData, by = "SNP")
if (nrow(eSNP)!=0) {
eSNPunique <- eSNP[!duplicated(eSNP[, "SNP"]), ]
snpSignal <- merge(signalPlink, eSNPunique[, c("SNP", "eSNP")], by = "SNP", all.x = TRUE)
snpSignal[, "eSNP"][is.na(snpSignal[, "eSNP"])] <- 0
snpLoss <- c(
length(eSNP[, "SNP"]),
length(unique(eSNP[, "SNP"])),
length(unique(intersect(eSNP[, "SNP"], signal[!is.na(signal[, "PVALUE"]), "SNP"])))
)
} else {
snpSignal <- signalPlink
snpSignal[, "eSNP"] <- 0
snpLoss <- c(0, 0, 0)
}
temp <- unique(snpSignal[, c("SNP", "PVALUE", "CHR", "POS", "MAF", "eSNP")])
data <- transform(temp, SNP = as.character(temp$SNP), PVALUE = as.numeric(temp$PVALUE), CHR = as.numeric(temp$CHR), POS = as.integer(temp$POS), MAF = as.numeric(temp$MAF), eSNP = as.numeric(temp$eSNP))
data[, "xSNP"] <- 0
if (any(duplicated(data[, "SNP"]))) {
cat(data[, "SNP"][duplicated(data[, "SNP"])], "\n")
stop('[Enrichment:readEnrichment] Duplicated SNPs in Signal.', call. = FALSE)
} else {
rownames(data) <- data[, "SNP"]
}
return(list(data = data, snpLoss = snpLoss))
}
.splitByChrom <- function (pattern, snpListFile, directory) {
if (missing(snpListFile)) {
stop('[Enrichment:readEnrichment] argument(s) missing.', call. = FALSE)
} else {}
snpList <- read.delim(file = snpListFile, header = FALSE, stringsAsFactors = FALSE,
na.string = c("NA", ""), check.names = FALSE, strip.white = TRUE)
switch(EXPR = as.character(ncol(snpList)),
"1" = {stop('[Enrichment:readEnrichment] at least two columns are needed: "Chromosome" and "rs name".', call. = FALSE)},
"2" = {colnames(snpList) <- c("CHR", "SNP")},
"3" = {colnames(snpList) <- c("CHR", "SNP", "TRANSCRIPT")},
colnames(snpList)[1:3] <- c("CHR", "SNP", "TRANSCRIPT")
)
if (ncol(snpList)>3) {
snpList <- snpList[, c("CHR", "SNP", "TRANSCRIPT")]
} else {}
filePathDetails <- unlist(strsplit(snpListFile, "/"))
if (is.null(directory)) {
filePath <- paste0(c(filePathDetails[-length(filePathDetails)], ""), collapse = "/")
} else {
dir.create(paste0(directory, "snpList/"), showWarnings = FALSE)
filePath <- paste0(directory, "snpList/")
}
by(snpList, snpList[, "CHR"], function (snpListChr) {
write.table(snpListChr, file = paste0(filePath, pattern, unique(snpListChr[, "CHR"]), "-", filePathDetails[length(filePathDetails)]), sep = "\t", col.names = FALSE, row.names = FALSE, quote = FALSE)
})
# cat(paste0('*** File "', filePathDetails[length(filePathDetails)], '" has been splitted into ', length(unique(snpList[, "CHR"])), ' chromosomes. ***\n\n'))
return(filePath)
}
readEnrichment <- function (pattern = "Chrom", signalFile, transcriptFile = FALSE, snpListDir, snpInfoDir, distThresh = 1000, sigThresh = 0.05, LD = FALSE, ldDir = NULL, mc.cores = 1) {
cat("############# Read Enrichment ##############\n")
if (missing(signalFile) | missing(snpListDir) | missing(snpInfoDir)) {
stop('[Enrichment:readEnrichment] argument(s) missing.', call. = FALSE)
} else {}
tmpDir <- paste0(gsub("\\\\", "/", tempdir()), "/snpEnrichment/")
if (length(list.files(tmpDir, pattern = "\\.signal")) != 22 & length(list.files(tmpDir, pattern = "\\.all")) != 22) {
stop('[Enrichment:readEnrichment] "initFiles" must be run before.', call. = FALSE)
} else {}
if (is.null(ldDir)) {
if (LD & length(list.files(tmpDir, pattern = "\\.ld")) != 22) {
stop('[Enrichment:readEnrichment] "writeLD" must be run before. (Or LD computation by PLINK.)', call. = FALSE)
} else {}
} else {
if (LD & length(list.files(ldDir, pattern = "\\.ld")) != 22) {
stop('[Enrichment:readEnrichment] "writeLD" must be run before. (Or LD computation by PLINK.)', call. = FALSE)
} else {}
}
snpInfoDir <- .checkFilePath(snpInfoDir)
.checkSnpInfoDir(snpInfoDir)
.checkSignalFile(signalFile)
if (length(list.files(snpListDir)) == 0) {
snpListDir <- .splitByChrom(pattern = pattern, snpListFile = snpListDir, directory = tmpDir)
} else {
snpListDir <- .checkFilePath(snpListDir)
.checkSnpListDir(snpListDir, pattern)
}
cat(" Read Chromosomes:\n ")
resParallel <- mclapply2(seq(22), mc.cores = min(22, mc.cores), FUN = function (iChr) {
files <- .readFiles(pattern = paste0(pattern, iChr), snpInfoDir = snpInfoDir, snpListDir = snpListDir, distThresh = distThresh)
if (LD) {
# newPattern <- unlist(strsplit(grep(paste0(pattern, iChr, "[^0-9]*.bim"), list.files(snpInfoDir), value = TRUE), ".bim"))[1]
newPattern <- gsub(".bim", "", grep(paste0(pattern, iChr, "[^0-9]"), list.files(snpInfoDir, pattern = ".bim"), value = TRUE))
linkageData <- .readLD(pattern = newPattern, snpInfoDir = snpInfoDir, ldDir = ldDir)
data <- files$data[order(files$data$POS), ]
snpSignal <- data[, "SNP"]
ldData <- linkageData[, 2]
names(ldData) <- linkageData[, 1]
eSNP <- data[data[, "eSNP"]==1, "SNP"]
ldSNP <- linkageData[linkageData[, "SNP_A"] %in% eSNP, "SNP_B"]
xSNP <- unique(c(eSNP, ldSNP))
data[data[, "SNP"]%in%xSNP, "xSNP"] <- 1
data <- data[!is.na(data[, "PVALUE"]), ]
resChr <- chromosome(Data = data, LD = ldData)
} else {
resChr <- chromosome(Data = files$data)
}
snpLoss <- files$snpLoss
if (all(transcriptFile == FALSE)) {
resCheckData <- resChr@Data
resCheckData <- transform(unique(resCheckData),
SNP = as.character(resCheckData$SNP),
PVALUE = as.numeric(resCheckData$PVALUE),
CHR = as.numeric(resCheckData$CHR),
POS = as.numeric(resCheckData$POS),
MAF = as.numeric(resCheckData$MAF),
eSNP = as.numeric(resCheckData$eSNP),
xSNP = as.numeric(resCheckData$xSNP)
)
signalLoss <- c(NA, length(unique(resChr@Data[, "SNP"])), length(unique(resChr@Data[!is.na(resChr@Data[, "PVALUE"]), "SNP"])))
} else {
resCheckData <- .checkTranscript(data = resChr@Data, transcriptFile = transcriptFile, distThresh = distThresh)
snpLoss <- c(snpLoss, sum(as.numeric(resCheckData[, "eSNP"])))
resCheckData <- transform(unique(resCheckData),
SNP = as.character(resCheckData$SNP),
PVALUE = as.numeric(resCheckData$PVALUE),
CHR = as.numeric(resCheckData$CHR),
POS = as.numeric(resCheckData$POS),
MAF = as.numeric(resCheckData$MAF),
eSNP = as.numeric(resCheckData$eSNP),
xSNP = as.numeric(resCheckData$xSNP)
)
signalLoss <- c(NA, length(unique(resChr@Data[, "SNP"])), length(unique(resChr@Data[!is.na(resChr@Data[, "PVALUE"]), "SNP"])), length(unique(resCheckData[!is.na(resCheckData[, "PVALUE"]), "SNP"])))
}
resChr <- chromosome(Data = resCheckData, LD = resChr@LD)
cat(iChr, " ", sep = "")
return(list(resChr, snpLoss, signalLoss))
})
names(resParallel) <- paste0("Chrom", seq(22))
result = enrichment()
result@Chromosomes <- lapply(resParallel, "[[", 1)
snpLoss <- t(sapply(resParallel, "[[", 2))
signalLoss <- rowSums(sapply(resParallel, "[[", 3))
loss <- rbind(Signal = signalLoss, Genome = apply(snpLoss, 2, sum), snpLoss)
colnames(loss) <- c("Rows", "Unique", "Intersect.Signal", "CIS")[seq(ncol(loss))]
SNPs <- result["List", seq(22)]
result@eSNP@List <- SNPs[["eSNP"]]
result@xSNP@List <- SNPs[["xSNP"]]
if (length(unlist(strsplit(readLines(signalFile, n = 1), split = "\t")))>1) {
signal <- read.delim(file = signalFile, header = TRUE, sep = "\t", colClasses = c("character", "numeric"), na.string = c("NA", ""), check.names = FALSE, strip.white = TRUE, col.names = c("SNP", "PVALUE"), stringsAsFactors = FALSE)
} else {
if (length(unlist(strsplit(readLines(signalFile, n = 1), split = " ")))>1) {
signal <- read.delim(file = signalFile, header = TRUE, sep = " ", colClasses = c("character", "numeric"), na.string = c("NA", ""), check.names = FALSE, strip.white = TRUE, col.names = c("SNP", "PVALUE"), stringsAsFactors = FALSE)
} else {
stop('[Enrichment:readEnrichment] only " " and "\t" are allowed as columns separator in Signal file.', call. = FALSE)
}
}
loss[1, c(1, 2)] <- c(length(signal[, "SNP"]), length(unique(signal[, "SNP"])))
result@Loss <- as.data.frame(loss)
sysCall <- sys.call(sys.parent())
argsSNP <- as.list(sysCall[-1])
formal <- as.list(names(formals(as.character(sysCall))))
names(formal) <- formal
if (is.null(names(argsSNP))) {
names(argsSNP) <- names(formal)[seq_along(argsSNP)]
} else {
emptyNames <- which(names(argsSNP)=="")
names(argsSNP)[emptyNames] <- names(formal)[emptyNames]
}
names(argsSNP)[grep("^$", names(argsSNP))] <- names(formal)[grep("^$", names(argsSNP))]
argsSNP <- c(argsSNP, lapply(formal[!names(formal) %in% names(argsSNP)], as.name))[names(formal)]
paramClass <- sapply(argsSNP, class)
for (iArg in names(formal)) {
if (iArg != "...") {
paramPos <- grep(iArg, names(formal), fixed = TRUE)
argTmp <- argsSNP[[paramPos]]
classTmp <- paramClass[[paramPos]]
switch(EXPR = classTmp,
"character" = {formal[[iArg]] <- argTmp},
"logical" = {formal[[iArg]] <- argTmp},
"numeric" = {formal[[iArg]] <- argTmp},
"integer" = {formal[[iArg]] <- argTmp},
"NULL" = {formal[[iArg]] <- "NULL"},
"call" = {formal[[iArg]] <- eval(argTmp)},
"name" = {
resEval <- eval(argTmp)
switch(EXPR = class(resEval),
"character" = {formal[[iArg]] <- resEval},
"logical" = {formal[[iArg]] <- resEval},
"numeric" = {formal[[iArg]] <- resEval},
"integer" = {formal[[iArg]] <- resEval},
"matrix" = {formal[[iArg]] <- argTmp},
"data.frame" = {formal[[iArg]] <- argTmp},
"Enrichment" = {formal[[iArg]] <- argTmp}
)
}
)
} else {}
}
# unlink(tmpDir, recursive = TRUE) # add parameter to ask if temporary directory should be cleaned
if (is.null(ldDir)) {
formal[["ldDir"]] <- tmpDir
} else {
formal[["ldDir"]] <- .checkFilePath(formal[["ldDir"]])
}
result@Call$readEnrichment <- formal
result <- computeER(object = result, sigThresh = sigThresh, mc.cores = mc.cores)
cat("\n########### Read Enrichment Done ###########\n\n")
return(result)
}
.enrichmentRatio <- function (table) {
tmp <- as.numeric(table)
return((tmp[1]*tmp[4])/(tmp[2]*tmp[3]))
}
.reSample <- function (object, nSample, empiricPvalue, sigThresh, MAFpool, mc.cores) {
nResampling <- nrow(object@eSNP@Resampling)
data <- object["Data"]
chrLD <- object["LD"]
isLD <- length(chrLD) != 0
eEnrichment <- object@eSNP@Table
eEnrichRatio <- object@eSNP@EnrichmentRatio
eSNPlist <- object@eSNP@List
data[, "MAFpool"] <- NA
data[, "MAFpool"] <- as.factor(cut(data[, "MAF"], breaks = MAFpool, labels = FALSE, include.lowest = TRUE))
nPool <- nlevels(data[, "MAFpool"])
eSNPlistPool <- table(data[data[, "eSNP"] == 1, "MAFpool"])
if (nResampling == 0) {
eSNPsample <- NULL
nSampleMin <- min(nSample, max(1000, nSample*10/100))
} else {
eSNPsample <- object@eSNP@Resampling
nSampleMin <- nSample
}
if (isLD) {
if (nResampling == 0) {
xSNPsample <- NULL
} else {
xSNPsample <- object@xSNP@Resampling
}
xEnrichment <- object@xSNP@Table
xEnrichRatio <- object@xSNP@EnrichmentRatio
xSNPlist <- object@xSNP@List
xSNPlistPool <- table(data[data[, "xSNP"] == 1, "MAFpool"], deparse.level = 0)
popSNP4Sample <- split(data[data[, "xSNP"] != 1, "SNP"], data[data[, "xSNP"] != 1, "MAFpool"])
} else {
popSNP4Sample <- split(data[data[, "eSNP"] != 1, "SNP"], data[data[, "eSNP"] != 1, "MAFpool"])
}
assoc <- factor(data[, "PVALUE"] < sigThresh, levels = c(FALSE, TRUE))
SNPlist <- data[, "SNP"]
rm(data)
cat(0, ".. ", sep = "")
resParallel <- mclapply2(X = seq(nSampleMin), mc.cores = mc.cores, FUN = function (i) {
eSNPlistRandom <- unlist(sapply(seq(nPool), function (g) {sample(popSNP4Sample[[g]], min(eSNPlistPool[g], length(popSNP4Sample[[g]])))}))
eSNPenrichStats <- table(assoc, factor(SNPlist%in%eSNPlistRandom, levels = c(FALSE, TRUE)), deparse.level = 0)
if (isLD) {
xSNPlistRandom <- intersect(SNPlist, unique(chrLD[which(names(chrLD)%in%eSNPlistRandom)]))
xSNPenrichStats <- table(assoc, factor(SNPlist%in%xSNPlistRandom, levels = c(FALSE, TRUE)), deparse.level = 0)
xTmp <- c(c(xSNPenrichStats), .enrichmentRatio(xSNPenrichStats))
} else {
xTmp <- NULL
}
eTmp <- c(c(eSNPenrichStats), .enrichmentRatio(eSNPenrichStats))
return(list(eSNP = eTmp, xSNP = xTmp))
})
eSNPsample <- rbind(eSNPsample, do.call("rbind", lapply(resParallel, function (l) {l$eSNP})))
eMeanEnrichRatio <- sum(eSNPsample[, 5])/length(eSNPsample[, 5])
Ze <- (eEnrichRatio-eMeanEnrichRatio)/sqrt((sum((eSNPsample[, 5]-eMeanEnrichRatio)^2))/(length(eSNPsample[, 5])-1))
if (isLD) {
xSNPsample <- rbind(xSNPsample, do.call("rbind", lapply(resParallel, function (l) {l$xSNP})))
xMeanEnrichRatio <- sum(xSNPsample[, 5])/length(xSNPsample[, 5])
Zx <- (xEnrichRatio-xMeanEnrichRatio)/sqrt((sum((xSNPsample[, 5]-xMeanEnrichRatio)^2))/(length(xSNPsample[, 5])-1))
} else {
Zx <- 1
}
iSample <- nSampleMin
cat(iSample, ".. ", sep = "")
catStep <- (nSample/10)
rm(resParallel)
while (iSample<nSample) {
resParallel <- mclapply2(X = seq(nSampleMin), mc.cores = mc.cores, FUN = function (i) {
eSNPlistRandom <- unlist(sapply(seq(nPool), function (g) {sample(popSNP4Sample[[g]], min(eSNPlistPool[g], length(popSNP4Sample[[g]])))}))
eSNPenrichStats <- table(assoc, factor(SNPlist%in%eSNPlistRandom, levels = c(FALSE, TRUE)), deparse.level = 0)
if (isLD) {
xSNPlistRandom <- intersect(SNPlist, unique(chrLD[which(names(chrLD)%in%eSNPlistRandom)]))
xSNPenrichStats <- table(assoc, factor(SNPlist%in%xSNPlistRandom, levels = c(FALSE, TRUE)), deparse.level = 0)
xTmp <- c(c(xSNPenrichStats), .enrichmentRatio(xSNPenrichStats))
} else {
xTmp <- NULL
}
eTmp <- c(c(eSNPenrichStats), .enrichmentRatio(eSNPenrichStats))
return(list(eSNP = eTmp, xSNP = xTmp))
})
eSNPsample <- rbind(eSNPsample, do.call("rbind", lapply(resParallel, function (l) {l$eSNP})))
eMeanEnrichRatio <- sum(eSNPsample[, 5])/length(eSNPsample[, 5])
Ze <- (eEnrichRatio-eMeanEnrichRatio)/sqrt((sum((eSNPsample[, 5]-eMeanEnrichRatio)^2))/(length(eSNPsample[, 5])-1))
if (isLD) {
xSNPsample <- rbind(xSNPsample, do.call("rbind", lapply(resParallel, function (l) {l$xSNP})))
xMeanEnrichRatio <- sum(xSNPsample[, 5])/length(xSNPsample[, 5])
Zx <- (xEnrichRatio-xMeanEnrichRatio)/sqrt((sum((xSNPsample[, 5]-xMeanEnrichRatio)^2))/(length(xSNPsample[, 5])-1))
} else {
Zx <- 1
}
iSample <- iSample + nSampleMin
if ((iSample-(iSample%/%catStep*catStep)) == 0) {cat(iSample, ".. ", sep = "")}
rm(resParallel)
}
if (empiricPvalue) {
whichPvalue <- 2
} else {
whichPvalue <- 1
}
empiricPvalue.eSNP <- sum(eEnrichRatio<eSNPsample[, 5])/length(eSNPsample[, 5])
statisticPvalue.eSNP <- pnorm(Ze, lower.tail = FALSE)
object@eSNP <- enrichSNP(List = eSNPlist, Table = eEnrichment, EnrichmentRatio = eEnrichRatio, Z = Ze, PValue = c(Distribution = statisticPvalue.eSNP, Empirical = empiricPvalue.eSNP)[whichPvalue], Resampling = eSNPsample)
if (isLD) {
empiricPvalue.xSNP <- sum(xEnrichRatio<xSNPsample[, 5])/length(xSNPsample[, 5])
statisticPvalue.xSNP <- pnorm(Zx, lower.tail = FALSE)
object@xSNP <- enrichSNP(List = xSNPlist, Table = xEnrichment, EnrichmentRatio = xEnrichRatio, Z = Zx, PValue = c(Distribution = statisticPvalue.xSNP, Empirical = empiricPvalue.xSNP)[whichPvalue], Resampling = xSNPsample)
} else {
object@xSNP <- enrichSNP()
}
return(object)
}
.compareEnrich <- function (object1, object2, nSample, empiricPvalue, sigThresh, MAFpool, mc.cores) {
DATA <- object1["Data"]
chrLD <- object1["LD"]
isLD <- length(chrLD) != 0
eEnrichment <- object1@eSNP@Table
eEnrichRatio <- object1@eSNP@EnrichmentRatio
eSNPlist <- object1@eSNP@List
eList <- union(object1@eSNP@List, object2@eSNP@List)
if (isLD) {
xList <- union(object1@xSNP@List, object2@xSNP@List)
data <- DATA[DATA[, "SNP"] %in% union(eList, xList), ]
} else {
data <- DATA[DATA[, "SNP"] %in% eList, ]
}
data[, "MAFpool"] <- NA
data[, "MAFpool"] <- as.factor(cut(data[, "MAF"], breaks = MAFpool, labels = FALSE, include.lowest = TRUE))
nPool <- nlevels(data$MAFpool)
eSNPlistPool <- table(data[eSNPlist, "MAFpool"])
eSNPsample <- NULL
nSampleMin <- min(nSample, max(1000, nSample*10/100))
if (isLD) {
xEnrichment <- object1@xSNP@Table
xEnrichRatio <- object1@xSNP@EnrichmentRatio
xSNPlist <- object1@xSNP@List
xSNPlistPool <- table(data[xSNPlist, "MAFpool"])
xSNPsample <- NULL
} else {}
popSNP4Sample <- split(eList, data[eList, "MAFpool"])
assoc_eSNP <- factor(data[eList, "PVALUE"] < sigThresh, levels = c(FALSE, TRUE))
if (isLD) {
assoc_xSNP <- factor(data[xList, "PVALUE"] < sigThresh, levels = c(FALSE, TRUE))
} else { }
rm(data)
cat("0.. ")
resParallel <- mclapply2(X = seq(nSampleMin), mc.cores = mc.cores, FUN = function (i) {
eSNPlistRandom <- unlist(sapply(seq(nPool), function (g) {sample(popSNP4Sample[[g]], min(eSNPlistPool[g], length(popSNP4Sample[[g]])))}))
eSNPenrichStats <- table(assoc_eSNP, factor(eList%in%eSNPlistRandom, levels = c(FALSE, TRUE)))
eTmp <- c(c(eSNPenrichStats), .enrichmentRatio(eSNPenrichStats))
if (isLD) {
xSNPlistRandom <- intersect(xList, unique(chrLD[which(names(chrLD)%in%eSNPlistRandom)]))
xSNPenrichStats <- table(assoc_xSNP, factor(xList%in%xSNPlistRandom, levels = c(FALSE, TRUE)))
xTmp <- c(c(xSNPenrichStats), .enrichmentRatio(xSNPenrichStats))
} else {
xTmp <- NULL
}
return(list(eSNP = eTmp, xSNP = xTmp))
})
eSNPsample <- rbind(eSNPsample, do.call("rbind", lapply(resParallel, function (l) {l$eSNP})))
eSNPanyDup <- eSNPsample[!duplicated(eSNPsample), ]
if (is.matrix(eSNPanyDup)) {
if (length(which(is.na(eSNPsample[, 5]))) != nrow(eSNPsample) & length(which(is.na(eSNPsample[, 5]))) > 0) {
eSNPsample[is.na(eSNPsample[, 5]), 5] <- rep(mean(eSNPsample[(!is.na(eSNPsample[, 5]) & is.finite(eSNPsample[, 5])), 5]), length(eSNPsample[is.na(eSNPsample[, 5]), 5]))
} else {}
if (length(which(is.infinite(eSNPsample[, 5]))) != nrow(eSNPsample) & length(which(is.infinite(eSNPsample[, 5]))) > 0) {
eSNPsample[is.infinite(eSNPsample[, 5]), 5] <- rep(max(eSNPsample[!is.infinite(eSNPsample[, 5]), 5])*1.05, length(eSNPsample[is.infinite(eSNPsample[, 5]), 5]))
} else {}
eMeanEnrichRatio <- sum(eSNPsample[, 5])/length(eSNPsample[, 5])
Ze <- (eEnrichRatio-eMeanEnrichRatio)/sqrt((sum((eSNPsample[, 5]-eMeanEnrichRatio)^2))/(length(eSNPsample[, 5])-1))
} else {
Ze <- as.numeric(NA)
}
if (isLD) {
xSNPsample <- rbind(xSNPsample, do.call("rbind", lapply(resParallel, function (l) {l$xSNP})))
xSNPanyDup <- xSNPsample[!duplicated(xSNPsample), ]
if (is.matrix(xSNPanyDup)) {
if (length(which(is.na(xSNPsample[, 5]))) != nrow(xSNPsample) & length(which(is.na(xSNPsample[, 5]))) > 0) {
xSNPsample[is.na(xSNPsample[, 5]), 5] <- rep(mean(xSNPsample[(!is.na(xSNPsample[, 5]) & is.finite(xSNPsample[, 5])), 5]), length(xSNPsample[is.na(xSNPsample[, 5]), 5]))
} else {}
if (length(which(is.infinite(xSNPsample[, 5]))) != nrow(xSNPsample) & length(which(is.infinite(xSNPsample[, 5]))) > 0) {
xSNPsample[is.infinite(xSNPsample[, 5]), 5] <- rep(max(xSNPsample[!is.infinite(xSNPsample[, 5]), 5])*1.05, length(xSNPsample[is.infinite(xSNPsample[, 5]), 5]))
} else {}
xMeanEnrichRatio <- sum(xSNPsample[, 5])/length(xSNPsample[, 5])
Zx <- (xEnrichRatio-xMeanEnrichRatio)/sqrt((sum((xSNPsample[, 5]-xMeanEnrichRatio)^2))/(length(xSNPsample[, 5])-1))
} else {
Zx <- as.numeric(NA)
}
} else {
Zx <- as.numeric(NA)
}
iSample <- nSampleMin
cat(iSample, ".. ", sep = "")
catStep <- (nSample/10)
rm(resParallel)
while (iSample<nSample) {
resParallel <- mclapply2(X = seq(nSampleMin), mc.cores = mc.cores, FUN = function (i) {
eSNPlistRandom <- unlist(sapply(seq(nPool), function (g) {sample(popSNP4Sample[[g]], min(eSNPlistPool[g], length(popSNP4Sample[[g]])))}))
eSNPenrichStats <- table(assoc_eSNP, factor(eList%in%eSNPlistRandom, levels = c(FALSE, TRUE)))
eTmp <- c(c(eSNPenrichStats), .enrichmentRatio(eSNPenrichStats))
if (isLD) {
xSNPlistRandom <- intersect(xList, unique(chrLD[which(names(chrLD)%in%eSNPlistRandom)]))
xSNPenrichStats <- table(assoc_xSNP, factor(xList%in%xSNPlistRandom, levels = c(FALSE, TRUE)))
xTmp <- c(c(xSNPenrichStats), .enrichmentRatio(xSNPenrichStats))
} else {
xTmp <- NULL
}
return(list(eSNP = eTmp, xSNP = xTmp))
})
eSNPsample <- rbind(eSNPsample, do.call("rbind", lapply(resParallel, function (l) {l$eSNP})))
eSNPanyDup <- eSNPsample[!duplicated(eSNPsample), ]
if (is.matrix(eSNPanyDup)) {
if (length(which(is.na(eSNPsample[, 5]))) != nrow(eSNPsample) & length(which(is.na(eSNPsample[, 5]))) > 0) {
eSNPsample[is.na(eSNPsample[, 5]), 5] <- rep(mean(eSNPsample[(!is.na(eSNPsample[, 5]) & is.finite(eSNPsample[, 5])), 5]), length(eSNPsample[is.na(eSNPsample[, 5]), 5]))
} else {}
if (length(which(is.infinite(eSNPsample[, 5]))) != nrow(eSNPsample) & length(which(is.infinite(eSNPsample[, 5]))) > 0) {
eSNPsample[is.infinite(eSNPsample[, 5]), 5] <- rep(max(eSNPsample[is.finite(eSNPsample[, 5]), 5])*1.05, length(eSNPsample[is.infinite(eSNPsample[, 5]), 5]))
} else {}
eMeanEnrichRatio <- sum(eSNPsample[, 5])/length(eSNPsample[, 5])
Ze <- (eEnrichRatio-eMeanEnrichRatio)/sqrt((sum((eSNPsample[, 5]-eMeanEnrichRatio)^2))/(length(eSNPsample[, 5])-1))
} else {
Ze <- as.numeric(NA)
}
if (isLD) {
xSNPsample <- rbind(xSNPsample, do.call("rbind", lapply(resParallel, function (l) {l$xSNP})))
xSNPanyDup <- xSNPsample[!duplicated(xSNPsample), ]
if (is.matrix(xSNPanyDup)) {
if (length(which(is.na(xSNPsample[, 5]))) != nrow(xSNPsample) & length(which(is.na(xSNPsample[, 5]))) > 0) {
xSNPsample[is.na(xSNPsample[, 5]), 5] <- rep(mean(xSNPsample[(!is.na(xSNPsample[, 5]) & is.finite(xSNPsample[, 5])), 5]), length(xSNPsample[is.na(xSNPsample[, 5]), 5]))
} else {}
if (length(which(is.infinite(xSNPsample[, 5]))) != nrow(xSNPsample) & length(which(is.infinite(xSNPsample[, 5]))) > 0) {
xSNPsample[is.infinite(xSNPsample[, 5]), 5] <- rep(max(xSNPsample[is.finite(xSNPsample[, 5]), 5])*1.05, length(xSNPsample[is.infinite(xSNPsample[, 5]), 5]))
} else {}
xMeanEnrichRatio <- sum(xSNPsample[, 5])/length(xSNPsample[, 5])
Zx <- (xEnrichRatio-xMeanEnrichRatio)/sqrt((sum((xSNPsample[, 5]-xMeanEnrichRatio)^2))/(length(xSNPsample[, 5])-1))
} else {
Zx <- as.numeric(NA)
}
} else {
Zx <- as.numeric(NA)
}
iSample <- iSample + nSampleMin
if ((iSample-(iSample%/%catStep*catStep)) == 0) {cat(iSample, ".. ", sep = "")}
rm(resParallel)
}
if (empiricPvalue) {
whichPvalue <- 2
} else {
whichPvalue <- 1
}
empiricPvalue.eSNP <- sum(eEnrichRatio<eSNPsample[, 5])/length(eSNPsample[, 5])
statisticPvalue.eSNP <- pnorm(Ze, lower.tail = FALSE)
object1@eSNP <- enrichSNP(List = eSNPlist, Table = eEnrichment, EnrichmentRatio = eEnrichRatio, Z = Ze, PValue = c(Distribution = statisticPvalue.eSNP, Empirical = empiricPvalue.eSNP)[whichPvalue], Resampling = eSNPsample)
if (isLD) {
empiricPvalue.xSNP <- sum(xEnrichRatio<xSNPsample[, 5])/length(xSNPsample[, 5])
statisticPvalue.xSNP <- pnorm(Zx, lower.tail = FALSE)
object1@xSNP <- enrichSNP(List = xSNPlist, Table = xEnrichment, EnrichmentRatio = xEnrichRatio, Z = Zx, PValue = c(Distribution = statisticPvalue.xSNP, Empirical = empiricPvalue.xSNP)[whichPvalue], Resampling = xSNPsample)
} else {
object1@xSNP <- enrichSNP()
}
if (is.enrichment(object1)) {
object1@Chromosomes <- lapply(object1@Chromosomes, reset, "Z")
object1@Chromosomes <- lapply(object1@Chromosomes, reset, "PValue")
object1@Chromosomes <- lapply(object1@Chromosomes, reset, "Resampling")
} else {}
return(object1)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.