BCFTOOLS VIEW 命令查看VCF文件

发表于 | 分类于 | | 热度 °C

BCFTOOLS VIEW 可以用于查看VCF或BCF文件，VCF/BCF互相转换，提取子集，或进行过滤，用法如下：

1
2
3
4
5
6
bcftools view [options] <in.vcf.gz> [region1 [...]]
# 只针对压缩的VCF文件，即vcf.gz，若不是压缩格式，先转换成压缩格式，比如现在有一个文件，xxx.vcf，使用bgzip进行压缩，
> bgzip xxx.vcf
# 得到xxx.vcf.gz，生成索引
> bcftools index -t xxx.vcf.gz
# 得到xxx.vcf.gz.tbi，即可开始使用BCFTOOLS VIEW进行操作

参数如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
About:   VCF/BCF conversion, view, subset and filter VCF/BCF files.
Usage:   bcftools view [options] <in.vcf.gz> [region1 [...]]

Output options:
    -G,   --drop-genotypes              drop individual genotype information (after subsetting if -s option set)
    -h/H, --header-only/--no-header     print the header only/suppress the header in VCF output
    -l,   --compression-level [0-9]     compression level: 0 uncompressed, 1 best speed, 9 best compression [-1]
          --no-version                  do not append version and command line to the header
    -o,   --output-file <file>          output file name [stdout]
    -O,   --output-type <b|u|z|v>       b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF [v]
    -r, --regions <region>              restrict to comma-separated list of regions
    -R, --regions-file <file>           restrict to regions listed in a file
    -t, --targets [^]<region>           similar to -r but streams rather than index-jumps. Exclude regions with "^" prefix
    -T, --targets-file [^]<file>        similar to -R but streams rather than index-jumps. Exclude regions with "^" prefix
        --threads <int>                 number of extra (de)compression threads [0]

Subset options:
    -a, --trim-alt-alleles        trim alternate alleles not seen in the subset
    -I, --no-update               do not (re)calculate INFO fields for the subset (currently INFO/AC and INFO/AN)
    -s, --samples [^]<list>       comma separated list of samples to include (or exclude with "^" prefix)
    -S, --samples-file [^]<file>  file of samples to include (or exclude with "^" prefix)
        --force-samples           only warn about unknown subset samples

Filter options:
    -c/C, --min-ac/--max-ac <int>[:<type>]      minimum/maximum count for non-reference (nref), 1st alternate (alt1), least frequent
                                                   (minor), most frequent (major) or sum of all but most frequent (nonmajor) alleles [nref]
    -f,   --apply-filters <list>                require at least one of the listed FILTER strings (e.g. "PASS,.")
    -g,   --genotype [^]<hom|het|miss>          require one or more hom/het/missing genotype or, if prefixed with "^", exclude sites with hom/het/missing genotypes
    -i/e, --include/--exclude <expr>            select/exclude sites for which the expression is true (see man page for details)
    -k/n, --known/--novel                       select known/novel sites only (ID is not/is '.')
    -m/M, --min-alleles/--max-alleles <int>     minimum/maximum number of alleles listed in REF and ALT (e.g. -m2 -M2 for biallelic sites)
    -p/P, --phased/--exclude-phased             select/exclude sites where all samples are phased
    -q/Q, --min-af/--max-af <float>[:<type>]    minimum/maximum frequency for non-reference (nref), 1st alternate (alt1), least frequent
                                                   (minor), most frequent (major) or sum of all but most frequent (nonmajor) alleles [nref]
    -u/U, --uncalled/--exclude-uncalled         select/exclude sites without a called genotype
    -v/V, --types/--exclude-types <list>        select/exclude comma-separated list of variant types: snps,indels,mnps,ref,bnd,other [null]
    -x/X, --private/--exclude-private           select/exclude sites where the non-reference alleles are exclusive (private) to the subset samples

