バリアントのフィルタリングを行うSnpSift

2021 9/29 タイトル修正（変異　＝＞　バリアント）, help追加

以前SnpEffというバリアントのアノテーションを行うことができるツールを紹介した（リンク）。このツールにはもう一つSnpSiftというツールが同梱されている。SnpSiftはバリアントコール結果のVCFファイルを扱うツールで、クオリティやp値など様々な指標に基づいて変異を仕分けることができる。大きなゲノムの解析では数万を超えるバリアントがコールされてくるのが普通だが、そういった時、統計の理屈上偽陽性は一定数出てきてしまう。そこで、はじめに信頼度などの条件に従って変異コール結果をフィルタリングできれば便利である。

SnpSiftはそんな時に使えるツールです。使い方をまとめてみます。

公式ページ

SnpSift

インストール

SnpEff

SnpEffのパッケージに入っているので、SnpEffをインストールしている人はすでに持っているはず。SnpEffはbrewでも導入できます。

VCFフォーマットについて簡単にまとめています。

#bioconda (link)
mamba install -c bioconda snpsift -y

> SnpSift

# SnpSift

SnpSift version 4.3t (build 2017-11-24 10:18), by Pablo Cingolani

Usage: java -jar SnpSift.jar [command] params...

Command is one of:

alleleMat : Create an allele matrix output.

annotate : Annotate 'ID' from a database (e.g. dbSnp). Assumes entries are sorted.

caseControl : Compare how many variants are in 'case' and in 'control' groups; calculate p-values.

ccs : Case control summary. Case and control summaries by region, allele frequency and variant's functional effect.

concordance : Concordance metrics between two VCF files.

covMat : Create an covariance matrix output (allele matrix as input).

dbnsfp : Annotate with multiple entries from dbNSFP.

extractFields : Extract fields from VCF file into tab separated format.

filter : Filter using arbitrary expressions

geneSets : Annotate using MSigDb gene sets (MSigDb includes: GO, KEGG, Reactome, BioCarta, etc.)

gt : Add Genotype to INFO fields and remove genotype fields when possible.

gtfilter : Filter genotype using arbitrary expressions.

gwasCat : Annotate using GWAS catalog

hwe : Calculate Hardy-Weimberg parameters and perform a godness of fit test.

intersect : Intersect intervals (genomic regions).

intervals : Keep variants that intersect with intervals.

intIdx : Keep variants that intersect with intervals. Index-based method: Used for large VCF file and a few intervals to retrieve

join : Join files by genomic region.

op : Annotate using an operator.

phastCons : Annotate using conservation scores (phastCons).

private : Annotate if a variant is private to a family or group.

rmRefGen : Remove reference genotypes.

rmInfo : Remove INFO fields.

sort : Sort VCF file/s (if multiple input VCFs, merge and sort).

split : Split VCF by chromosome.

tstv : Calculate transiton to transversion ratio.

varType : Annotate variant type (SNP,MNP,INS,DEL or MIXED).

vcfCheck : Check that VCF file is well formed.

vcf2tped : Convert VCF to TPED.

Options common to all SnpSift commands:

-d : Debug.

-download : Download database, if not available locally. Default: true.

-noDownload : Do not download a database, if not available locally.

-noLog : Do not report usage statistics to server.

-h : Help.

-v : Verbose.

実行方法

filter - 条件に従ってフィルタリングする。

quality30以上のコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter " ( QUAL >= 30 )" > filtered.vcf

quality20以上のindelとquality30以上のコールを取り出す。(パイプ"|"でA or Bを表す）

cat variants.vcf | java -jar SnpSift.jar filter "(( exists INDEL ) & (QUAL >= 20)) | (QUAL >= 30 )" > filtered.vcf

カバレッジ (DP) 25以上のコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "(DP >= 25)" > filtered.vcf

３つ以上のサンプルでhomozygoteのコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "(countHom() > 3)" > filtered.vcf

heterozygoteでカバレッジ25以上のコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "((countHet() > 0) && (DP >= 25))" > filtered.vcf

VCFの標準フォーマットに従っていれば、色々な取り出し方が可能である。例えば下のように

f:id:kazumaxneo:20170819204540j:plain

カラムの説明で#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Zがコールの１行前につけられている。上のVCFからCHROMフィールドが"chr"だけ取り出すなら

cat variants.vcf | java -jar SnpSift.jar filter "( CHROM = 'chr' )" > filtered.vcf

フィルターがPASSしているものだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "(FILTER = 'PASS')"  > filtered.vcf

Positionが 123456-654321のコールだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter "( POS > 123456 ) & ( POS < 654321 )" > filtered.vcf

-n オプションを使い、chrとマッチしないものだけ取り出す。

cat variants.vcf | java -jar SnpSift.jar filter -n "( CHROM = 'chr' )" > filtered.vcf

-n　Inverse. Show lines that do not match filter expression

genotypeフィールドでフィルタリングする。

２サンプルのコール結果なら、例えばgenotypeフィールドは

"GT:PL:GQ 1/1:255,66,0:63 0/1:245,0,255:99"

のようになっている。1番サンプルがGQ>60以上の部位だけ取り出すなら

cat variants.vcf | java -jar SnpSift.jar filter "( GEN[0].GQ > 60 )"  > filtered.vcf

サンプルはindexで管理されているので、２番サンプルならGEN[1]。全サンプルならGEN[*]で表せる。

PLのように複数フィールドを持っているなら（上の例の１番サンプルならPLは255,66,0）、PLの後にも[index]をつける。

cat variants.vcf | java -jar SnpSift.jar filter "( GEN[0].PL[2] = 0 )"  > filtered.vcf

他にも多種多様なフィルタリングが可能になっている。公式に膨大な例があるので、そちらを確認してみてください。

Annotate - 別のVCFファイルのアノテーション情報を移す。

dbSNPや1000ゲノムなど大規模データのアノテーション情報を移したりするのに使う。

SnpSift annotate dbSnp132.vcf.gz variants.vcf.gz > variants_annotated.vcf

例えばアノテーション情報を移す前のVCFが以下なら、

22      16346045    .            T    C    0.0    FAIL    NS=244

変異が参照先のVCFと合致すると

22      16346045    rs56234788   T    C    0.0    FAIL    NS=244;RSPOS=16346045;GMAF=0.162248628884826;dbSNPBuildID=129;SSR=0;SAO=0;VP=050100000000000100000100;WGT=0;VC=SNV;SLO;GNO

のようにIDフィールドとINFOフィールドに追記される。-idオプションをつけると、IDフィールドの書き込みだけ行われる。

他にも様々なコマンドが用意されている。時間を見つけて順番にまとめていきます。

CaseControl

RmRefGen

TsTv

ExtractFields

varType

gwasCat

Split

Concordance

Private

Vcf2Tped