今天推文的主要內(nèi)容參考 鏈接 https://wurmlab.github.io/genomicscourse/2016-SIB/practicals/population_genetics/popgen，示例vcf格式數(shù)據(jù)下載自https://github.com/wurmlab/genomicscourse/tree/master/2016-SIB/data/popgen/vcf， 大家可以自己到鏈接下載示例數(shù)據(jù)，也可以給這篇推文點(diǎn)贊留言獲取數(shù)據(jù)

首先是使用bcftools軟件操作vcf文件

• 將vcf文件按照染色體拆分

bcftools view snp.vcf.gz scaffold_1 > popgenome-vcf/scaffold_1
bcftools view snp.vcf.gz scaffold_2 > popgenome-vcf/scaffold_2

如果當(dāng)前目錄下只有vcf格式文件，會(huì)遇到報(bào)錯(cuò)Failed to open .vcf.gz: could not load index，可以參考 https://www.cnblogs.com/chenwenyan/p/11945445.html

 tabix -p vcf snp.vcf.gz

如果當(dāng)前目錄下沒(méi)有popgenome-vcf這個(gè)目錄，還需要新建目錄

mkdir popgenome-vcf

今天參考的文章里寫(xiě)道 In theory, the r PopGenome can read VCF files directly, using the readVCF function. However, because our samples are haploid, we need to use a different function, r readData, which requires a folder with a separate VCF for each scaffold. 這個(gè)是為什么呢？

接下來(lái)是在R語(yǔ)言里的操作

讀入數(shù)據(jù)

#install.packages("PopGenome")
library(PopGenome)
getwd()
setwd("VCF/")
snp<-readData("popgenome-vcf",format = "VCF")

統(tǒng)計(jì)一些基本信息

get.sum.data(snp)

image.png

這里可以直接統(tǒng)計(jì) 轉(zhuǎn)換和顛換的比例

獲取樣本名稱(chēng)并分組

pops<-get.individuals(snp)[[1]]
pop1<-pops[grep("B\\.bam",pops)]
pop2<-pops[grep("b\\.bam",pops)]
pop1
pop2

給數(shù)據(jù)劃分類(lèi)群

snp<-set.populations(snp,list(pop1,pop2))
snp@populations

計(jì)算FST

snp<-F_ST.stats(snp)
get.F_ST(snp)

image.png

這里的指標(biāo)都是什么意思呢？

計(jì)算核苷酸多樣性

get.diversity(snp)[[1]]

image.png

這里的指標(biāo)也看不懂是什么意思呀

設(shè)置劃窗計(jì)算指標(biāo)

win_snp<-sliding.window.transform(snp,
                                  width = 10000,
                                  jump = 2000,type = 2)
win_snp<-F_ST.stats(win_snp)



win_snp@nucleotide.F_ST
win_snp@nuc.diversity.within

接下來(lái)用折線圖來(lái)展示結(jié)果

FST

library(ggplot2)
win_fst <- data.frame(x=1:dim(win_snp@nucleotide.F_ST)[1],
                      y=win_snp@nucleotide.F_ST[,1])
head(win_fst)
p1<-ggplot(win_fst,aes(x=x,y=y))+
  geom_point()+
  geom_line()+
  theme_bw()+
  theme(panel.grid = element_blank())+
  scale_x_continuous(breaks = win_fst$x,
                     labels = win_fst$x)+
  labs(x=NULL,y=NULL,title = "FST")
ggsave("FST.pdf",p1,width = 15,height = 4)

image.png

核苷酸多樣性

bb_div  <- win_snp@nuc.diversity.within[,1] # diversity among B (bb = "big B")
lb_div  <- win_snp@nuc.diversity.within[,2] # diversity among B (lb = "little b")
bb_div
df1<-data.frame(x=1:length(bb_div),y=bb_div)
df2<-data.frame(x=1:length(lb_div),y=lb_div)
p2<-ggplot()+
  geom_line(data=df1,aes(x=x,y=y),color="red")+
  geom_point(data=df1,aes(x=x,y=y),size=2,color="red")+
  geom_line(data=df2,aes(x=x,y=y),color="blue")+
  geom_point(data=df2,aes(x=x,y=y),size=2,color="blue")+
  theme_bw()+labs(x=NULL,y=NULL)
ggsave("diversity.pdf",p2,width = 15,height = 4)

image.png

歡迎大家關(guān)注我的公眾號(hào)
小明的數(shù)據(jù)分析筆記本