这个例子来源于一篇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格式的文件接下来改怎么处理
这里用到的是线粒体基因组的序列，线粒体相当于单倍体，如果是核基因组两倍体会有不一样的地方吗？
慢慢学习吧，希望可以找到答案！
推文的示例数据和代码大家可以直接找到开头提到的论文附件，或者直接给推文打赏
1元，入股打赏了没有收到我的回复，可以加我的微信mingyan24催我
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