一、查找条目

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
> bcftools view -e 'GT="1|1"' xxx.vcf.gz | wc -l # 查找不符合 GT="1|1" 的条目，并输出数量
> bcftools view -i 'GT="1|1"' xxx.vcf.gz | wc -l # 查找符合 GT="1|1" 的条目，并输出数量
# -i/e，前者是包含，后者是排除
> bcftools view -i 'AF=0.500' xxx.vcf.gz | wc -l # 查找符合 AF=0.500 的条目，并输出数量
> bcftools view -e 'AF=0.500' xxx.vcf.gz | wc -l # 查找不符合 AF=0.500 的条目，并输出数量
# 上面是针对VCF文件第8列，即INFO列，也可针对其他列，如QUAL或REF等
> bcftools view -i 'QUAL=59.28' xxx.vcf.gz | wc -l # 查找符合 QUAL=59.28 的条目，并输出数量
> bcftools view -e 'QUAL=59.28' xxx.vcf.gz | wc -l # 查找不符合 QUAL=59.28 的条目，并输出数量
# 特别注意，等号之后的字符串要加""，如REF="G"
> bcftools view -i 'REF="G"' xxx.vcf.gz | wc -l # 查找符合 QUAL=59.28 的条目，并输出数量
> bcftools view -e 'REF="G"' xxx.vcf.gz | wc -l # 查找不符合 QUAL=59.28 的条目，并输出数量

# 若需展示条目，可传递给less/more，
> bcftools view -i -H 'GT="1|1"' xxx.vcf.gz | less

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample1
chr1    187485  .       G       A       2007.06 .       AC=2;AF=1;AN=2;CNN_1D=-0.124;DP=76;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=39.48;QD=27.12;SOR=0.992    GT:AD:DP:GQ:PGT:PID:PL:PS       1|1:0,74:74:99:1|1:187485_G_A:2021,221,0:187485
chr1    827209  .       G       C       255.97  .       AC=2;AF=1;AN=2;CNN_1D=1.611;DP=6;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=49.2;QD=27.24;SOR=0.693       GT:AD:DP:GQ:PGT:PID:PL:PS       1|1:0,6:6:18:1|1:827209_G_C:270,18,0:827209
chr1    827212  .       C       G       255.97  .       AC=2;AF=1;AN=2;CNN_1D=2.44;DP=6;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=49.2;QD=28.2;SOR=0.693 GT:AD:DP:GQ:PGT:PID:PL:PS       1|1:0,6:6:18:1|1:827209_G_C:270,18,0:827209
chr1    827221  .       T       C       255.97  .       AC=2;AF=1;AN=2;CNN_1D=1.913;DP=6;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=49.2;QD=25;SOR=0.693  GT:AD:DP:GQ:PGT:PID:PL:PS       1|1:0,6:6:18:1|1:827209_G_C:270,18,0:827209

二、提取INDEL/SNP

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
> bcftools view -v indels xxx.vcf.gz | wc -l
# 11993
> bcftools view -v indels xxx.vcf.gz | less
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample1
chr1    15903   .       G       GC      59.28   .       AC=2;AF=1;AN=2;CNN_1D=-0.8;DP=3;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=30.13;QD=29.64;SOR=2.303     GT:AD:DP:GQ:PL  1/1:0,2:2:6:71,6,0
chr1    16734   .       TG      T       31.6    .       AC=1;AF=0.5;AN=2;BaseQRankSum=-0.674;CNN_1D=-6.239;DP=10;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=32.17;MQRankSum=-0.14;QD=3.51;ReadPosRankSum=-0.921;SOR=0.132       GT:AD:DP:GQ:PL  0/1:7,2:9:39:39,0,197
chr1    19190   .       GC      G       32.6    .       AC=1;AF=0.5;AN=2;BaseQRankSum=-0.669;CNN_1D=-7.533;DP=13;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=31.81;MQRankSum=-2.12;QD=2.51;ReadPosRankSum=0.099;SOR=0.33 GT:AD:DP:GQ:PL  0/1:11,2:13:40:40,0,330
chr1    63735   .       CCTA    C       367.61  .       AC=1;AF=0.5;AN=2;BaseQRankSum=0.385;CNN_1D=0.827;DP=12;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=46.67;MQRankSum=0.385;QD=25.36;ReadPosRankSum=0.674;SOR=0.307 GT:AD:DP:GQ:PL  0/1:1,9:10:15:375,0,15
chr1    189392  .       ACC     A       2306.03 .       AC=2;AF=1;AN=2;BaseQRankSum=-1.809;CNN_1D=0.211;DP=57;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=35.61;MQRankSum=1.12;QD=28.73;ReadPosRankSum=1.48;SOR=0.504      GT:AD:DP:GQ:PL  1/1:1,53:54:99:2320,118,0

