使用vcftools根据vcf文件计算种群核苷酸多样性

用户7010445

发布于 2020-03-03 14:54:18

8.5K2

发布于 2020-03-03 14:54:18

找到了一份种群基因组学数据分析的教程，原文用的数据是2015年发表在science上的一篇论文Genomic islands of speciation separate cichlid ecomorphs in an East African crater lake。这份教程利用这篇文章的数据分析了部分内容。

教程地址 https://github.com/singhal/popgen_tutorial/blob/master/pop_gen_tutorial.rst
vcf文件下载链接 http://datadryad.org/bitstream/handle/10255/dryad.101389/Massoko_Dryad_VCF_final.vcf.gz

争取把教程中的内容都重复完，然后多看几遍原论文（高水平论文看起来还真有些吃力呢！）

下载数据

wget http://datadryad.org/bitstream/handle/10255/dryad.101389/Massoko_Dryad_VCF_final.vcf.gz

为了减小计算压力，教程中的处理方式是只保留36个样本（正常数据中好像是有146个样本，解压出来的vcf文件有11G），并且删除了inde，只保留snp位点。但是原文中保留的36个个体的文本文件inds_to_keep.txt我现在找不到，需要自己重新构造一份需要保留的个体的样本名。处理方式是：

首先使用bcftools工具将所有的样本名重定向到一个文件里

bcftools query -l Massoko_Dryad_VCF_final.vcf.gz > inds_to_keep.txt

我选择的是每个群体保留六个样本（样本前缀名一直我就认为他们是来自同一个群体），最后我保留了38个个体这一步大家可以自行进行处理或者给我留言获得inds_to_keep.txt文件。

提取指定的样本并删除indel

vcftools --gzvcf Massoko_Dryad_VCF_final.vcf.gz --keep inds_to_keep.txt --stdout --recode --recode-INFO-all --remove-indels | bgzip  > Massoko_Dryad_VCF_final_subset_noIndels.vcf.gz

为了减小计算压力，进一步对文件进行处理（这一步使用到的两个参数自己还不太明白是什么意思，这一步完全照搬原教程）

vcftools --gzvcf Massoko_Dryad_VCF_final_subset_noIndels.vcf.gz --maf 0.05 --max-maf 0.95 --stdout --recode --recode-INFO-all | bgzip > Massoko_Dryad_VCF_final_subset_noIndels_maf05.vcf.gz
vcftools --gzvcf Massoko_Dryad_VCF_final_subset_noIndels_maf05.vcf.gz --thin 1000 --stdout --recode --recode-INFO-all | bgzip > Massoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf.gz

这里不明白的参数 --maf --max-maf通常会设置最小等位基因频率来过滤vcf文件，但这里设置最大等位基因频率是什么意思？ --thin 1000

接下来计算两个不同群体的核苷酸多样性

获得两个不同群体所有的样本名，存入文件中

bcftools query -l Massoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf | grep "littoral" > littoral.txt
bcftools query -l Massoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf | grep "benthic" > benthic.txt

计算群体核苷酸多样性

vcftools --vcf Massoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf --keep littoral.txt --window-pi 100000 --out littoral_pi
vcftools --vcf Massoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf --keep benthic.txt --window-pi 100000 --out benthic_pi

--window-pi 指定窗口的长度 --out 指定输出文件的前缀名

将结果文件导出，使用ggplot2做折线图和箱线图

箱线图

bb<-read.table("../../vcf_handling/Fish_Populations/benthic_pi.windowed.pi",header=T)
ll<-read.table("../../vcf_handling/Fish_Populations/littoral_pi.windowed.pi",header=T)
dim(bb)
head(bb)
bb$indiv<-"benthic"
dim(ll)
head(ll)
ll$indiv<-"littoral"
df<-rbind(bb,ll)
colnames(df)
dim(df)
library(ggplot2)
options(scipen=200)
ggplot(df,aes(x=indiv,y=PI,fill=indiv))+
  geom_boxplot()+theme_bw()

image.png

折线图

ggplot(df,aes(x=BIN_START,y=PI,group=indiv,color=indiv))+
  geom_line()+theme_bw()+
  theme(legend.position = "top",
        legend.title = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title = element_blank())

image.png

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2019-08-30，如有侵权请联系 cloudcommunity@tencent.com 删除

数据分析

本文分享自小明的数据分析笔记本微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

数据分析

登录后参与评论

0 条评论

热度