In vsbuffalo/ProgenyArray: Some statistical genetics for progeny array experiments

load_all("~/Projects/ProgenyArray")

library(rmarkdown)
library(Rcpp)
library(parallel)
library(caret)
library(gridExtra)
library(reshape2)
library(animation)
options(mc.cores=4)
set.seed(42)

Some Notes

Initialize based on k-means (k=2), and then shrink these probabilities.

(1) Remove all loci incomplete in all individuals early on.

(2) For E-step, we'll have to use all complete loci only. Note that for loci that are homozygous in the other parent, we can trivially impute progeny genotype to supply to this model -- this may be a valuable approach later on.

(3) $\pi_k$ for each of the $k \in {1, 2}$ haplotypes may be set to $1/2$ or let vary as a free parameter. As a free parameter, this could be a useful diagnostic, and potentially indicate when a region is homozygous in the shared parent. Responsibilities are also a valuable diagnostic.

EM

#load_all("~/Projects/ProgenyArray")
set.seed(0)

sf <- sibFamily(40, 1000, ehet=0.8, ehom=0.1)
pa <- ProgenyArray(sf$progeny, sf$genos, mothers=sf$metadata$mother)
pa <- calcFreqs(pa)
pa <- inferParents(pa, 0.8, 0.1)
rngs <- GRanges("1", IRanges(1:1000, width=1))
chrom <- GRanges("1", IRanges(1, 1000))
pa@ranges <- rngs

ph <- phaseSibFamily(pa, 1, "1", chrom)

# res <- emHap(pa, 1L, "parent_1", 1:nrow(progenyGenotypes(pa)), 0.8, 0.1)

Plots

haps2dataframe <- function(x) {
  d <- data.frame(x=as.vector(row(x)), y=as.vector(col(x)), z=as.vector(x))
  d[order(d$x), ]
}

emPlot <- function(res, sim) {
  real_haps <- sim$haplos[[1]]
  r1 <- confusionMatrix(res$haplos[, 1], real_haps[, 1])$overall['Accuracy']
  r2 <- confusionMatrix(res$haplos[, 1], real_haps[, 2])$overall['Accuracy']
  hap_order <- 1:2
  if (r2 > r1)
    hap_order <- 2:1
  #ani.record(reset = TRUE)
  var_sites <- !apply(real_haps == 0, 1, all) # all zero sites
  #orig <- haps2dataframe(data.frame(rep(FALSE, nrow(real_haps)), rep(FALSE, nrow(real_haps))))
  #levelplot(z ~ x + y, orig, col.regions=c("red","white"))
  iters <- as.integer(as.character(res$all_lls$iter))
  saveGIF({
    mapply(function(tt, i) {
         #print(table(tt[var_sites, ] == real_haps[var_sites, hap_order]))
         #p <- levelplot(z ~ x + y, haps2dataframe((tt[var_sites, ] == real_haps[var_sites, hap_order])), col.regions=c("red","white"))
         d <- haps2dataframe(1L + (tt[var_sites, ] == real_haps[var_sites, hap_order]))
         d$z <- factor(d$z, levels=1:2)
         #p <- ggplot(d) + geom_tile(aes(x=x,y=y, fill=z))
         p1 <- levelplot(z ~ x + y, d, col.regions=c("red","grey"), colorkey=FALSE, xlab="locus", ylab="haplotype\n(incorrect alleles are red)",
                         scales=list(y=list(at=1:2, labels=c("h1", "h2"))))
         p2 <- xyplot(ll ~ iter | ind, res$all_lls[iters <= i , ], group=hap, type='l', xlim=c(0, max(iters)),
                      ylab="log likelihood", xlab="iteration", ylim=c(1.1*min(res$all_lls$ll), 0.9*max(res$all_lls$ll)))
         grid.arrange(p1, p2, ncol=2)
         #ani.record()
    }, res$all_thetas, 1:max(iters))
    #oopts = ani.options(interval = 0.5)
  }, ani.width = 1000, ani.height=800)
}

EM plot of 80% heterozygous and 10% homozygous error.

set.seed(1)

sf <- sibFamily(40, 1000, ehet=0.8, ehom=0.1)
pa <- ProgenyArray(sf$progeny, sf$genos, mothers=sf$metadata$mother)
pa <- calcFreqs(pa)
pa <- inferParents(pa, 0.8, 0.1)

res <- emHap(pa, 1L, "parent_1", 1:nrow(progenyGenotypes(pa)), 0.8, 0.1)

emPlot(res, sf)

EM plot with low error.

sf <- sibFamily(40, 1000, ehet=0.1, ehom=0.1)
pa <- ProgenyArray(sf$progeny, sf$genos, mothers=sf$metadata$mother)
pa <- calcFreqs(pa)
pa <- inferParents(pa, 0.1, 0.1)

res <- emHap(pa, 1L, "parent_1", 1:nrow(progenyGenotypes(pa)), 0.1, 0.1)

emPlot(res, sf)

EM with homozygous mother and low error.

sf <- sibFamily(40, 1000, ehet=0.1, ehom=0.1, mom_inbred=TRUE)
pa <- ProgenyArray(sf$progeny, sf$genos, mothers=sf$metadata$mother)
pa <- calcFreqs(pa)
pa <- inferParents(pa, 0.1, 0.1)

res <- emHap(pa, 1L, "parent_1", 1:nrow(progenyGenotypes(pa)), 0.1, 0.1)

emPlot(res, sf)

table(sf$haplos[[1]][,1], res$haplos[,1])

EM with homozygous mother and high error.

sf <- sibFamily(30, 1000, ehet=0.8, ehom=0.1, mom_inbred=TRUE)
pa <- ProgenyArray(sf$progeny, sf$genos, mothers=sf$metadata$mother)
#pa <- calcFreqs(pa)
freqs(pa) <- sf$freqs
pa <- inferParents(pa, 0.8, 0.1)

res <- emHap(pa, 1L, "parent_1", 1:nrow(progenyGenotypes(pa)), 0.8, 0.2)

emPlot(res, sf)

table(sf$haplos[[1]][,1], res$haplos[,1])

Note: in re-running and experimenting with this, it seems that the model has an issue in choosing whether to say that $\mathbf{\pi} = {0.5, 0.5}$ and haplotypes are identical, or weighting one $\pi_k$ more heavily, and letting the other haplotype pick up error.

Label switching and identifiability

There are two potential identifiability issues: (1) label switching or degeneracy, and (2) a more serious issue when mixing two distributions gives us distribution in the same family (see Shalizi's discussion on this on page 3). This later identifiability is called generic identifiability.

We have (2) in the case of a homozygous mother, becauese both haplotype densities are the same (with equal mixing probability). As an example, suppose the mother's haplotypes are homozygous in a region. This means that $\theta_1 = \theta_2$, $\mathbf{\pi} = {0.5, 0.5}$ or $\theta_1$ is the haplotype and $\mathbf{\pi} = {1, 0}$.

vsbuffalo/ProgenyArray documentation built on May 3, 2019, 6:41 p.m.