args <- commandArgs(trailingOnly = TRUE)

if (length(args) != 6) {
   cat("*****
Converter of PLINK to OmicABEL format. This converter is intended 
for use with small example data sets, for large data sets please 
use other procedures (e.g. direct manipulations of input using the 
DatABEL library)
*****
usage: R --slave -f convertGenA2OmicA.R --args tpedfile tfamfile phenofile \
                              covFrom-covTo phenoFrom-phenoTo outfilebase

tpedfile:        PLINK TPED file
tfamfile:        PLINK (t)fam-file
phenofile:       phenotype file (following GenA format)
covFrom-covTo:   which covariates to include, x1-x2 where x1 is the index 
                 of the first column of covariates in the phenotype file 
                 and x2 is the index of the last one (all covariates in 
                 between will be included). Use 0-0 for no covariates
pheFrom-pheTo:   which phenotypes to analyze, x1-x2 where x1 is the index 
                 of the first column of phenotypes for analysis  
                 and x2 is the index of the last one (all covariates in 
                 between will be included). Use 0-0 for no covariates
outfilebase:     the 'base' name of the converted files
   ")
   stop("command line arguments are not valid")
}

inPed <- args[1]
inFam <- args[2]
inPhe <- args[3]
XLFromTo <- strsplit(args[4],"-")[[1]]
if (length(XLFromTo) != 2) {
   stop("covFrom-covTo argument should be of form i1-i2 where i1 and i2 are integers representing the range of covariates in the phenofile (use 0-0 for no covariates)")
}
cat("Covariates:",XLFromTo[1],"to",XLFromTo[2],"\n")
pheFromTo <- strsplit(args[5],"-")[[1]]
if (length(pheFromTo) != 2) {
   stop("phenoFrom-phenoTo argument should be of form i1-i2 where i1 and i2 are integers representing the range of phenotypes in the phenofile")
}
cat("Phenotypes:",pheFromTo[1],"to",pheFromTo[2],"\n")
outRaw <- paste(args[6],".raw",sep="")
outRData <- paste(args[6],".RData",sep="")
outRel <- paste(args[6],".grel",sep="")
outDAphe <- paste(args[6],".phe",sep="")
outXL <- paste(args[6],".XL",sep="")
outXR <- paste(args[6],".XR",sep="")

library(GenABEL)
#if (! require( library(GenABEL) ) ) 
#   stop("GenABEL required")

library(DatABEL)
#if (! require( library(DatABEL) ) ) 
#   stop("DatABEL required")

unlink(paste(args[6],".*",sep=""))

# arrange GenA data
convert.snp.tped(tpedfile=inPed,tfamfile=inFam,outfile=outRaw)
df <- load.gwaa.data(genofile=outRaw,phenofile=inPhe)
# compute gkin
gkin <- ibs(df,weigh="freq")
# save RData
save(df, gkin, file = outRData)
# arrange grel
grel4OmicA <- gkin
grel4OmicA[upper.tri(grel4OmicA)] <- t(gkin)[upper.tri(gkin)]
grel4OmicA <- grel4OmicA * 2
# arrange XL
  XL <- matrix(rep(1,nids(df)),ncol=1)
  rownames(XL) <- idnames(df)
  colnames(XL) <- "const"
if (XLFromTo[1]>0) {
  XL <- cbind(XL,as.numeric(phdata(df)[,XLFromTo[1]:XLFromTo[2]]))
  colnames(XL) <- c("const",names(phdata(df))[XLFromTo[1]:XLFromTo[2]])
}
# arrange phes
phes <- as.matrix(phdata(df)[, pheFromTo[1]:pheFromTo[2] ] )

# export data to OmicA format
tmp <- matrix2databel(from=grel4OmicA,filename=outRel,type="DOUBLE")
tmp <- matrix2databel(from=as.numeric(gtdata(df)),filename=outXR,type="DOUBLE")
tmp <- matrix2databel(from=XL,filename=outXL,type="DOUBLE")
tmp <- matrix2databel(from=phes,filename=outDAphe,type="DOUBLE")

cat("*** You can now run CLAK-GWAS with command (replace capitals with your values): ***
      CLAK-GWAS -var eigen -nths NUMTHREADS -cov",outXL,"-phi",outRel," \
      		    -snp",outXR,"-pheno",outDAphe,"-out MYRESFILE\n")