NGSデータから取得した変異情報(VCF形式)をGWAS解析のPLINKに読み込みます。
SNPアレイデータを解析に使うことが多いかなと思いますが、
試してみたのでご紹介します。
作業環境
- CentOS 7
- bcftools 1.8
- bgzip (htslib) 1.9
- PLINK v1.90b6.7 64-bit
やり方
２ステップになります。
- 検体ごとに分かれたvcfを1つに統合
- PLINKでバイナリ化
1. 検体ごとに分かれたvcfを1つに統合
1つのファイルに多サンプルのデータが入った、multi sample VCFをつくります。
まとめたいそれぞれのVCFは 圧縮、インデックス作成をしておきましょう。
インデックスといっても、tabix ではなく bcftools index のほうです。 *.csi というファイルができます。
bcftools view -O z -o sampleA.vcf.gz sampleA.vcf #または bgzip sampleA.vcf bcftools index sampleA.vcf.gz
次に
統合したいファイル名（or ファイルへのパス）が書かれたテキストファイルを作ります。
### "vcflist.txt" ### sampleA.vcf.gz sampleB.vcf.gz :
統合します。
bcftools merge -l vcflist.txt -m both -O z -o <出力ファイル名>.vcf.gz
1. のおまけ
統合しようとするとたまにエラーで怒られることですが、
もし重複するサンプル名があるという場合は、そのままでは統合ができないので
該当のファイルは取り除くか、あらかじめbcftoolsでサンプル名を変えておきましょう。
### サンプル名が書かれたテキストファイル ### newname
bcftools reheader -s <サンプル名が書かれたテキストファイル> sampleA_2.vcf -o <出力ファイル名>.vcf
または統合した後で一気に全部のサンプル名を変えるのもありです。
bcftools merge --force-samples -l vcflist.txt -m both -O z -o temp.vcf.gz bcftools reheader -s <サンプル名のテキストファイル> -o <出力ファイル名>.vcf temp.vcf.gz
2. PLINKでバイナリ化
解析したい表現型(case/control や量的形質)が書かれたテキストファイルも用意しましょう。
### 表現型のテキストファイル例 ### FamilyA sampleA 1 100 FamilyB sampleB 0 90 FamilyC sampleC 0 80 ：
↑３列目以降がphenotype の値となる。1-2行目はfamily ID と個体ID。血縁関係がなければfamily ID は個体IDと同じでよいだろう。
では読み込みし、PLINKに必要なバイナリファイルを作ります。
VCFは圧縮/非圧縮ともに読み込めるようです。
また、ALT が複数ある行は最初の１個け読み込まれるようです。
plink --make-bed --vcf <入力ファイル>.vcf.gz --maf 0.05 \
--pheno <表現型のテキストファイル> --out <出力先>
・--maf: マイナーアリル頻度がこれ未満の変異は除外する
・--pheno ：スペースorタブ区切りのテキストファイル。３列目以降がphenotype の値となる。1-2行目はfamily ID と個体ID
・--out <出力先>.<拡張子> というファイルができる。 なのでここはdirectory/prefix の形で書くと良い。
このほかにもオプションで、VCFのフィルターや変異クオリティを用いたフィルタリングを行えます。
実際に行うときはデータに合わせてオプションを検討しましょう。*1
このようなファイルが出来ます。
$ ls directory prefix.bed prefix.bim prefix.fam prefix.log prefix.nosex
このあとは SNPアレイデータを読み込んだのと同じように解析ができそうです。
