bcftoolsのフィルタリングにおける条件式について

(本記事の内容はすべてをテストして調べたわけではないので、正しいことを保証できません。最終的には公式のマニュアルを参照されることをおすすめします。）

bcftoolsの多くのサブコマンドで使うことのできるフィルタリング用のオプションとして、 -i, --include と -e, --exclude があります。

これらのオプションは1つの引数をとり、その引数は EXPRESSION です。 EXPRESSION とはプログラミング用語で「式」のことです。式の詳細な説明はWikipediaなどにありますが、ここでは「決まった法則に従い、VCFの各レコードに対しTRUEまたはFALSEを返す条件式」という理解で十分かと思います。より厳密には各レコードに1つのTRUE/FALSEの返り値に加え、各サンプルに対してもTRUE/FALSEが返されます

もしTRUEが返されたレコードについては、 -i, --include を指定したときはそのレコードが出力され、-e, --exclude の場合はその逆となります。

文法

EXPRESSION にはあらかじめ決められた文法に従った式を記述する必要があります。この文法は多くのプログラミング言語の文法と大差がないので普段プログラミングする方にとっては簡単ですが、そうでない場合はややこしいかもしれません。

定数の書き方

定数とは条件式を評価する上でレコードに関わらず常に一定の値を示す具体的な値で、おもに閾値や検索対象文字列などを表現するために使われる部分です。以下に表記例を示します。

数字: 1, 1.0, 1e-4
文字列: "String"
ファイル名: @file_name

算術演算子

算術演算子としては +,*,-,/ の4つが使用可能です。意味は左からそれぞれ、加算、乗算、減算、除算です。

比較演算子

比較演算子としては、 == (= でもよい), >, >=, <=, <, != が使えます。意味は左からそれぞれ、イコール、大なり、大なりイコール、小なりイコール、小なり、ノットイコールです。多くのプログラミング言語では = は代入を表しますが、bcftoolsではイコールとして扱われます。

正規表現

正規表現も使うことができ、~ で正規表現にマッチしたらTRUEを返す式を表現できます。逆は !~ です。デフォルトではCase sensitiveなので、例えば、 CHROM ~ "Chr1" は Chr1 にはヒットしますが、CHR1 にはヒットしません。 Case insensitiveにするにはCHROM ~ "Chr1/i" とします。

かっこ

()が使えます。マニュアルには詳細がありませんでしたが、おそらくは計算順序を変えるために使用します。各演算子には優先順位があり、それを変えたい時に使うものだと思います。（演算子の優先順位についても特にマニュアルに言及はなさそうですが。）

論理演算子

論理演算子には &&, &, ||, | の4つがあります。&& 及び & は「かつ」、 || 及び | は「または」を表現します。

2つと1つの違いは片方または両方の被演算子がサンプルごとの結果を含むときに発生します。例えば FORMAT/DP>=10 や SMPL_MAX(FORMAT/AD)>=10 などです。逆に QUAL>=50 や MAX(INFO/AD)>10 などは1サイトにつき1つの結果のみを含み、サンプルに対する結果を含みません。要は FORMAT タグがらみかどうかです。

以下に簡単に各論理演算子の挙動をまとめました。siteにはそのサイト（レコード）がパスするかどうか（TRUE/FALSE）がスカラーの形で入っており、sampleにはサンプルごとのTRUE/FALSEがベクトル形式で入っているイメージで疑似コード的にまとめています。例えば QUAL>=50 の結果は site に、FORMAT/DP>=10 の結果は sample に入っているイメージです。（かなりややこしいのであまり自信はありません。）

`site & sample`

site = site and any(sample)
sample[i] = site and sample[i]

例：QUAL>=100 & FORMAT/DP>=30

QUALが100以上でかつ少なくとも1サンプルでDPが30以上の行でTRUE
QUALが100以上でかつ少なくともDPが30以上のサンプルでTRUE

`site && sample`

site & sample と同じ

`sample & sample`

site = any(sample and sample) ← 同一サンプルでTRUEであることが必要
sample[i] = sample[i] and sample[i] ← 同一サンプルでTRUEであることが必要

例：FORMAT/GQ>=50 & FORMAT/DP>=30

GQが50以上でかつDPが30以上のサンプルを少なくとも1つ含む行でTRUE
GQが50以上でかつDPが30以上のサンプルでTRUE

`sample && sample`

site = any(sample) and any(sample) ← いずれかのサンプルでTRUEであればOK
sample[i] = (any(sample) and any(sample)) and (sample[i] or sample[i]) ← そのサイトがパスする上で、サンプルについてはOR条件

個人的にはかなりややこしいので filter の -S, --set-GTs と使わない方がいいと思っています。

例：FORMAT/GQ>=50 && FORMAT/DP>=30

GQが50以上のサンプルが1以上、かつ、DPが30以上のサンプルが1つ以上の行でTRUE
上記条件が満たされれば、GQが50以上またはDPが30以上のサンプルはTRUE、満たさない場合は全サンプルでFALSE

`site | sample`

site = site or any(sample)
sample[i] = sample[i] ← もとのsiteがTRUE/FALSEに関係なくもとのsampleのコピー

例：QUAL>=100 | FORMAT/DP>=30

QUALが100以上または少なくとも1サンプルでDPが30以上の行でTRUE
DPが30以上のサンプルでTRUE

`site || sample`

site = site or any(sample)
sample[i] = site or sample[i] ← もとのサイトがTRUEなら全サンプルTRUE

例：QUAL>=100 || FORMAT/DP>=30

QUALが100以上または少なくとも1サンプルでDPが30以上の行でTRUE
QUALが100以上なら全サンプル、そうでないときはDPが30以上のサンプルでTRUE

`sample | sample`

site = any(sample) or any(sample)
sample[i] = sample[i] or sample[i] ← 同一サンプルに対してOR条件

例：FORMAT/GQ>=50 | FORMAT/DP>=30

GQが50以上のサンプルが1つ以上、または、DPが30以上のサンプルが1つ以上の行でTRUE
GQが50以上またはDPが30以上のサンプルでTRUE

`sample || sample`

site = any(sample) or any(sample)
sample[i] = any(sample[i]) or any(sample[i]) ← いずれかのサンプルでTRUEならそれがサンプル全体に波及

例：FORMAT/GQ>=50 || FORMAT/DP>=30

GQが50以上のサンプルが1つ以上、または、DPが30以上のサンプルが1つ以上の行でTRUE
GQが50以上のサンプルが1つ以上、または、DPが30以上のサンプルが1つ以上の行ならば全サンプルでTRUE、そうでなければ全サンプルでFALSE

条件を適用できる項目

CHROMなどの列、INFOタグ、FORMATタグについて条件を適用することができます。これらは文字列ではないのでダブルクオーテーションで囲む必要がありません。 INFO/ や FORMAT/ (FMT/) の部分は指している対象が明瞭であれば省略することができますが、個人的には省略しない方が後で意味を理解しやすいのでおすすめです。

INFO/DP or DP                                # INFO列のDPタグ
FORMAT/DV, FMT/DV, or DV                     # FORMAT列（サンプル列）のDVタグ(サンプル数分の要素を持つベクトルのイメージ）
FILTER, QUAL, ID, CHROM, POS, REF, ALT[0]　  # ALT[0] はALT列のコンマ区切りで1番目に記載されているアレル

Flag型のデータに対する条件式

Flag型のデータはそのタグが存在するかどうかしかありません。このタイプのタグでよく見るのは INFO/INDELだと思います。このタグが存在することを条件にするなら、 INFO/INDEL=1 、存在しないことを条件とするならINFO/INDEL=0 といった具合に 0,1で表現します。

欠損値の表現

欠損値は文字列 "." で表現されます。VCFの中身を眺めてみてもたまに値が . になっているときがあると思います。一番よく見るのはID列とALT列だと思います。 (例：-e 'ALT="."' でALT列が欠損の行を除く）

また、INFO/DP>=20 といったような条件式において、INFO/DP の実際の値が欠損値 . だった場合、結果は必ずFALSEになる模様です。（つまり、-i 'INFO/DP>=20 とした場合、結果には INFO列が DP=.; となっている行は含まれず、 -e INFO/DP>=20 とした場合は DP=.; となっている行は除かれません。）

ただし遺伝子型（ FORMAT/GT ）については欠損値は1種類ではない（マニュアルでは ".|.", "./.", ".", "0|." の4つが挙げられています。）ので、これらすべてにマッチする簡便な方法として GT="mis" とすることで欠損値にマッチさせることができます。そのほかの方法としては正規表現を用いて GT~"\." とする方法があります。ここで \ が入っているのは . が正規表現上特別な意味を持つのでそれをエスケープするためです。

遺伝子型に対する条件の便利な記法

遺伝子型(FORMAT/GT)については単なる文字列マッチや正規表現だけでなく、特別な意味を持つ文字列を使った条件式の立て方があり、以下に列記します。

GT="ref"  # REF型の半数体または二倍体（おそらく、0/0,0|0,0にマッチ）
GT="alt"  # alternate（おそらくGT="REF"でマッチしないかつ欠損でないものにマッチ。ヘテロにもマッチする。）
GT="mis"  # 欠損値にマッチ
GT="hom"  # ホモ接合体
GT="het"  # ヘテロ接合体
GT="hap"  # 半数体
GT="RR"   # REF-REFのホモ
GT="AA"   # ALT-ALTのホモ
GT="RA" or GT="AR"  # REF-ALTのヘテロ
GT="Aa" or GT="aA"  # ALT-ALTのヘテロ（おそらくALTが2つ以上のときに有効）
GT="R"    # REFの半数体
GT="A"    # ALTの半数体

TYPE変数

VCFのタグや列名として明示的に記されている訳ではないが使用できる変数として、 TYPEがあります。これはレコードのREFとALTを元にその変異がSNPなのかInDelなのかといった情報が格納されています。有効な値はindel,snp,mnp,ref,bnd,other,overlapの 7つのようです。refはおそらくREFのみでALTのデータがないレコードにマッチします。また、bndはこれもおそらくですが、リアレンジメントのレコードにマッチすると思います。

以下はリアレンジメントのレコード例

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
2   321681  .   G   G]17:198982]    .   .   .
17  198982  .   A   A]17:321681]    .   .   .

