1. 程式人生 > >昨天寫的關於處理LBLGXE的一段R程式

昨天寫的關於處理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)