collr2: QCDH

require(devtools)
install_bitbucket("kindlychung/txtutils")
# install_bitbucket("kindlychung/manqq2")
# install_github("kindlychung/gmailR")
# install_bitbucket("kindlychung/autoGmail")
install_bitbucket("kindlychung/collr2")


library(collr2)
setwd("~/Desktop/mmp3")
plinkArgs    = getPlinkParam(allow_no_sex = "", pheno = "mmp13.phe", pheno_name = "Page", covar = "mmp13.phe", covar_name = "Sex,Cage", linear = "hide-covar")
initGwasArgs = getPlinkParam(maf=.01, hwe=1e-4, allow_no_sex = "", pheno = "mmp13.phe", pheno_name = "Page", assoc= "")
crt1("~/Desktop/mmp3", "test4", plinkArgs, initGwasArgs, FALSE, 0.1, 5)

debug(permute)
pv = permute(plinkArgs = plinkArgs, initGwasArgs = initGwasArgs, phenoFileOrig = "mmp13.phe", pFilter = .2, nShiftMax = 10, n=10)
sort(pv)



require(collr2)
setwd("/media/data1/kaiyin/RS123_1KG/")
plinkArgs = getPlinkParam(
    allow_no_sex = "", missing_phenotype = 9999,
    pheno = "RS123.1kg.pheno/dermatology.csv",
    pheno_name = "EasilyBurn",
    covar = "RS123.1kg.pheno/dermatology.csv",
    covar_name = "sex,age,SC,SunProtect",
    logistic = "hide-covar",
    one = ""
)
initGwasArgs = getPlinkParam(
    allow_no_sex = "",
    missing_phenotype = 9999,
    pheno = "RS123.1kg.pheno/dermatology.csv",
    pheno_name = "EasilyBurn",
    one = "",
    assoc=""
)
crt1(
    wDir = "/media/data1/kaiyin/RS123_1KG",
    taskName = "easilyburn2_s200_p0.01",
    plinkArgs = plinkArgs,
    initGwasArgs = initGwasArgs,
    initGwas = TRUE,
    pFilter = 1e-2,
    nMaxShift = 200
)

require(collr2)
setwd("/media/data1/kaiyin/RS123_1KG")
plinkArgs = getPlinkParam(
    allow_no_sex = "", missing_phenotype = 9999,
    pheno = "RS123.1kg.pheno/dermatology.csv",
    pheno_name = "AK",
    covar = "RS123.1kg.pheno/dermatology.csv",
    covar_name = "sex,age,SC",
    linear = "hide-covar"
)
initGwasArgs = getPlinkParam(
    allow_no_sex = "", missing_phenotype = 9999,
    pheno = "RS123.1kg.pheno/dermatology.csv",
    pheno_name = "AK",
    assoc = ""
)
crt1(
    wDir = "/media/data1/kaiyin/RS123_1KG",
    taskName = "AK_s200_p0.01",
    plinkArgs = plinkArgs,
    initGwasArgs = initGwasArgs,
    initGwas = TRUE,
    pFilter = 1e-2,
    nMaxShift = 200
)


#############################
## sagging
#############################
require(collr2)
setwd("/media/data1/NL_AU_UK/NL_AU_UK_sagging_replicate/chr5frag/")
plinkArgs = getPlinkParam(
    allow_no_sex = "",
    linear = "hide-covar"
)
initGwasArgs = getPlinkParam(
    allow_no_sex = "",
    assoc=""
)
routine2(
    pheno = "NL_AU_UK_sagging1.pheno",
    pheno_name = "sag",
    covar_name = "sex,age",
    plinkArgs = plinkArgs,
    initGwasArgs = initGwasArgs,
    pFilter = 1e-1,
    nMaxShift = 400
)

require(collr2)
setwd("/media/data1/kaiyin/RS123_1KG")
plinkArgs = getPlinkParam(
    allow_no_sex = "", missing_phenotype = 9999,
    pheno = "RS123.1kg.pheno/dermatology.csv",
    pheno_name = "Sagging",
    covar = "RS123.1kg.pheno/dermatology.csv",
    covar_name = "sex,age,SC",
    linear = "hide-covar"
)
initGwasArgs = getPlinkParam(
    allow_no_sex = "", missing_phenotype = 9999,
    pheno = "RS123.1kg.pheno/dermatology.csv",
    pheno_name = "Sagging",
    assoc = ""
)
crt1(
    wDir = "/media/data1/kaiyin/RS123_1KG",
    taskName = "Sagging_s200_p0.01",
    plinkArgs = plinkArgs,
    initGwasArgs = initGwasArgs,
    initGwas = TRUE,
    pFilter = 1e-2,
    nMaxShift = 200
)


require(collr2)
setwd("/media/data1/kaiyin/RS123_1KG")
plinkArgs = getPlinkParam(
    allow_no_sex = "", missing_phenotype = 9999,
    pheno = "RS123.1kg.pheno/dermatology.csv",
    pheno_name = "Log_PS",
    covar = "RS123.1kg.pheno/dermatology.csv",
    covar_name = "sex,age,SC",
    linear = "hide-covar"
)
initGwasArgs = getPlinkParam(
    allow_no_sex = "", missing_phenotype = 9999,
    pheno = "RS123.1kg.pheno/dermatology.csv",
    pheno_name = "Log_PS",
    assoc = ""
)
crt1(
    wDir = "/media/data1/kaiyin/RS123_1KG",
    taskName = "Log_PS_s200_p0.01",
    plinkArgs = plinkArgs,
    initGwasArgs = initGwasArgs,
    initGwas = TRUE,
    pFilter = 1e-2,
    nMaxShift = 200
)






