vcf文件与vcftools(二)

vcftools是用来处理vcf和bcf文件的工具集,功能涵盖了过滤,数据格式转换,一些指标的统计计算,vcf文件之间变异信息的差异比较等等。

1.基本用法
vcftools [ --vcf FILE | --gzvcf FILE | --bcf FILE] \
[ --out OUTPUT PREFIX ] [ FILTERING OPTIONS ] [ OUTPUT OPTIONS ]
2.使用

官网提供了几个常用示例,下面用我的vcf文件测试一下,我的vcf文件共249个样本,27012个SNP标记

2.1 对来自chr1的每一个位点统计其基因频率
vcftools --gzvcf combined200.vcf.gz --freq --chr chr1 --out chr1_analysis

因为没有加任何过滤选项,所以没有过滤任何记录

$ less chr1_analysis.log | tail -n 4
After filtering, kept 249 out of 249 Individuals
Outputting Frequency Statistics...
After filtering, kept 3453 out of a possible 27012 Sites #chr1上原本就有3453个SNP位点
Run Time = 44.00 seconds

第三列表示在所有样本中,该位点出现了几种碱基,大多数情况下为2,也可能是3、4;第四列表示该位点出现的碱基总数,一个样本贡献2个。

$ head -n 5 chr1_analysis.frq 
CHROM   POS N_ALLELES   N_CHR   {ALLELE:FREQ}
chr1    1146710 2   494 G:0.983806  A:0.0161943
chr1    1146714 2   494 G:0.989879  A:0.0101215
chr1    1146763 2   494 G:0.995951  A:0.00404858
chr1    1146783 2   494 C:0.995951  T:0.00404858
2.2 过滤掉vcf文件中的InDel记录,只保留SNP记录
vcftools --gzvcf combined200.vcf.gz --remove-indels --recode --recode-INFO-all --out SNPs_only

--recode 
表示过滤之后会生成一个新文件,以.recode.vcf为后缀
--recode-INFO-all
因为在过滤之后,原先存在的INFO列的注释信息可能不对,
比如剔除了一些样本,那么AN就需要重新计算。
所以过滤之后默认情况下是不含有INFO列的。
该选项表示将原来vcf文件中所有INFO信息保留下来。
--recode-INFO <string> 
比如--recode-INFO AC表示只将AC保留
--keep-only-indels和--remove-indels作用相反

根据vcf文件的FILTER列来过滤
--remove-filtered-all 过滤掉FILTER列不是“PASS”的记录,这一步grep也能做。

vcftools --vcf combined3000_head3k.vcf --remove-filtered-all --recode --out SNPs_filter

--max-missing
--max-missing的取值是0-1,为1时表示某个位点上所有的样本必须都有基因型,一个样本的基因型都不能缺。所以这个选项可以理解为:能分型的样本占总样本的比例至少为多少。

vcftools --gzvcf combined200.vcf.gz --max-missing 0.99 --recode --out output_noMissing

有过滤操作都要加上--recode

2.3 比较两个vcf文件的变异位点
vcftools --gzvcf combined200.vcf.gz --diff output_noMissing.recode.vcf --diff-site --out in1_v_in2

下面的表头的1、2分别表示文件1、文件2,第4列1或2表示仅文件1或文件2有,B表示都有

$ head in1_v_in2.diff.sites_in_files 
CHROM   POS1    POS2    IN_FILE REF1    REF2    ALT1    ALT2
chr1    1146710 1146710 B   G   G   A   A
chr1    1146714 1146714 B   G   G   A   A
chr1    1146763 1146763 B   G   G   A   A
chr1    1146783 1146783 B   C   C   T   T
chr1    1146785 1146785 B   G   G   C   C
chr1    1146852 1146852 B   C   C   T   T
chr1    1146853 1146853 B   G   G   A   A
chr1    1147062 .   1   G   .   A   .
chr1    1147243 .   1   A   .   C   .
2.4 结合管道操作符

用--stdout代替--out即可结合管道进行其他处理

2.5 计算Hardy-Weinberg p-value

Hardy-Weinberg遗传平衡检验用到的统计方法是卡方测验,详细过程参考:http://www.dxy.cn/bbs/thread/223387#223387以及https://wenku.baidu.com/view/9820fc2258fb770bf78a5571.html

vcftools --gzvcf combined200.vcf.gz --hardy --out each_site_hardy   
# Only using biallelic SNPs

第6列即为所求

$ lsx each_site_hardy.hwe | head -n 2
CHR POS OBS(HOM1/HET/HOM2)  E(HOM1/HET/HOM2)    ChiSq_HWE   P_HWE   P_HET_DEFICIT   P_HET_EXCESS
chr1    1146710 239/8/0 239.06/7.87/0.06    6.692747e-02    1.000000e+00    1.000000e+00    9.440689e-01
3.其它参数说明
3.1 基本选项

--vcf | --gzvcf | --bcf 输入格式,三选一即可
--out <输出文件前缀>
--stdout 和 -c 输出到标准输出,能通过管道与其它命令结合

3.2 过滤SNP的选项

POSITION FILTERING
可以提供一些位置信息,或是能表示位置的文件,具体参数见官方文档
ALLELE FILTERING
--maf 和 --max-maf用来限定最小等位基因频率(MAF)的范围
--mac 和 --max-mac和上面类似,用来限定最小等位基因的数量
--min-alleles 和 --max-alleles用来限定等位基因的数量,比如这条命令只会保留vcf文件ALT列只有一种碱基的情况。

vcftools --gzvcf combined200.vcf.gz --min-alleles 2 --max-alleles 2

GENOTYPE VALUE FILTERING
--min-meanDP 和 --max-meanDP用来限定所有样本DP的平均值,DP表示某一个样本某一位点所有allele的总深度
--hwe 2.5 计算Hardy-Weinberg p-value讲到如何求p值,这个参数就是根据p值来过滤的,小于阈值则被过滤掉
--max-missing 前面已经举例;--max-missing-count 某个位点缺失样本个数多于某个阈值则过滤掉
--phased 某个位点如果含有未定相的基因型则过滤
--minQ 根据vcf文件的QUAL列来过滤,比如

vcftools --gzvcf combined200.vcf.gz --minQ 100 --recode --out minQ_100
3.3 过滤样本的选项

--indv 和 --remove-indv保留或去除哪一个样本,后接样本ID
--keep 和 --remove保留或去除哪些样本,后接文件名,文件中一行一个样本ID

3.4 输出选项

等位基因统计
--freq 和 --counts用于统计等位基因的频率和数量
深度统计
每个样本所有位点的平均深度、单个位点所有样本的深度加和(或平均)、每个样本每个位点基因型的深度。感觉没什么用
LD统计
Ti/Tv统计
NUCLEOTIDE DIVERGENCE STATISTICS
FST STATISTICS
OTHER STATISTICS
格式转换

3.5 比较选项

reference

The variant call format and VCFtools: https://academic.oup.com/bioinformatics/article/27/15/2156/402296
VCFtools--A set of tools written in Perl and C++ for working with VCF files:
https://vcftools.github.io/man_latest.html#EXAMPLES

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,175评论 5 466
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,674评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,151评论 0 328
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,597评论 1 269
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,505评论 5 359
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 47,969评论 1 275
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,455评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,118评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,227评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,213评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,214评论 1 328
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,928评论 3 316
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,512评论 3 302
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,616评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,848评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,228评论 2 344
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,772评论 2 339