次世代シーケンシング技術（NGS）は、ゲノム上に存在するほぼすべての遺伝子バリアントの発見を可能にする。しかし、これらのバリアントの一部は、NGSやバリアントコーラーの限界により、シーケンスの質が低い場合がある。多数のシーケンスされた個体を解析する遺伝学的研究では、偽の発見の原因となる可能性があるため、品質の悪いバリアントを検出して除去することが重要である。本論文では、従来のフィルタリングアプローチと機械学習アプローチを組み合わせて、NGSデータから同定されたバリアントの品質管理を行う統計ツール「ForestQC」を紹介する。本ソフトウェアは、シークエンシングデプス、ジェノタイピングの品質、GC含有量などのシークエンス品質に関する情報を用いて、特定のバリアントが偽陽性である可能性が高いかどうかを予測する。ForestQCを評価するために、2つの全ゲノムシーケンスデータセットに適用した。1つのデータセットは家族からなる血縁者で構成され、もう1つのデータセットは血縁者ではない個人で構成されている。その結果、ForestQCは、GATKのVQSRなど、バリアントの品質管理を行うために広く用いられている手法よりも、解析に含めるべきバリアントの品質を大幅に向上させることができた。また、ForestQCは非常に効率的であるため、大規模なシーケンスデータセットに適用することができる。シーケンスの品質情報を用いて学習した機械学習アルゴリズムとフィルタリング手法を組み合わせることは、シーケンスデータから遺伝的バリアントの品質管理を行うための実用的なアプローチであると結論づけた。
ForestQCは、生のVCFファイルを入力とし、高品質または低品質のバリアントを決定する。本手法は、事前に定義された一連のフィルターによって高品質なバリアントと低品質なバリアントを決定するフィルターアプローチと、機械学習を用いてバリアントが高品質か低品質かを分類する分類アプローチを組み合わせている。論文図1に示すように、本手法ではまず、GWASでQCを行う際に一般的に使用されるいくつかのフィルターについて、各バリアントの統計量を計算する。これらの統計値は、ABHet、HWE p-value、遺伝子型欠損率、家族ベースのデータセットのメンデルエラー率、および任意のユーザー定義の統計値で構成されている（詳細はMaterials and Methodsに記載）。ForestQCは、これらの統計情報をフィルターとして使用して、3つのバリアントのセットを特定する。1) すべてのフィルターを通過した高品質なバリアントのセット、2) いずれかのフィルターに失敗した低品質なバリアントのセット、3) 高品質でも低品質でもない不確定なバリアントのセット。フィルターには厳格なしきい値を使用しており（S1およびS2表）、高品質のバリアントは品質が良く、低品質のバリアントは実際に偽陽性であるか、または明らかにシーケンス品質が悪いことを強く確信している。ForestQCの次のステップは、フィルタリングステップで検出した高品質バリアントと低品質バリアントを使って、ランダムフォレスト機械学習モデルをトレーニングすることである。ForestQCでは、高品質なバリアントと低品質なバリアントの7つのシーケンス品質メトリクスがランダムフォレストモデルを学習するための特徴として使用されている。その中には、シークエンシングデプスに関連する3つの指標、遺伝子型の品質に関連する3つの指標、およびGC含有量に関連する1つの指標が含まれる。最終的に、モデルは各未決定バリアントが高品質であるか低品質であるかを予測する。ランダムフォレスト分類法で予測された高品質バリアントと、フィルタリングステップで検出された高品質バリアントを組み合わせて、ForestQCで決定された高品質バリアントの完全なセットとする。低品質のバリアントの識別にも同じ手順が適用される。
インストール
ubuntuでcondaの仮想環境を使って導入した（python3.9）。
mamba create -n forestqc_env -y
conda activate forestqc_env
mamba install forestqc -c avallonking -y
#pandasとscikit-learnが無かったので追加で導入した
> ForestQC -h
$ ForestQC -h
ForestQC v1.1.5.4 by Jae Hoon Sul Lab at UCLA
--Quality control on genetic variants from next-generation sequencing data using random forest
usage: ForestQC [-h] {stat,split,classify,set_outlier,compute_gc} ...
positional arguments:
{stat,split,classify,set_outlier,compute_gc}
ForestQC command help
stat variants statistics help
split variants splitting help
classify variants classification help
set_outlier set outlier help
compute_gc compute GC content help
optional arguments:
-h, --help show this help message and exit
実行方法
１、hg19以外のリファレンスゲノムを使用している場合、最初にGCコンテンツテーブルを生成する（レポジトリのexamples/gc_content_hg19.tsvにhg19のGCコンテンツテーブルが用意されている）。
ForestQC compute_gc -i genome.fasta -o GC_table
（GCの計算は数分で終わる）
２、vcfファイルから統計量を計算する必要がある。各バリアントのすべての統計情報を含むファイルが出力される。
ForestQC stat -i input.vcf -o stat_out -c GC_table
vcfファイルには、ゲノムの全ての領域からのバリアントが含まれている必要がある。また、VCFには、それぞれのサイトについて各個人の情報が含まれている必要がある（ほかにも家系情報を提供しない場合の制限などがある。詳細はGithub参照）。
３、データセットを高品質、低品質、"undertermined" の各バリアントに分ける
ForestQC split -i stat_out -o split_out
出力ファイルはgood<suffix>, bad<suffix>, gray<suffix>の3つのモデルファイルになる。
４、ランダムフォレストモデルを訓練して、グレーバリアントの仕分けに適用する。
ForestQC classify -g good.xxx -b bad.xxx -y gray.xxx
出力ファイルは predicted_good_<suffix> と predictred_bad_<suffix>となる（”-o”でファイル名指定）。
５、最終的な良いバリアントと悪いバリアントのセットを得るために、ランダムフォレストモデルによって予測されたバリアント（ステップ4）と、フィルターによって分離されたバリアント（ステップ3）を組み合わせる
cat good.ForestQC_split predicted_good.variants.tsv > total_good_variants
cat bad.ForestQC_split predicted_bad.variants.tsv > total_bad_variants
６、最終的なForestQCをパスしたバリアントを含むVCFファイルを得るため、badのバリアントの位置を取得し、元のVCFファイルからそれらを取り除く。awkで特定のフィールドだけ取り出し、vcftoolsの指定されたポジションと異なるポジションだけ出力する機能を使う。
awk -F "\t" 'NR>1{print $2"\t"$3}' total_bad_variants > bad_variants_positions
vcftools --gzvcf [original_gziped_vcf_file] --exclude-positions [bad_variants_positions] --recode --recode-INFO-all -c | gzip -c > [output_vcf]
（レポジトリではchrのヘッダ名が1や2などの場合にsedに渡す例あり）
試したデータ（GATKのベストプラクティスガイドラインに従ってコールしたraw vcf）では、95％以上のポジションがフィルタリングされました。
引用
ForestQC: Quality control on genetic variants from next-generation sequencing data using random forest
Jiajin Li,Brandon Jew,Lingyu Zhan,Sungoo Hwang,Giovanni Coppola,Nelson B. Freimer,Jae Hoon Sul
PLoS Comput Biol. 2019 Dec 18;15(12):e1007556
関連
