這個例子來源于一篇plos的論文

論文題目是

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

image.png

論文提供了完整的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)

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

能運行完代碼，但是還有很多疑問，

• 首先是單倍型的圖怎們看
• 怎么獲取畫圖數(shù)據(jù)然后用ggplot2來畫圖

還有的論文中會得到一個表格

image.png

怎么才能得到這個單倍型的序列。

先在的群體大部分都是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é)文獻(xiàn)閱讀筆記；3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！