R/ScoreBL.R

In wenlongren/ScoreBL: ScoreBL:a fast score test integrated with Bayesian Lasso for multi-locus GWAS

ScoreBLtest <- function(genfile,famfile,bimfile,IniPthreshold,LOD){
  XX.gen <- as.matrix(fread(file = genfile, sep = "\t",header = FALSE))
  XX.fam <- fread(file = famfile, sep = "\t",header = TRUE)
  XX.bim <- fread(file = bimfile, sep = "\t",header = TRUE)
  XX.gen <- t(XX.gen)
  x <- as.bed.matrix(XX.gen, XX.fam, XX.bim)
  standardize(x) <- "p"
  Yphe <- as.matrix(x@ped[,6])
  Xfix <- matrix(1,dim(x)[1],1)
  KM <- (1/dim(XX.gen)[2])*GRM(x)
  nmaker <- dim(XX.gen)[2]
  retainP <- 0.5/nmaker
  LODthred <- LOD
  iniP <- IniPthreshold
  resLODTotal2 <- numeric()
  resEach <- numeric()
  resEtotal <- numeric()
  resLodTot2 <- numeric()
  resScore <- association.test(x = x,Y = Yphe,X = Xfix, method = "lmm", response = "quantitative", test = "score", K = KM)
  initialChose <- resScore[which(resScore[,8]<=iniP),]
  if(length(initialChose) > 0){
    retainSNP <-  resScore[which(resScore[,8]<=retainP),]
    if(length(retainSNP) > 0){
      diffSNP <- setdiff(initialChose[,3],retainSNP[,3])
      initialChose <- initialChose[match(diffSNP,initialChose[,3]),]
      retainGen <- as.matrix(XX.gen[,retainSNP[,3]])
      bb1 <- ebayes_EM(Xfix,retainGen,as.matrix(Yphe))
      lod1 <- likelihood(Xfix,retainGen,as.matrix(Yphe),bb1)
      retainSave <- cbind(retainSNP,lod1,bb1)
      retainSave <- retainSave[which(retainSave[,9]>=LODthred),]
      colnames(retainSave) <- c("chr","pos","id","A1","A2","freq","score","p","lod","lrtb")
      choseGene <- XX.gen[,diffSNP]
      bb2 <- ebayes_EM(Xfix,choseGene,as.matrix(Yphe))
      lod2 <- likelihood(Xfix,choseGene,as.matrix(Yphe),bb2)
      resEach <- cbind(initialChose,lod2,bb2)
      resEach <-resEach[which(resEach[,9]>=LODthred),]
      colnames(resEach) <- c("chr","pos","id","A1","A2","freq","score","p","lod","lrtb")
      resEtotal <- rbind(retainSave,resEach)
    }else{
      choseGene <- XX.gen[,initialChose[,3]]
      bb <- ebayes_EM(Xfix,choseGene,as.matrix(Yphe))
      lod <- likelihood(Xfix,choseGene,as.matrix(Yphe),bb)
      resEach <- cbind(initialChose,lod,bb)
      resEach <-resEach[which(resEach[,9]>=LODthred),]
      colnames(resEach) <- c("chr","pos","id","lsb1.0","lsb3.0","freq","score","p","lod","lrtb")
      resEtotal <- rbind(resEtotal,resEach)
    }
  }
  return (resEtotal)
}