使用SnpSift filter对VCF文件进行筛选

2020-05-11 10:31:04 浏览数 (1)

当完成突变位点注释之后，我们会得到一个巨大的VCF文件，文件大小从几十M到几十G不等。在数量如此多的突变位点中，我们只会根据注释结果从中挑选部分感兴趣的突变位点，这就要求对VCF文件进行过滤。如此大的文件用Excel 操作是不现实的，脚本语言处理大文件时效果也不尽人意，所以SnpEff的开发团队专门开发了一款工具，叫做SnpSift, 用来对VCF文件进行过滤。

对VCF进行过滤是SnpSift最基本的功能，除此之外，还提供了许多非常实用的功能。软件是集成在SnpEff软件包中的， 在下载的SnpEff安装包中，就会有这款软件。

基本用法如下

代码语言：javascript复制
cat variants.vcf | java -jar SnpSift.jar filter "(CHROM = 'chr1')" > filtered.vcf

上面的例子是删选位于1号染色体上的突变位点， CHROM = 'chr1'就是筛选的条件表达式，用双引号包裹起来。`snpsift 提供了灵活的筛选策略，最核心的地方，就是条件表达式的写法，有以下几种用法

1. 单一值的变量

在VCF文件中，常用的字段都可以作为变量名称

1. CHROM
代码语言：javascript复制
"( CHROM = 'chr1' )"
2. POS
代码语言：javascript复制
"( POS > 123456 ) & ( POS < 654321 )"
3. REF
代码语言：javascript复制
"( REF = 'A' )"
4. ALT
代码语言：javascript复制
"( ALT = 'T' )"
5. QUAL
代码语言：javascript复制
"( QUAL > 30 )"
6. FILTER
代码语言：javascript复制
"( na FILTER ) | (FILTER = 'PASS')"
7. INFO
代码语言：javascript复制
( DP > 10 ) & ( AF1 = 0 )
2. 多个值的变量

当一个变量有多个值时，可以使用下标进行访问，下标从0开始，比如CI95=0.04167,0.5417，对这个变量的值进行筛选可以采用如下的写法

代码语言：javascript复制
"( CI95[0] > 0.1 ) & (CI95[1] <= 0.3)"

也可以采用通配符一次匹配多个下标，ANY和*的作用相同，至少一个值符合要求即可，写法如下

代码语言：javascript复制
"( CI95[ANY] > 0.1 )"
"( CI95[*] > 0.1 )"

ALL和?的作用相同，所有值都符合要求，写法如下

代码语言：javascript复制
"( CI95[ALL] > 0.1 )"
"( CI95[?] > 0.1 )"
3. 基因型信息

在表示基因型信息时，会看到如下的字段

代码语言：javascript复制
GT:PL:GQ 1/1:255,66,0:63 0/1:245,0,255:99

代表了两个样本的基因型相关信息，在筛选时，对于不同的样本即可以采用数字下标表示，也可以用样本名称表示

代码语言：javascript复制
"( GEN[0].GQ > 60 )"
"( GEN[HG00097].GQ > 60 )"

用GEN特指基因型信息，在对应样本时，可以用数字下标，也可以直接用样本名称，对于每个样本，会有GT, GQ, BL等信息，和样本之间用.点号连接。对于多个样本的基因型，也可以使用上面提到的通配符进行筛选。

4. 集合

当筛选一个变量的多种取值时，可以采用|逻辑或操作符进行连接，比如下面的示例，筛选位于1到3号染色体上的突变位点

代码语言：javascript复制
"((CHROM = 'chr1') |  (CHROM = 'chr2') | (CHROM = 'chr3'))"

但是当取值很多时，比如5个，10个，甚至更多，如果用|符号连接，对应的表达式就会特别的冗长，此时可以将这些值写到一个文件中，作为一个集合，每一行是一个值

代码语言：javascript复制
chr1
chr2
chr3

在筛选时，通过-s参数指定这个文件，用法如下

cat variants.vcf | java -jar SnpSift.jar filter -s set.txt > filtered.vcf

5. 操作符

在书写条件表达式时，用到的操作符如下

可以分成如下几类

1. 比较操作符

=, >, >=, <, <=

2. 正则匹配

=~和!~

3. 逻辑操作符

&,|,!

4. 特殊操作符

exists判断某个字段的值是否存在；has表示包含的意思，有些字段会有多个取值的情况，比如snpEFF在注释突变类型时，会有以下结果

代码语言：javascript复制
intron_variant&nc_transcript_variant

代表突变位点属于这两种突变类型，对这种突变类型过滤时，采用has操作符更加合适，写法如下

代码语言：javascript复制
"ANN[*].EFFECT has 'missense_variant'"
6. 函数

函数主要用来操作基因型信息，详细的列表如下

对于基因型，提出了以下4个概念

  1. homozygous
  2. heteroygous
  3. variant
  4. reference

homozygous代表纯合，比如0/0和1/1; heteroygous代表杂合，比如0/1;variant代表发生了突变，包括0/1, 1/1; reference代表没有发生突变，专指0/0;函数可以分成两大类型 第一种类型判断基因型的类型，is开头，包括以下4种

  1. isHom
  2. isHet
  3. isVariant
  4. isRef

第二种类型统计特定基因型的个数，`count开头，包括以下4种

  1. countHom
  2. countHet
  3. countVariant
  4. countRef

以上的规则对于所有的VCF文件都通用，对于SnpEff注释产生的VCF文件，包含以下几个特殊字段的信息

1. ANN

snpEff 注释的结果存储在ANN字段中，ANN字段中每个sub field对应的变量名称如下

  • ALLELE (alias GENOTYPE)
  • EFFECT (alias ANNOTATION):
  • GENE: Gene name
  • GENEID: Gene ID
  • FEATURE
  • FEATUREID (alias TRID: Transcript ID)
  • BIOTYPE:
  • RANK
  • HGVS_C
  • HGVS_P
  • CDNA_POS (alias POS_CDNA)
  • CDNA_LEN (alias LEN_CDNA)
  • CDS_POS (alias POS_CDS)
  • CDS_LEN (alias LEN_CDS)
  • AA_POS (alias POS_AA)
  • AA_LEN (alias LEN_AA)
  • DISTANCE
  • ERRORS (alias WARNING, INFOS)

条件表达式的写法如下：

代码语言：javascript复制
"ANN[*].EFFECT has 'missense_variant'"

由于ANN字段可以由多个feature构成，所以也是采用了数组的形式

2. EFF

旧版本的snpEff会将注释信息保存在EFF字段中， 每个sub field对应的变量名称如下

  • EFFECT
  • IMPACT
  • FUNCLASS
  • CODON
  • AA
  • GENE
  • BIOTYPE
  • CODING
  • TRID
  • RANK

条件表达式的写法如下：

代码语言：javascript复制
"( EFF[0].EFFECT = 'NON_SYNONYMOUS_CODING' )"
3. LOF 和 NMD

这两个字段的变量相似，有以下几种

  • LOF.GENE | NMD.GENE
  • LOF.GENEID | NMD.GENEID
  • LOF.NUMTR | NMD.NUMTR
  • LOF.PERC | NMD.PERC

条件表达式的写法如下

代码语言：javascript复制
"(exists LOF[*].PERC) & (LOF[*].PERC > 0.9)"

通过snpsift, 可以方便灵活的对VCF文件进行过滤，在后续的文章中，还会介绍更多该软件实用的功能。

cdn 编程算法

0 人点赞