eigenstratは、principal components analysis(主成分解析)の手法をもちいて集団の構造化を解析し、関連検定の結果を補正するソフトウェア。
convertfプログラムでファイル形式をEIGENSTRAT形式とし、smartpca.perlプログラムで構造化の解析を図示するところまで解説する。
おおまかにはbinとEIGENSTRATとCONVERTFのディレクトリが重要。作業はCONVERTFまたはEIGENSTRATのディレクトリ内で行う。まずはCONVERTFディレクトリに移動しておく。
cd CONVERTF
入力ファイル形式は5つが許されている。linkage関連のソフトの標準的な入力フォーマットであり、Haploviewの入力にも使用可能なPED形式がやりやすいと思われる。convertfプログラムを用いて、PED形式のファイルをEIGENSTRAT形式に変換する。
PED形式は三つのファイルから構成される。遺伝子型データは.pedファイル、SNP座位情報(chromosome locationなど)は.pedsnpファイル、サンプル個体データは.pedindファイルに入力する。
これらのファイルの具体的な例を見たいときには、CONVERTFディレクトリ内にあるexample.ped、example.pedsnp、example.pedindファイルをみればよい。すなわち
more example.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 |
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 |
.pedファイルの最初の6-7列目が入力されたもの。
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.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
そのあと次の命令が実行されるが、これは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
../CONVERTF/ind2pheno.perl
最後に次の命令で終わる。
gc.perl example.chisq example.chisq.GC
これは、EIGENSTRAT補正が終わった後のファイルをさらにGC補正をかけてくれるというスグレモノらしい。
convertfのときと同様、example.perlファイルをひらき、exampleという文字列をすべてtestなどの自分のファイル名の文字列に変更する。
A few explanation about UNIX commands.
cd xxx
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