require(collr2)
setwd("/media/data1/kaiyin/RS123_1KG/")
globalMinVec = numeric(0)
phenoFileOrig = "RS123.1kg.pheno/dermatology.csv"
phenoFile = "~/Desktop/shuffle_pheno.csv"
for(i in 1:100) {
    i = 1
    taskName = sprintf("easilyburn2_permutation_%d", i)
    shuffleLines(phenoFileOrig, phenoFile, "EasilyBurn")
    plinkArgs = getPlinkParam(allow_no_sex = "", missing_phenotype = 9999, pheno = phenoFile, covar = phenoFile, covar_name = "sex,age,SC,SunProtect", logistic = "hide-covar", one = "", pheno_name = "EasilyBurn")
    initGwasArgs = getPlinkParam(allow_no_sex = "", missing_phenotype = 9999, pheno = phenoFile, one = "", pheno_name = "EasilyBurn", assoc="", maf = 0.01, hwe = 1e-4, geno=0.02)
    hubcollr = collrinfo()
    bedcollinfo(hubcollr, genbed=FALSE)
    taskinfo(hubcollr, taskName, plinkArgs, initGwasArgs, TRUE, "bin")
    taskBedsPlinkOut(hubcollr, taskName, hubcollr[[taskName]]$fullGwasOut, 1e-3, 1, 200)
    taskAnalyze(hubcollr, taskName)
    readcoll.task(hubcollr, taskName)
    globalMin = min(hubcollr[[taskName]]$minPvals)
    if(globalMin > 1)      globalMin = 1
    if(globalMin < 1e-300) globalMin = 1e-300
    globalMinVec = c(globalMinVec, globalMin)
    message("*****************************************************")
    print(globalMinVec)
    message("*****************************************************")
    unlink(hubcollr[[taskName]]$taskPath, recursive = TRUE, force = TRUE)
}
perm




# erf4 huid07, winszie=200
setwd("/media/data1/kaiyin/erf4/")
require(collr2)
plinkArgs = getPlinkParam(allow_no_sex = "", pheno = "erf4.pheno/erf_all_phdata.txt", pheno_name = "huid07", covar = "erf4.pheno/erf_all_phdata.txt", covar_name = "sex,age", linear = "hide-covar")
initGwasArgs = getPlinkParam(maf=.01, hwe=1e-4, allow_no_sex = "", pheno = "erf4.pheno/erf_all_phdata.txt", pheno_name = "huid07", assoc= "")
crt1(".", "huid07b", plinkArgs, initGwasArgs, TRUE, "con", 0.001, 200)


shuffleLines = function(phenoFileOrig, phenoFile, shuffleVar) {
    tab = read.table(phenoFileOrig, head=TRUE)
    # print(head(tab, 10))
    tab[, shuffleVar] = sample(tab[, shuffleVar])
    # print(head(tab, 10))
    write.table(tab, phenoFile, quote = FALSE, row.names = FALSE)
}
permute = function(plinkArgs, initGwasArgs, phenoFileOrig, phenoFileShuffle = "~/Desktop/shuffle_pheno.csv", wDir=".", n=100) {
    setwd(wDir)
    globalMinVec = numeric(0)
    plinkArgs$pheno = phenoFileShuffle
    if(!is.null(plinkArgs$covar))  plinkArgs$covar = phenoFileShuffle
    initGwasArgs$pheno = phenoFileShuffle
    taskName = "permutation_analysis"
    for(i in 1:n) {
        shuffleLines(phenoFileOrig, phenoFileShuffle, plinkArgs$pheno_name)
        hubcollr = collrinfo()
        bedcollinfo(hubcollr, genbed=FALSE)
        taskinfo(hubcollr, taskName, plinkArgs, initGwasArgs, TRUE)
        taskBedsPlinkOut(hubcollr, taskName, hubcollr[[taskName]]$fullGwasOut, 1e-3, 1, 200)
        taskAnalyze(hubcollr, taskName)
        readcoll.task(hubcollr, taskName)
        globalMin = min(hubcollr[[taskName]]$minPvals)
        if(globalMin > 1)      globalMin = 1
        if(globalMin < 1e-300) globalMin = 1e-300
        globalMinVec = c(globalMinVec, globalMin)
        message("*****************************************************")
        print(globalMinVec)
        message("*****************************************************")
        unlink(hubcollr[[taskName]]$taskPath, recursive = TRUE, force = TRUE)
    }
    globalMinVec
}


require(collr2)
setwd("/media/data1/kaiyin/erf4")
phenoFileOrig = "erf4.pheno/erf_all_phdata.txt"


require(collr2)
plinkArgs = getPlinkParam(allow_no_sex = "", pheno = "erf4.pheno/erf_all_phdata.txt", pheno_name = "huid07", covar = "erf4.pheno/erf_all_phdata.txt", covar_name = "sex,age", linear = "hide-covar")
initGwasArgs = getPlinkParam(maf=.01, hwe=1e-4, allow_no_sex = "", pheno = "erf4.pheno/erf_all_phdata.txt", pheno_name = "huid07", assoc= "")