使用 vcftools 从 vcf 文件计算群体核苷酸多样性
最编程
2024-03-02 16:11:27
...
找到了一份种群基因组学数据分析的教程，原文用的数据是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