前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >使用vcftools根据vcf文件计算种群核苷酸多样性

使用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

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

下载数据
代码语言:javascript
复制
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工具将所有的样本名重定向到一个文件里
代码语言:javascript
复制
bcftools query -l Massoko_Dryad_VCF_final.vcf.gz > inds_to_keep.txt

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

提取指定的样本并删除indel
代码语言:javascript
复制
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

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

代码语言:javascript
复制
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

接下来计算两个不同群体的核苷酸多样性
  • 获得两个不同群体所有的样本名,存入文件中
代码语言:javascript
复制
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
计算群体核苷酸多样性
代码语言:javascript
复制
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做折线图和箱线图
  • 箱线图
代码语言:javascript
复制
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

  • 折线图
代码语言:javascript
复制
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 条评论
热度
最新
推荐阅读
目录
  • 下载数据
  • 提取指定的样本并删除indel
  • 接下来计算两个不同群体的核苷酸多样性
  • 计算群体核苷酸多样性
  • 将结果文件导出,使用ggplot2做折线图和箱线图
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档