また、レコードによってはSNPかつInDelという場合もあります。そのため、例えば TYPE~"snp" とした場合は少なくともSNPであればマッチし、 TYPE="snp" とした場合はどのALTをみてもSNPでありそれ以外ではない場合にマッチします。

配列状のデータの要素へのアクセス

アレル頻度(INFO/AF)などタグによってはその中に複数のデータを内包している場合があります。（アレル頻度の場合は各レコードにアレルの数分のデータが含まれます。）また、FORMATタグにはサンプル分のデータが含まれます。さらにサンプルごとのアレル深度(FORMAT/AD)となるとサンプル x アレル数分のデータが含まれます。これらのデータについては何番目の要素に条件式を適用するかを明示的に指定することができます。 0-based なインデックスを[]内に入れることで表現できます。* を指定することですべての要素をとりだし、また、-を用いて範囲指定ができます。 :を使うことで何サンプル目の何番目のデータといったところまで表現できます。

また、複数のデータに対して条件式を適用すると、1レコードにつき複数のTRUEまたはFALSEが返されることになりますが、その場合は1つでもTRUEがあればそのレコードに対する返り値はTRUEになります。これはRに慣れている方であれば、any() 関数が暗黙のうちに呼ばれていると考えると理解しやすいかもしれません。(例：-i 'INFO/AF>0.1' とした場合は、INFO/AF内の複数データのうちいずれかが>0.1を満たしていればTRUE）

関数

関数として MAX, MIN, AVG, MEAN, MEDIAN, STDEV, SUM, STRLEN, ABS, COUNT が使えます。これらは、FORMATタグやINFOタグにおける複数データを含むものに対して適用することでそれらのデータの要約を取ることができます。また SMPL_MAX, SMPL_MIN, SMPL_AVG, SMPL_MEAN, SMPL_MEDIAN, SMPL_STDEV, SMPL_SUM という関数もあり、これらはサンプルごとに要約します。（結果サンプル数分の返り値を持ちます。）

binom test

binom() を使うことで二項検定ができるようです。おそらくですが、例えば'binom(FORMAT/AD)'はGTを参照して対応するアレルのADを引っ張ってきて、 p=0.5、x=AD_allele1、n=(AD_allele1+AD_allele2) の二項検定のp値を返すのだと思います。（この場合ホモの場合は必ずp値は1になると思います。）

必要に応じて計算される量

フィルタリングでよく使う統計量はVCF中に存在しなくても適宜計算されて利用することができます。利用できる変数は以下の通りです。

N_ALT, N_SAMPLES, AC, MAC, AF, MAF, AN, N_MISSING, F_MISSING, ILEN

N_PASSとF_PASS

これらは関数として用い、引数に EXPRESSION を指定することでそれぞれ条件式を満たすサンプル数、サンプル割合を計算してくれます。例えば N_PASS(FORMAT/GT="hom") はホモの数を返します。

その他注意点

マニュアルには

String comparisons and regular expressions are case-insensitive

とありますが、使ってみている感じcase-sensitiveであるように感じます。

文字列に,を含むときはそれらで分割される文字列に対してOR演算子で結合した形で一致判定が行われます。

EXPRESSION に含める文字は往々にしてShellに対しても特別な意味を持つので、 Shellに解釈されないように常にシングルクオーテーションで囲むことが推奨されます。

使用例

# QUALが100以上の行を抽出
bcftools view -i 'QUAL>=100' input.vcf
# FORMAT/DPが30以上かつFORMAT/GQが50以上を残し、条件を満たさないサンプルを欠損とする
bcftools filter -S. -i 'FORMAT/DP>=30 & FORMAT/GQ>=50' input.vcf
# ホモを示すサンプルが一つもない行を抽出
bcftools view -i 'N_PASS(GT="hom")=0' input.vcf
# MAFが0.05以上かつ欠損率が0.1未満の行を抽出
bcftools view -i 'MAF>=0.05 && F_MISSING<0.1' input.vcf
# MAFが0.05未満の行に対しFILTER列に"lowmaf"と記入
bcftools filter -s lowmaf -e 'MAF<0.05' input.vcf
# 1番目と2番目のサンプルがホモでかつ異なる遺伝子型を示す行を抽出
bcftools view -i '(GT[0]="RR" && GT[1]="AA") || GT[0]="AA" && GT[1]="RR")' input.vcf
# biallelicなSNPの行を抽出
bcftools view -i 'N_ALT=1 && TYPE="snp"' input.vcf
# PLの最小が0より大きいサンプルが1つでも存在する行を除く
bcftools view -e 'SMPL_MIN(FORMAT/PL)>0' input.vcf

（蛇足）

EXPRESSION を使わないフィルタリング

# biallelicなSNPの行を抽出
bcftools view -m2 -M2 -v snp input.vcf
# 近傍1kbにInDelを含まない行を抽出
bcftools filter -g 1000 input.vcf
# 指定した領域に含まれる行を抽出
bcftools view -t 1:1-10000 input.vcf
# Minor Allele Frequencyが0.05以上の行を抽出
bcftools view -q 0.05:minor input.vcf
# 全サンプルでコールされていない行を除く
bcftools view -U input.vcf

・・・