计算杂合度函数calHe（配合polysat包使用）

Posted 2021-03-06 liulele

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 = "	",quote = FALSE)
  return(hets)
}