2021 9/29 タイトル修正(変異 => バリアント), help追加
以前SnpEffというバリアントのアノテーションを行うことができるツールを紹介した(リンク)。このツールにはもう一つSnpSiftというツールが同梱されている。SnpSiftはバリアントコール結果のVCFファイルを扱うツールで、クオリティやp値など様々な指標に基づいて変異を仕分けることができる。大きなゲノムの解析では数万を超えるバリアントがコールされてくるのが普通だが、そういった時、統計の理屈上偽陽性は一定数出てきてしまう。そこで、はじめに信頼度などの条件に従って変異コール結果をフィルタリングできれば便利である。
SnpSiftはそんな時に使えるツールです。使い方をまとめてみます。
公式ページ
インストール
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
3つ以上のサンプルで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の標準フォーマットに従っていれば、色々な取り出し方が可能である。例えば下のように
カラムの説明で#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Zがコールの1行前につけられている。上の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フィールドでフィルタリングする。
2サンプルのコール結果なら、例えば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で管理されているので、2番サンプルならGEN[1]。全サンプルならGEN[*]で表せる。
PLのように複数フィールドを持っているなら(上の例の1番サンプルなら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
Intersect
RmInfo
GeneSets
GT
VcfCheck
似たツールにvcflibがある。
https://github.com/vcflib/vcflib
引用
Using Drosophila melanogaster as a model for genotoxic chemical mutational studies with a new program, SnpSift
Cingolani P, Patel VM, Coon M, Nguyen T, Land SJ, Ruden DM, Lu X
Front Genet. 2012 Mar 15;3:35.