计算群体pi

1、提取3k材料的名称

'''

bcftools query -l 3k.snp.vcf > 3k.list

'''

#2、挑选要计算pi的subset

#3、使用vcftools 计算pi #以30kb的窗口进行计算

vcftools --vcf 3k.snp.vcf  --keep group1.list --window-pi 30000  --out group1.int

导入R

#读入数据

p1 <- read.table("group1.int.windowed.pi", header = T)

#按照群体命名

p1$pop <- "group1"

###合并三个文件

###df <- rbind(p1,p2,p3)

#画图

library(ggplot2)

options(scipen = 200)

ggplot(p1, aes(x=pop, y = PI, fill = pop)) +

  geom_boxplot() +theme_classic()

如果是展示多个，可以使用grid或者facet_wrap进行展示