R语言做单倍型网络（haplotype network）的一个小例子

用户7010445

发布于 2021-12-01 18:14:34

2.4K0

发布于 2021-12-01 18:14:34

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

这个例子来源于一篇plos的论文

论文题目是

A workflow with R: Phylogenetic analyses and visualizations using mitochondrial cytochrome b gene sequences

image.png

论文提供了完整的R语言代码和示例数据

今天的推文试着重复一下里面单倍型网络的代码

单倍型到底是个啥还是没有搞明白

首先是示例数据集

120个熊蜂 Bombus terrestris dalmatinus
mitochondrial cyt b sequences (373 bp)
8个群体

读取fasta格式的DNA序列

library(ape)
nbin<-read.FASTA("pone.0243927.s002.fas")
class(nbin)

计算单倍型

library(pegas)
h<-pegas::haplotype(nbin,strict=FALSE,trailingGapsAsN=TRUE)
h
hname<-paste("H",1:nrow(h),sep="")
hname
rownames(h)<-hname
h

函数用到的是pegas::haplotype但是用到的参数还不知道是啥意思

计算单倍型网络

net<-pegas::haploNet(h,d=NULL,getProb = TRUE)
net
ind.hap<-with(
  utils::stack(setNames(attr(h, "index"), rownames(h))),
  table(hap=ind, individuals=names(nbin))
)
ind.hap
plot(net, size=attr(net, "freq"), scale.ratio = 2, cex = 0.6, labels=TRUE, pie = ind.hap, show.mutation=1, font=2, fast=TRUE)
legend(x= 57,y=15, colnames(ind.hap), fill=rainbow(ncol(ind.hap)), cex=0.52, ncol=6, x.intersp=0.2, text.width=11)

image.png

这个是针对个体的

还有一个针对群体的

h<-pegas::haplotype(nbin, 
                    strict = FALSE, 
                    trailingGapsAsN = TRUE)
hname<-paste("H", 1:nrow(h), sep = "")
rownames(h)<-paste(hname)
net<-haploNet(h, d = NULL, 
              getProb = TRUE) 
net
labels(nbin)
names(nbin)
ind.hap<-with(
  utils::stack(setNames(attr(h, "index"), rownames(h))),
  table(hap=ind, individuals=labels(nbin)[values])
)

ind.hap

bg<-c(rep("dodgerblue4", 15), 
      rep("olivedrab4",15), 
      rep("royalblue2", 15), 
      rep("red",15), 
      rep("olivedrab3",15), 
      rep("skyblue1", 15), 
      rep("olivedrab1", 15),  
      rep("darkseagreen1", 15))


hapcol<-c("Aksu", 
          "Demre", 
          "Kumluca", 
          "Firm", 
          "Bayatbadem", 
          "Geyikbayir", 
          "Phaselis", 
          "Termessos")


ubg<-c(rep("dodgerblue4",1), 
       rep("royalblue2",1), 
       rep("skyblue1",1), 
       rep("red",1), 
       rep("olivedrab4",1), 
       rep("olivedrab3",1),
       rep("olivedrab1",1), 
       rep("darkseagreen1",1))


plot(net, size=attr(net, "freq"), 
     bg = bg, 
     scale.ratio = 2, 
     cex = 0.7, 
     labels=TRUE, 
     pie = ind.hap, 
     show.mutation=1, 
     font=2, 
     fast=TRUE)
legend(x=-45,y=60, 
       hapcol, 
       fill=ubg,
       cex=0.8, 
       ncol=1, 
       bty="n",
       x.intersp = 0.2)