计算杂合度函数calHe(配合polysat包使用)
calHe <- function(fre, filename = "hets.txt"){ pops = row.names(fre) loci = unique(as.matrix(as.data.frame(strsplit(names(fre), split = ".", fixed = TRUE), stringsAsFactors = FALSE))[1, ]) loci <- loci[loci != "Genomes"] hets <- array(0, dim = c(length(pops),length(loci)), dimnames = list(pops, loci)) for(P in pops){ freqs <- fre[P, ] genomes <- freqs$Genomes names(genomes) <- row.names(freqs) thesegenomes <- genomes for (L in loci) { thesefreqs <- freqs[, grep(paste("^", L, "\\.", sep = ""), names(freqs)), drop = FALSE] thesefreqs <- thesefreqs[, names(thesefreqs) != paste(L, "Genomes", sep = "."), drop = FALSE] hsByPop <- apply(as.matrix(thesefreqs), 1, function(x) 1 - sum(x^2)) avgfreq <- unlist(lapply(thesefreqs, weighted.mean, w = thesegenomes)) hets[P,L] <- 1 - sum(avgfreq^2) } } hets<-data.frame(hets) hets$mean <- apply(hets,1,mean) write.table(hets, filename, sep = "\t",quote = FALSE) return(hets) }