guess genotypes from allele dosages

guess_geno.R

# guess genotypes from allele dosages
#
# alleleprob = n_ind x n_alleles matrix of allele dosages (rows sum to 1)
# epsilon = tolerance value

guess_geno_from_alleleprob <-
    function(alleleprob, epsilon=0.01)
    {
        stopifnot(is.matrix(alleleprob)) # expect ind x allele matrix
        stopifnot(is.numeric(epsilon), length(epsilon)==1, epsilon >= 0, epsilon <= 1)

        # allele codes
        alleles <- colnames(alleleprob)
        if(is.null(alleles)) alleles <- LETTERS[seq_len(ncol(alleleprob))]

        top2pr <- t(apply(alleleprob, 1, function(a) sort(a, decreasing=TRUE)[1:2]))
        top2alle <- t(apply(alleleprob, 1, function(a) order(a, decreasing=TRUE)[1:2]))

        guess <- setNames(rep(NA, nrow(alleleprob)), rownames(alleleprob))

        hom <- top2pr[,1] > 1-epsilon
        if(any(hom)) {
            a <- alleles[top2alle[hom,1]]
            guess[hom] <- paste0(a, a)
        }

        het <- top2pr[,2] > (1-epsilon)/2
        if(any(het)) {
            ahet <- cbind(top2alle[het,1], top2alle[het,2])
            ahet <- t(apply(ahet, 1, sort))
            a <- alleles[ahet[,1]]
            b <- alleles[ahet[,2]]
            guess[het] <- paste0(a, b)
        }

        guess

    }


### Example
#
# library(qtl2)
#
# tmpfile <- tempfile()
# file <- paste0("https://raw.githubusercontent.com/rqtl/",
#                "qtl2data/master/DOex/DOex_alleleprobs.rds")
# download.file(file, tmpfile)
# apr <- readRDS(tmpfile)
#
# p <- pull_genoprobpos(apr, marker="UNC020114284")
# g <- guess_geno_from_alleleprob(p)
# mean(is.na(g)) # 49% missing
# table(g)