昨天写的关于处理LBLGXE的一段R程序

LBLprocess <- function(inputFile, rootDir){
dt <- read.table(inputFile, colClasses = "character")
cov <- sample(c(0, 1), dim(dt)[1], replace = T)
vec <- matrix(0,  dim(dt)[1], 202)
for(i in 1:dim(dt)[1]){
x <- strsplit(dt[i,1], "")
for(j in 1:(length(x[[1]])-1)){
vec[i, 2*j+1] <- as.numeric(x[[1]][j])
if(vec[i, 2*j+1] == 1){
prob <- runif(1, 0.0, 1.0)
if(prob >= 0.8) vec[i, 2*(j+1)] <- 1
if(prob < 0.8) vec[i, 2*(j+1)] <- 0
}
}
if(x[[1]][length(x[[1]])] == 'A'){
vec[i, 1] = 0
}
else{
vec[i, 1] = 1
}
}
vec[,2] <- cov

dvec <- as.data.frame(vec)
#save(dvec, file = "LBL/LBLdata.rda")
for(i in 1:20){
col <- c(1, 2, ((i-1)*10+3):(i*10+2))
head <- c("affected", "cov", "M1.1", "M1.2", "M2.1", "M2.2", "M3.1", "M3.2", "M4.1", "M4.2", "M5.1", "M5.2")
LBL.data <- dvec[, col]
path <- paste(rootDir, "LBLdata", i, ".rda", sep = "")
#LBL.data <- as.data.frame(LBLmat)
names(LBL.data) <- head
save(LBL.data, file = path)
}
}

comb <- function(str){
if(length(str) == 0){
res = "00000"
return(res)
}
if(length(str) == 1){
return(str)
}
else{
ch <- strsplit(str[1], "")
for(i in 2:length(str)){
for(j in 1:5){
x <- strsplit(str[i], "")
if(x[[1]][j] == '1' | ch[[1]][j] == 1){
ch[[1]][j] = '1'
}
}
}
s <- paste(ch[[1]][1], ch[[1]][2], ch[[1]][3], ch[[1]][4], ch[[1]][5], sep = "")
return(s)
}
}

LBLprocess("LBL/00000000001", "LBL/")
library(hapassoc)
library(dummies)
library(LBLGXE)
final  <- ""
for(i in 1:20){
path <- paste("LBL/LBLdata", i, ".rda", sep = "")
load(path)
out.LBL<-LBL(LBL.data, numSNPs=5, burn.in=100, num.it=1000)
bf <- out.LBL$BF
s <- names(bf)
zo <- c()
for(j in 1:length(bf)){
if(nchar(s[j]) == 6){
if(nchar(bf[j]) == 4){
next

}
if(as.numeric(bf[j]) >= 2.0){
zo[length(zo)+1] <- substr(s[j], 2, 6)
}
}
}
s <- comb(zo)
final <- paste(final, s, sep = "")
}

write.table(final, "LBL/LBLfinal.txt", row.names = F, col.names = F, quote = F)