# 可以直接输出文件，
> bcftools view -v indels xxx.vcf.gz > indel.vcf

# 同理，提取SNP,
> bcftools view -v snps xxx.vcf.gz | wc -l
# 58314
> bcftools view -v snps xxx.vcf.gz | less
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample1
chr1    16495   .       G       C       36.65   .       AC=1;AF=0.5;AN=2;BaseQRankSum=-0.967;CNN_1D=-4.292;DP=3;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=22;MQRankSum=0;QD=12.22;ReadPosRankSum=0.967;SOR=1.179       GT:AD:DP:GQ:PL  0/1:1,2:3:18:44,0,18
chr1    17614   .       G       A       54.64   .       AC=1;AF=0.5;AN=2;BaseQRankSum=0;CNN_1D=-0.693;DP=5;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=47.88;MQRankSum=-1.645;QD=10.93;ReadPosRankSum=1.04;SOR=1.179     GT:AD:DP:GQ:PL  0/1:2,3:5:35:62,0,35
chr1    19322   .       C       T       44.64   .       AC=1;AF=0.5;AN=2;BaseQRankSum=0.703;CNN_1D=-3.724;DP=8;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=38.66;MQRankSum=-0.956;QD=5.58;ReadPosRankSum=0.414;SOR=0.818 GT:AD:DP:GQ:PGT:PID:PL:PS       0|1:5,3:8:52:0|1:19322_C_T:52,0,140:19322
chr1    19342   .       G       A       58.64   .       AC=1;AF=0.5;AN=2;BaseQRankSum=1.65;CNN_1D=-3.127;DP=5;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=41.16;MQRankSum=-0.524;QD=11.73;ReadPosRankSum=-1.036;SOR=1.609        GT:AD:DP:GQ:PGT:PID:PL:PS       0|1:3,2:5:66:1|0:19322_C_T:66,0,116:19322
chr1    29443   .       A       G       41.32   .       AC=2;AF=1;AN=2;CNN_1D=-3.532;DP=3;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=22;QD=20.66;SOR=2.303      GT:AD:DP:GQ:PL  1/1:0,2:2:6:53,6,0
# 直接输出文件
> bcftools view -v snps xxx.vcf.gz > snp.vcf

