biostar handbook(十一)|基因组变异的表示形式

2018-02-06 本文已影响268人 xuzhougeng

VCF文件格式

在biostar handbook(十)|如何进行变异检测部分我们最后以VCF格式存放找到的变异。尽管大部分情况下，我们都不需要直接和VCF文件打交道，通常就是将其作为输入提供给后续的分析。但是，你对VCF的格式越熟悉，你就能使用bcftools,vcftools或其他工具提取你任意需要的数据。

VCF(Variant Call Format)可以用来存放找到的变异信息，包括三个部分，元信息，标题行和数据行。举个例子

案例

元信息(meta-information)行以"##"起始，首先是VCF的版本信息，后面的INFO定义和解释INFO列出现的ID的含义，FILTER解释说明做了过滤的类型，而FORMAT则解释FORMAT列出现的ID的含义和数据结构，SAMPLE则是告知有哪些样本，都比较的浅显易懂。

标题列固定8列，CHROM(染色体ID)， POS(变异所在位置)， ID(已有注释), REF(参考碱基)， ALT(变异碱基)， QUAL(质量)，FILTER(过滤方式)，INFO(总体信息列)，第9列为FORMAT，定义了后面每个样本的数据存放结构。

对于数据行，我比较在意的是如何在VCF存放和解读变异信息，所以也重点介绍这部分。

对于SNP和较小的indels

对于比较小的SNP，或者是插入缺失的碱基在20bp以内的indel，表示比较容易，读起来也不太费劲。

比如说在参考基因组和样本的基因组上某个位置上有如下情况

案例	序列	说明
Ref	a t C g a	参考序列
1	a t G g a	C突变成G
2	a t - g a	C缺失
3	a t CAg a	C后插入A

如果只有一个样本，表示方法位:

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO   # 解读
20      3   .   C   G   .       PASS    DP=100 # ref第三位是C，而ALT第三位是G
20      2   .   TC  T   .       PASS    DP=100 # ref第二位开始时TC，而ALT第二位开始只有T，说明ALT的第三位C缺失
20      2   .   TC  TCA .       PASS    DP=100 # ref第二位是TC，而ALT第二位开始时TCA，说明多了一个A

如果同时表示三个样本

#CHROM  POS ID  REF ALT      QUAL    FILTER  INFO   # 解读
20      2   .   TC  TG,T,TCA .       PASS    DP=100 # 表示在该位点上有三个突变信息

那么已知ref为atCga，根据VCF信息进行重组

#CHROM POS ID REF ALT QUAL FILTER INFO # 重组结果
20     3   .   C   T   .   PASS DP=100 #  atTga
20     3   .   C CTAG  .   PASS DP=100 #  atCTAGga
20     2   .   TCG T   .   PASS DP=100 #  aTa

结构变异(structure variants)

这里的定义SV，指的是插入缺失大于20bp，小于12kb的情况，先随意感受下VCF是如何处理这种情况。

结构变异

即为了表示SV,需要专门定义INFO和FORMAT。根据定义就能对这6个变异进行解读

一个准确的缺失, 发生在2827694-2827708
一个不太准确的缺失(DEL)，长度大概在205bp(SVLEN=-205)
一个不太准确的ALU缺失(DEL:ME:ALU),长度大概为209bp(SVLEN=-209)
一个不太准确的LA插入(INS:ME:L1)
也差不多

大规模重排

暂时不在讨论范围之内。

BCFtools

BCFtools是一套处理VCF和BCF格式的工具。它有提供许多子命令实现不同功能，我个人用的比较多有以下几个：

mpileup + call: 根据参考基因组寻找变异位点
view: 选取，过滤以及VCF/BCF之间的格式转换。这个命令完成了convert,filter,
query: VCF/BCF格式输出为更适合人类阅读的格式
merge: 将多个VCF/BCF文件整合成一个
isec: 求不同VCF/BCF文件的交集，合集和补集

通用参数

在单独介绍每个命令之前，需要了解一下所有子命令都可以用的参数：

文件输出： -o, --ouput FILE，默认输出到标准输出，通过该选项指定文件。 -O, --output-type b|u|z|v：输出为压缩的BCF(b)，未压缩的BCF(u), 压缩的VCF(z), 未压缩的VCF(v). 使用-Ou能够让bcftools命令间的操作更快。

FILE：输入文件，可以是VCF或BCF，以及这些文件对应的BGZIP压缩形式。如果是-, 则认为是标准输入。有些工具需要tabix或CSI的索引文件。

mpileup和call是一套组合，最基本的用法为:

bcftools mpileup -Ou -f reference.fa alignments.bam | bcftools call -mv -Ob -o calls.bcf
# -m: 允许多倍体
# -v：表示只输出和基因组不同的位点

根据不同情况可以添加call的参数，比如说 bcftools call -P 1.1e-3。

view一般要配合-f LIST参数共同使用，能有效对VCF/BCF文件内容进行筛选。比如说仅选择非indel, 且ref的reads数小于1，深度在20和100之间。

bcftools view -i 'TYPE!="indel" && (DP4[0]+DP4[1])<1 && DP >20 && DP < 100'

有效表达式包括：

数值常量，字符常量和文件名: 1,1.0,1e-4，"string"，@filename
算术运算符：+, *, -, /
比较操作符：==, >, >=, <, <=, !=
正则表达式：INFO/HAYSTACK ~"needle/i"
括号： ()
逻辑运算符： &&, ||
INFO标签和FORMAT标签以及列名

INFO/DP 或 DP
FORMAT/DV, FMT/DV, 或 DV
FILTER, QUAL, ID, POS, REF, ALT[0]

1或0用于判断flag是否存在
"."则是判断是否有缺失值
样本基因型：纯合("hom")，杂合("het")，单倍体("hap")，alt-alt纯合("AA")，ref-alt杂合("RA")，alt-alt杂合("Aa")，单倍体参考("R")，单倍体替换("A")
REF/ALT列的变异类型(TYPE)： indel, snp, mnp, ref, bnd, other
数组下标： (DP[0]+DP[1])/(DP[2]+DP[3]) > 0.3。其中*表示任意，-表示返回
FORMAT和INFO标签的函数: MAX, MIN, AVG, SUM, STRLEN, ABS
运行过程中新增变量：N_ALT, N_SAMPLES, AC, MAC, AF, MAF, AN, N_MISSING, F_MISSING

query可以将VCF/BCF文件转换成更加人类可读的格式，依赖于-f FORMAT参数。

其中FORMAT可以是：

所有的列名：%CHROM, %POS, %ID, %REF, %ALT, %QUAL, %FILTER
INFO列的其中一个：%INFO/标签（如INFO/DP4），此外标签是多值结果，能用{}进一步选取，例如DP4{1}
基因型: "%GT", "%TGT"
换行符和制表符:"\n","\t"

实际操作

后续操作需要下载案例数据

curl -O http://data.biostarhandbook.com/variant/subset_hg19.vcf.gz
curl -O http://data.biostarhandbook.com/variant/subset_hg19.vcf.gz.tbi

从VCF中按照自定义格式提取数据

bcftools query -f '%CHROM %POS %REF %ALT \n' subset_hg19.vcf.gz | head -3
# 结果
19 400410 CA C
19 400666 G C
19 400742 C T

列出存放的所有样本

bcftools query -f subset_hg19.vcf.gz

从指定区域提取所有变异位点

bcftools query -f '19:400300-400800' -f '%CHROM\t%POS\t%REF%ALT\n' subset_hg19.vcf.gz | head -3
19  400410  CAC
19  400666  GC
19  400742  CT

这里是按照特定格式提取，如果希望输出也是VCF文件，则用filter或view命令。

从指定区域外提取所有变异

bcftools view -H -t ^'19:400300-400800' subset_hg19.vcf.gz | head -3
# 结果如下
19  400819  rs71335241  C   G   100 PASS    AC=0;AF=0.225839;AN=12;NS=2504;DP=10365;EAS_AF=0.2897;AMR_AF=0.2349;AFR_AF=0.2088;EUR_AF=0.161;SAS_AF=0.2434;AA=N|||;VT=SNP GT  0|0 0|0 0|0 0|0 0|0 0|0
19  400908  rs183189417 G   T   100 PASS    AC=1;AF=0.0632987;AN=12;NS=2504;DP=13162;EAS_AF=0.002;AMR_AF=0.1153;AFR_AF=0.0726;EUR_AF=0.0885;SAS_AF=0.0511;AA=-|||;VT=SNP    GT  0|0 0|0 0|0 0|0 0|0 0|1
19  400926  rs28420134  C   T   100 PASS    AC=1;AF=0.0259585;AN=12;NS=2504;DP=13731;EAS_AF=0.005;AMR_AF=0.0879;AFR_AF=0.003;EUR_AF=0.0457;SAS_AF=0.0143;AA=C|||;VT=SNP GT  0|0 0|0 0|0 0|0 0|1 0|0

根据样本中的基因型信息提取

