parallel_varcomp.R

get_varcomp = function(x) {
    library(hierfstat)
    loci = data.frame(x)
    res <- varcomp(cbind(levels, loci),diploid=T)$overall
}

finish_varcomp = function(m) {
    tot <- apply(m, 2, sum, na.rm = TRUE)
    nblevels <- length(tot)
    f <- matrix(rep(0, (nblevels - 1)^2), ncol = (nblevels -
                                                      1))
    for (i in 1:(nblevels - 1)) {
        for (j in i:(nblevels - 1)) {
            f[i, j] <- sum(tot[i:j])/sum(tot[i:nblevels])
        }
    }
    row.names(m) <- lnames
    print(names(tot))
    tf <- t(f)
    row.names(tf) <- fnames
    f <- t(tf)
    row.names(f) <- c("Total", fnames[-length(fnames)])
    return(list(loc = m, overall = tot, F = f))
}


library(hierfstat)
library(data.table)
library(snow)
data = data.frame(fread("hierfstat_samtools_imputed.txt", header=T, sep="\t"))
levels = data.frame(data[,1:2])
loci = data[,3:ncol(data)]
lnames=names(loci)
fnames=c(names(levels), "Ind")
cl = makeSOCKcluster(50)
clusterExport(cl, "levels", envir=environment())
res = matrix(parCapply(cl, loci, get_varcomp), nrow=length(names(loci)),byrow=T)
res = finish_varcomp(res)
saveRDS(res, "hierfstat_samtools_imputed.rds")
stopCluster(cl)
print("Done!")