這個例子來源于一篇plos的論文
論文題目是
A workflow with R: Phylogenetic analyses and visualizations using mitochondrial cytochrome b gene sequences
論文提供了完整的R語言代碼和示例數(shù)據(jù)
今天的推文試著重復(fù)一下里面單倍型網(wǎng)絡(luò)的代碼
單倍型到底是個啥還是沒有搞明白
首先是示例數(shù)據(jù)集
- 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
函數(shù)用到的是pegas::haplotype但是用到的參數(shù)還不知道是啥意思
計算單倍型網(wǎng)絡(luò)
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)
這個是針對個體的
還有一個針對群體的
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)
能運行完代碼吆倦，但是還有很多疑問，
- 首先是單倍型的圖怎們看
- 怎么獲取畫圖數(shù)據(jù)然后用ggplot2來畫圖
還有的論文中會得到一個表格
怎么才能得到這個單倍型的序列。
先在的群體大部分都是snp數(shù)據(jù)取逾，對應(yīng)的vcf文件径密，如果拿到vcf格式的文件接下來改怎么處理
這里用到的是線粒體基因組的序列，線粒體相當(dāng)于單倍體，如果是核基因組兩倍體會有不一樣的地方嗎唤蔗？
慢慢學(xué)習(xí)吧义锥，希望可以找到答案柳沙！
推文的示例數(shù)據(jù)和代碼大家可以直接找到開頭提到的論文附件，或者直接給推文打賞
1元拌倍，入股打賞了沒有收到我的回復(fù)赂鲤，可以加我的微信mingyan24催我
歡迎大家關(guān)注我的公眾號
小明的數(shù)據(jù)分析筆記本
小明的數(shù)據(jù)分析筆記本 公眾號 主要分享：1噪径、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子；2蛤袒、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)熄云、基因組學(xué)、群體遺傳學(xué)文獻閱讀筆記妙真；3缴允、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！
