I think I found a mistake in gl.report.heterozygosity with the position of one of the na.rm arguments. Currnetly it doesnt remove missing data when calculating mean heterozygosity per locus with colMeans(), only when averaging across loci with mean(). Therefore it will only remove loci for which there are read errors for all indivudals (highly unlikley) but not NAs 'peppered' through dataset.
See comments below.
function (gl) {
x <- gl
sgl <- seppop(x)
hs <- data.frame(lapply(sgl, function(x) mean(colMeans(as.matrix(x,
na.rm = TRUE) == 1), na.rm = TRUE))) ## should be arg for colMeans() i.e: na.rm = TRUE) == 1, na.rm = TRUE))))
sums <- data.frame(lapply(sgl, function(x) mean(colSums(as.matrix(x,
na.rm = TRUE) == 1), na.rm = TRUE))) ## and here
n <- t(sums/hs)
n <- cbind(row.names(n), n)
d1 <- data.frame(names = names(hs), hs = as.numeric(hs[1,
]))
d2 <- data.frame(n)
names(d2) <- c("names", "freq")
dd <- join(d1, d2, by = "names")
dd <- dd[order(dd$hs), ]
par(mai = c(2.5, 1, 0.5, 0.2))
barplot(dd$hs, names.arg = paste(dd$names, dd$freq, sep = " | "),
las = 2, cex.names = 1, space = 0, border = F, col = rainbow(nrow(dd)),
main = "Observed Heterozygosity")
return(dd)
}
You will also need to fix up how n is calculated or displayed in the plot....