「eigenstrat」の編集履歴（バックアップ）一覧はこちら
「eigenstrat」(2007/10/10 (水) 23:32:45) の最新版変更点
追加された行は緑色になります。
削除された行は赤色になります。
eigenstratは、principal components analysis（主成分解析）の手法をもちいて集団の構造化を解析し、関連検定の結果を補正するソフトウェア。
convertfプログラムでファイル形式をEIGENSTRAT形式とし、smartpca.perlプログラムで構造化の解析を図示するところまで解説する。
=convertf=
入力ファイル形式は5つが許されている。linkage関連のソフトの標準的な入力フォーマットであり、Haploviewの入力にも使用可能なPED形式がやりやすいと思われる。convertfプログラムを用いて、PED形式のファイルをEIGENSTRAT形式に変換する。
==PED形式==
PED形式は三つのファイルから構成される。遺伝子型データは.pedファイル、SNP座位情報（chromosome locationなど）は.pedsnpファイル、サンプル個体データは.pedindファイルに入力する。
===.pedファイル===
ダウンロードされたファイルのCONVERTFファイルに入っているexample.pedファイルは次のようになっている。一行目はここでの説明用の行で、ファイルに含まれている情報ではない。
{| style="text-align:center;"
|-
! 家族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ファイルは次のようである。例によって一行目はファイルに含まれていない解説用の行。
{| style="text-align:center;"
|-
! 染色体 !! 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のあとの引数が命令のファイルとなる。このファイルの中身は次のようなもの
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
よくわからないが、実際に使うときはこれをいじくればいいのではないかと思われる。
これで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 は、出力される<math>\chi^2</math>検定の結果ファイル。結果は一列目がArmitageの<math>\chi^2</math>統計量、二列目がEIGENSTRATの<math>\chi^2</math>統計量となっている。
*データが大きいとき（20億ジェノタイプ以上のとき）は、eigenstrat.big.perlを使うらしい。
最後に次の命令で終わる。
gc.perl example.chisq example.chisq.GC
これは、EIGENSTRAT補正が終わった後のファイルをさらにGC補正をかけてくれるというスグレモノらしい。
==結果のファイル==
*example.logに、おおまかな結果のまとめが出力される。
*e
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ファイルは次のようになっている。一行目はここでの説明用の行で、ファイルに含まれている情報ではない。
{| style="text-align:center;"
|-
! 家族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ファイルは次のようである。例によって一行目はファイルに含まれていない解説用の行。
{| style="text-align:center;"
|-
! 染色体 !! 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 は、出力される<math>\chi^2</math>検定の結果ファイル。結果は一列目がArmitageの<math>\chi^2</math>統計量、二列目がEIGENSTRATの<math>\chi^2</math>統計量となっている。
*データが大きいとき（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