# 这个VCF总条目数为70294，INDEL=11993，SNP=58314，合计70307，是否有既是SNP又是INDEL的条目？
> bcftools view -v indels xxx.vcf.gz | bcftools view -v snps | wc -l
# 511
> bcftools view -v indels xxx.vcf.gz | bcftools view -v snps | less

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample1
chr1    78005376        .       T       C,TTC   268.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=0.499;DP=9;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=29.78;SOR=1.402 GT:AD:DP:GQ:PL  1/2:0,6,3:9:84:285,101,84,175,0,196
chr10   58717368        .       A       C,AC    329.06  .       AC=1,1;AF=0.5,0.5;AN=2;BaseQRankSum=2.17;CNN_1D=0.755;DP=13;ExcessHet=3.0103;FS=3.358;MLEAC=1,1;MLEAF=0.5,0.5;MQ=58.98;MQRankSum=0;QD=25.31;ReadPosRankSum=1.61;SOR=0.707       GT:AD:DP:GQ:PL  1/2:1,7,5:13:99:346,129,102,168,0,202
chr13   35044829        .       TAG     T,GAG   110.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=0.437;DP=4;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=27.51;SOR=0.693 GT:AD:DP:GQ:PL  1/2:0,2,2:4:49:127,49,78,84,0,78
chr13   40558927        .       G       T,GT    150.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=-0.053;DP=5;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=27.85;SOR=0.693        GT:AD:DP:GQ:PL  1/2:0,2,2:4:62:167,72,62,74,0,64
chr17   67072015        .       A       AC,C    395.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=0.428;DP=10;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=28.72;SOR=1.863        GT:AD:DP:GQ:PL  1/2:0,4,4:8:89:412,119,89,119,0,89
chr17   76292181        .       T       C,TGCTGAACTGCACCAGGTTGGACCAAACCAC       425.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=3.964;DP=19;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=26.76;SOR=1.27 GT:AD:DP:GQ:PL  1/2:0,7,4:11:99:442,150,230,184,0,200
chr2    97151827        .       C       CT,T    441.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=0.512;DP=12;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=29.84;QD=31.8;SOR=1.022      GT:AD:DP:GQ:PL  1/2:0,6,6:12:99:458,217,193,220,0,234
chr2    197008003       .       T       TTCTC,C 235.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=3.73;DP=6;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=31.03;SOR=0.693  GT:AD:DP:GQ:PL  1/2:0,4,2:6:31:252,49,31,99,0,81
chr22   24755923        .       CA      C,AA    209.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=1.317;DP=8;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=29.87;SOR=0.941 GT:AD:DP:GQ:PL  1/2:0,4,3:7:62:226,82,62,103,0,119
chr3    149515881       .       T       G,TTG   280.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=1.794;DP=9;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=31.12;SOR=0.892 GT:AD:DP:GQ:PL  1/2:0,4,5:9:96:297,208,194,96,0,132
chr6    125929252       .       T       TAAAA,A 193.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=1.255;DP=5;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=28.92;SOR=1.022 GT:AD:DP:GQ:PL  1/2:0,2,3:5:69:210,97,88,80,0,69
chr7    98315625        .       A       ATAAT,T 371.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=1.314;DP=10;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=33.02;SOR=0.941        GT:AD:DP:GQ:PL  1/2:0,5,2:7:28:388,58,28,148,0,120
chr8    19404624        .       T       A,TATA  295.06  .       AC=1,1;AF=0.5,0.5;AN=2;CNN_1D=1.322;DP=11;ExcessHet=3.0103;FS=0;MLEAC=1,1;MLEAF=0.5,0.5;MQ=60;QD=31.75;SOR=1.179        GT:AD:DP:GQ:PL  1/2:0,2,4:8:59:312,193,178,59,0,67
(END)
# 果如预料，在13个条目中，二倍体中的两个等位基因皆发生突变，而且突变后的ALT不同，因此，会同时被认定是INDEL或SNP
# 注意，上面提到INDEL+SNP=70307，而文件的总条目是70294，多出13条，代表有13条记录既是INDEL又是SNP，我们确实看到13条相应的条目，但为何使用“wc -l”命令查看时却有511条记录呢？这是因为VCF文件头部有注释行，正确的查看方式是使用"grep -v "^#""命令，略过注释行

> bcftools view -v indels xxx.vcf.gz | bcftools view -v snps | grep -v "^#" | wc -l
# 13

相应地，我们看一下真实的INDEL与SNP的条目数，
> bcftools view -v snps xxx.vcf.gz | grep -v "^#" | wc -l
# 57817，相比58314，少497行，即有497行注释
bcftools view -v indels xxx.vcf.gz | grep -v "^#" | wc -l
# 11496，相比11993，少少497行，即有497行注释
# 另外，同时是INDEL和SNP的行有13行，如果加上注释为511行，注释497行，还有一行是标题行
  • 本文作者：括囊无誉
  • 本文链接： WES/bcftools_view/
  • 版权声明： 本博客所有文章均为原创作品，转载请注明出处！
------ 本文结束 ------
坚持原创文章分享，您的支持将鼓励我继续创作！