# 通过表达式
bcftools view -e 'GT="." | GT="0|0"' subset_hg19.vcf.gz | bcftools query -f '%POS[\t%GT\t]\n' | head -3
402556  0|1     0|1     1|1     1|0     0|1     1|1
402707  0|1     0|1     1|1     1|0     0|1     1|1
402723  0|1     0|1     1|1     1|0     0|1     1|1
# 或者是-g/--genotype
## 选择至少有一个样本是杂合，且所有样本都不包含缺失位点信息
bcftools view -g het subset_hg19.vcf.gz | bcftools view -g ^miss | bcftools query -f '%POS[\t%GT]\n' | head -3

仅提取INDEL, 可用-v/--type或-i/--include, 当然这两者有细微区别。

bcftools view -v indels subset_hg19.vcf.gz | bcftools query -f '%POS\t%TYPE\n' | wc -l
bcftools view -i 'TYPE="indel"' subset_hg19.vcf.gz | bcftools query -f '%POS\t%TYPE\n' | wc -l

仅选择或不选择某几个样本

bcftools view -s HG00115,HG00118 subset_hg19.vcf.gz | bcftools query -H -f '%POS[\t%GT]\n' | head -n 4
bcftools view -s ^HG00115,HG00118 subset_hg19.vcf.gz | bcftools query -H -f '%POS[\t%GT]\n' | head -n 4

选择等位基因大于或者低于一定值的变异，即比较AC(alternate alleles count)

# 大于5
bcftools view  -c 5 subset_hg19.vcf.gz | bcftools query -f '%POS[\t%GT]\n' | head
## 结果如下
400666  1|0 0|1 0|1 0|0 0|0 1|1
401818  0|1 0|1 1|1 1|0 0|0 1|1
401907  0|1 0|1 1|0 1|0 0|0 0|1
# 低于5
bcftools view  -C 5 subset_hg19.vcf.gz | bcftools query -f '%POS[\t%GT]\n' | head -3
## 结果如下
400410  0|0 0|0 0|0 0|0 0|0 0|0
400666  1|0 0|1 0|1 0|0 0|0 1|1
400742  0|0 0|0 0|0 0|0 0|0 0|0

根据变异质量和覆盖深度选择

bcftools query -i 'QUAL>50 && DP>5000' -f '%POS\t%QUAL\t%DP\n' subset_hg19.vcf.gz | head -3
400410  100 7773
400666  100 8445
400742  100 15699

对于多个VCF文件，则需要用到merge和isec。

# 合并列表中的样本
bcftools merge -l samplelist > multi-sample.vcf
# 提取在所有样本都出现的变异
bcftools isec -p outdir -n=3 sample1.vcf.gz sample2.vcf.gz sample3.vcf.gz
# 提取至少在两个样本出现的变异
bcftools isec -p outdir -n+2 sample1.vcf.gz sample2.vcf.gz sample3.vcf.gz
# 提取仅仅在一个样本中出现的变异
bcftools isec -p outdir -C sample1.vcf.gz sample2.vcf.gz sample3.vcf.gz

输出结果就能直接导入到R,Python进行分析。

VCFtools

VCFtools: 用于描述性统计数据，计算数据，过滤数据以及数据格式转换。

基本用法：

vcftools [--vcf VCF文件 | --gzvcf gz压缩的VCF文件 --bcf BCF文件] [--out OUTPUT PREFFI]

他能做的事情：

输出第一条染色体的所有位点等位基因频率
从输入文件中仅保留SNP位点
输出两个vcf文件的比较结果
标准输出不含有filer tag的位点，并且以gzip压缩
计算每个位点的hardy-weinberg p-value，这些位点不包括缺失的基因型
计算一系列核酸多态性

常用参数如下：

和输入输出有关

--vcf, --gzvcf, --bcf：根据输入文件格式进行选择
--out, --stdout, -c --temp: 选择合适的输出方式

位点筛选(site filtering)有关参数

# 根据位置过滤
--chr/--not-chr Chr1: 选择染色体
--from-bp/--to-bp: 选择碱基范围
# 根据第三列ID进行过滤，不常用
--snp rsID
# 根据变异类型
--keep-only-indels: 仅保留INDEL
--remove-indels: 仅保留SNP
# 根据FILTER列进行过滤
--remove-filtered-all: 移除filter tag位点
--keep-filtered/--remove-filtered
# 根据INFO列进行过滤
--keep-INFO/--remove-INFO 目标类型
# 根据等位基因过滤
--maf/--max-maf # Minor Allele Frequency
--mac/--max-mac # Minor Allele Count

样本样本参数（我不常用）
基因型过滤参数（没怎么用）
输出VCF选项

--recode
--recode-INFO-all
--diff-site: 比较位点
--hardy: hardy-weinberg p值
--max-missing: 基因型缺失
--site-pi