eigenstrat


※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

eigenstratは、principal components analysis(主成分解析)の手法をもちいて集団の構造化を解析し、関連検定の結果を補正するソフトウェア。

convertfプログラムでファイル形式をEIGENSTRAT形式とし、smartpca.perlプログラムで構造化の解析を図示するところまで解説する。

ファイル構成

おおまかにはbinとEIGENSTRATとCONVERTFのディレクトリが重要。作業はCONVERTFまたはEIGENSTRATのディレクトリ内で行う。まずはCONVERTFディレクトリに移動しておく。

cd CONVERTF

convertf

入力ファイル形式は5つが許されている。linkage関連のソフトの標準的な入力フォーマットであり、Haploviewの入力にも使用可能なPED形式がやりやすいと思われる。convertfプログラムを用いて、PED形式のファイルをEIGENSTRAT形式に変換する。

PED形式

PED形式は三つのファイルから構成される。遺伝子型データは.pedファイル、SNP座位情報(chromosome locationなど)は.pedsnpファイル、サンプル個体データは.pedindファイルに入力する。

これらのファイルの具体的な例を見たいときには、CONVERTFディレクトリ内にあるexample.ped、example.pedsnp、example.pedindファイルをみればよい。すなわち

more example.ped

など。

.pedファイル

ダウンロードされたファイルのCONVERTFファイルに入っているexample.pedファイルは次のようになっている。一行目はここでの説明用の行で、ファイルに含まれている情報ではない。

家族ID サンプルID 親ID1 親ID2 性別 表現型 optional
1 SAMPLE0 0 0 2 2 1 1 2 1 1 2 2 1 1 2 2 1 1 2 2
2 SAMPLE1 0 0 1 2 1 1 2 1 2 1 2 1 1 1 2 1 1 2 2
3 SAMPLE2 0 0 2 1 1 1 2 2 2 1 2 1 2 1 2 1 2 1 2
4 SAMPLE3 0 0 1 1 1 1 1 1 2 1 1 2 2 1 1 1 2 1 2
5 SAMPLE4 0 0 2 1 1 1 1 2 2 1 2 2 2 1 1 1 2 1 1
  • EIGENSTRATはunrelated samplesを扱うソフトなので、1,3,4列目は無視される。なんでもいい。
  • 性別は、男が1、女が2
  • 表現型は、コントロールが1、ケースが2 (PEDファイルでは量的形質やカテゴリを表すこともある)
  • そのあとのSNP情報は、0ACGTを01234とする。ここで0はミッシングデータ。

.pedsnpファイル

example.pedsnpファイルは次のようである。例によって一行目はファイルに含まれていない解説用の行。

染色体 SNP名 genetic position (Morgan) physical position (bp) reference (optional) variant alleles(optional)
11 rs0000 0.000000 0
11 rs1111 0.001000 100000
11 rs2222 0.002000 200000
11 rs3333 0.003000 300000
11 rs4444 0.004000 400000
11 rs5555 0.005000 500000
11 rs6666 0.006000 600000
  • X染色体はXを入力する。

.pedindファイル

.pedファイルの最初の6-7列目が入力されたもの。

convertfを使用

CONVERTFディレクトリに移動する。コマンドは、ダウンロード時についてくる例を使えば次のようになる。

../bin/convertf -p par.PED.EIGENSTRAT

-pのあとの引数が命令のファイルとなる。このファイルpar.PED.EIGENSTRATの中身は次のようなもの

genotypename:    example.ped
snpname:         example.pedsnp # or example.map, either works 
indivname:       example.pedind # or example.ped, either works
outputformat:    EIGENSTRAT
genotypeoutname: example.eigenstratgeno
snpoutname:      example.snp
indivoutname:    example.ind

よくわからないが、実際に使うときはこれをいじくればいいのではないかと思われる。たとえばtest.ped,...などとすべてtestという名前をつけたなら、

vi par.PED.EIGENSTRAT

にてテキストエディタのviを起動し、

:%s/example/test/g

としてすべて置換すればいいだろう。もとのpar.PED.EIGENSTRATを残しておきたいなら

cp par.PED.EIGENSTRAT par.test

