首页 > 其他 > 详细

计算杂合度函数calHe(配合polysat包使用)

时间:2020-03-11 23:20:48      阅读:141      评论:0      收藏:0      [点我收藏+]
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)
}

 

计算杂合度函数calHe(配合polysat包使用)

原文:https://www.cnblogs.com/liulele/p/12466350.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!