提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  其他

当前栏目

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

网络 语言
2023-03-20 15:38:10 时间

这个例子来源于一篇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)

image.png

能运行完代码，但是还有很多疑问，

  • 首先是单倍型的图怎们看
  • 怎么获取画图数据然后用ggplot2来画图

还有的论文中会得到一个表格

image.png

怎么才能得到这个单倍型的序列。

先在的群体大部分都是snp数据，对应的vcf文件，如果拿到vcf格式的文件接下来改怎么处理

这里用到的是线粒体基因组的序列，线粒体相当于单倍体，如果是核基因组两倍体会有不一样的地方吗？

慢慢学习吧，希望可以找到答案！

相关文章