などとしてコピーをとり、それを編集すればよい。

これでPED形式をEIGENSTRAT形式に変換することができた。

smartpca, eigenstrat

smartpca.perlは主成分解析を行い、遺伝子型行列の分散・共分散行列から固有値の高い順に固有ベクトル(=主成分)を計算し、固有ベクトルと各個体の遺伝子型データの内積をプロットすることにより集団の構造化を検出するものである。

・・・といってもよくわからない。コマンドの例が用意されている。EIGENSTRATディレクトリに移動しよう。次のように入力する。

example.perl

そうすると次のようなコマンドが実行される。

smartpca.perl -i example.geno  -a example.snp  -b example.ind -k 2  -o example.pca  -p example.plot  -e example.eval  
-l example.log  -m 5  -t 2  -s 6.0
  • -i example.geno は、EIGENSTRAT形式の入力ファイルである。
  • -a example.snp は、入力するsnpファイル
  • -b example.ind は、入力するindivファイル
  • -k 2 は、出力する主成分(principal components)の数で、デフォルトは10
  • -o example.pca は、主成分の出力ファイル。除外値として除外された個体は、このファイルの中で0. 0と表される。
  • -p example.plot は、トップ二つの主成分を用いた出力プロットファイルである。
  • -e example.eval は、全ての固有値の出力ファイル
  • -l example.log は、出力のログファイル
  • -m 5 は、除外値さがしの繰り返し数で、デフォルトは5。除外をしない場合は-m 0とする。
  • -t 2 は、除外値さがしのときに用いる主成分の数で、デフォルトは10
  • -s 6.0 は、除外値の決定に用いる標準偏差。デフォルトは6.0

そのあと次の命令が実行されるが、これはsmartpca.perlが呼び出しているようだ。これがsmartpcaのコアとなるプログラムのようだ。

smartpca -p example.pca.par >example.log

その次の命令もsmartpca.perlの呼び出しのようだ。これはトップの二つの主成分からのプロットを表示するperlプログラムなのだそうだ。

ploteig -i example.pca.evec -c 1:2 -p Case:Control -x -o example.plot.xtxt

その次がようやくお出ましのeigenstratプログラム。

eigenstrat -i example.geno  -j example.pheno  -p example.pca  -l 1  -o example.chisq
  • -i example.geno は、EIGENSTRAT形式の入力ファイル。
  • -j example.pheno は、表現型の入力ファイルのようだが、下の命令によってPED形式のindファイルから変換して作成できるようだ。
../CONVERTF/ind2pheno.perl
  • -p example.pca は、smartpca.perlの結果の出力ファイル。
  • -l 1 は、構造化の補正に用いる主成分の数で、デフォルトは10。
  • -o example.chisq は、出力される\chi^2検定の結果ファイル。結果は一列目がArmitageの\chi^2統計量、二列目がEIGENSTRATの\chi^2統計量となっている。
  • データが大きいとき(20億ジェノタイプ以上のとき)は、eigenstrat.big.perlを使うらしい。

最後に次の命令で終わる。

gc.perl example.chisq example.chisq.GC

これは、EIGENSTRAT補正が終わった後のファイルをさらにGC補正をかけてくれるというスグレモノらしい。

結果のファイル

  • example.logに、おおまかな結果のまとめが出力される。
  • なんとか.pdfに、構造化の結果が可視化されて出力される。

実際に行うときには

convertfのときと同様、example.perlファイルをひらき、exampleという文字列をすべてtestなどの自分のファイル名の文字列に変更する。

おまけ appendix

A few explanation about UNIX commands.

move to other directories

cd xxx

move a file to other directory

maybe necessary when using .eigenstratgeno, .snp, or .ind file produced by convertf (in directory /CONVERTF) in directory /EIGENSTRAT.

mv xxx.eigenstratgeno ../EIGENSTRAT/xxx.eigenstratgeno

It is also possible to remain the file in /CONVERTF and make the copy of the file in /EIGENSTRAT

cp xxx.eigenstratgeno ../EIGENSTRAT/xxx.eigenstratgeno
ツールボックス

下から選んでください:

新しいページを作成する
ヘルプ / FAQ もご覧ください。