2015-03-02

SNPでIBD〜MCMC_IBDfinderを使う

IBD SNP

これがペイパー
個人のSNPフェーズ不明データから、個人ペアごとに処理してIBD範囲を推定する方法は多いが(Relate,PLINK,BEAGLE,GERMLINEとか(詳細は知らないのですが、MCMC_IBDfinderのペイパーの記載によれば、ということで))、MCMC_IBDfinderは３人以上いたら、いたで、それを統合的に処理してIBD範囲を教えてくれるらしい
近親婚にも対応しているらしい
アプリのダウンロード先はこちら
- ダウンロードして、解凍して、makeする
  - Windowsでは、こちらなどのようにmakeをとってきて(コンパイル用のg++があって)も、エラーを吐いたので、Macで実行。Macでなら、ダウンロード→クリックして解凍→srcディレクトリまで下りて行って、

make

- - すれば実行可能
入力ファイルは３種類
- ジェノタイプファイル
  - タブ区切り
  - 人数分の行、ローカス分の列
  - 値は0,1,2、これがSNPジェノタイプ
- マーカー情報ファイル
  - タブ区切り
  - ２行
    - 第１行は、マーカーのセンチモルガン単位の位置情報
    - 第２行は、マーカーのアレル頻度(帰属集団のそれを取ってくる)
- 初期IBDファイル
  - タブ区切り
  - 人数分の行、ローカスの数ｘ２の列
  - ある個人の行では、第１SNPの第１染色体、第１SNPの第２染色体、第２SNPの第１染色体、第２SNPの第に染色体、…、最後のSNPの第２染色体、という並び方
  - 推定結果のときに、このサイズで整数値の列が帰ってくる。その整数が同じ、ということは、IBDな関係、ということになる
  - 全部１、全部０とかのファイルでよい
  - 例
    - 例えば、以下の場合、３行、１２列だから、３人、６SAPのIBD関連情報となる
    - 第１列は、３人とも１。これは、第１SNPの第１染色体は、３人とも、「１という同じ祖に紐づく、IBD関係」という意味
    - 第２列は、上２人が２、下１人は0。第１SNPの第２染色体は、上二人はIBDだけれど、下一人は関係ない、ということ

1	2	1	2	1	2	1	2	1	2	1	2
1	2	1	2	1	2	1	2	1	2	1	2
1	0	1	0	1	0	1	0	1	0	1	0

実行。ダウンロードするとついてくるサンプルデータファイルで実行

$ ../src/MCMC_IBDfinder -g examplegenotypefile.txt -p exampleposmaffile.txt -z exampleinitzfile0.txt

- このままだと、画面にドバっと出るので

$ ../src/MCMC_IBDfinder -g examplegenotypefile.txt -p exampleposmaffile.txt -z exampleinitzfile0.txt > testout.txt

- とかすれば、以下のようなtestout.txtファイルができる
  - 前半に実行条件など
  - 後半に、初期IBDファイルと同じ構成での出力が、何回も繰り返されている
  - 繰り返しは、人数３の場合は、人数分の行とその前に、繰り返し処理での推定状態を表す行
  - 推定状態を表す行は">"で始まり、６つの値"the sample/iteration number,

the sampled value of , the sampled value of , the logarithm of the transition
probability of the sample, the logarithm of the emission probability of the sample and the time since the sampling started measured in seconds"

- - 結局、最後の人数分の行が推定結果

You are running MCMC_IBDfinder version 1.0 with the following parameter values:

Data parameters
-g	Name of file with genotypes		examplegenotypefile.txt
-p	Name of file with pos maf info		exampleposmaffile.txt

Initialisation parameters
-z	Initial value of z (filename)		exampleinitzfile0.txt
-l	Initial value of lambda			0.200000
-r	Initial value of rho 			0.200000
-s	Seed for random numbers			0

Run parameters
-e	Genotyping error rate on allele level	0.010000
-f	Index of first SNP			0
-k	Max number of IBD sets			2
-ni	Number of iterations			10
-ns	Number of SNPs to include in MCMC	All

Move parameters
-nl	Number of moves of type l		0
-nr	Number of moves of type r		0
-nlw	Number of moves of type lw		1
-nrw	Number of moves of type rw		1
-nzreg	Number of moves of type zreg		350
-nz	Number of moves of type z		150
-nzbor	Number of moves of type zbor		200
-nzmreg	Number of moves of type zmreg		0
-nzxorr	Number of moves of type zxorr		150
-nznes	Number of moves of type znes		0
-nzsnes	Number of moves of type zsnes		0
-nzcp	Number of moves of type zcp		50
-nzbc	Number of moves of type zbc		100

Output parameters
-o	Output destination (filename/screen)	Screen
-t	Thinning t (print every t iteration)	1

MCMC has been started...

> 0	0.200000	0.200000	-8.627707	-91.637074	0
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	
> 1	0.174258	0.188314	-29.077655	-52.712640	0
z 1	2	1	2	1	2	1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	1	0	1	0	
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	1	0	1	0	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	1	0	1	0	
> 2	0.169135	0.212046	-29.260509	-51.973668	0
z 1	2	1	2	1	2	1	2	1	2	1	2	1	0	1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	1	0	1	0	1	0	
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	1	0	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	1	0	1	0	1	0	
> 3	0.216820	0.239780	-29.032881	-51.196142	0
z 1	2	1	0	1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	1	0	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	1	0	1	0	1	0	1	0	
> 4	0.222808	0.189791	-21.661145	-51.820420	0
z 1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	1	0	1	0	1	0	1	0	1	0	
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	1	0	1	0	1	0	1	0	
> 5	0.245059	0.149260	-21.379536	-52.611477	0
z 1	0	1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	1	0	1	0	
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	1	0	1	0	1	0	1	0	1	0	1	0	
> 6	0.264767	0.181247	-21.564392	-50.431730	0
z 1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	0	
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	1	0	1	0	1	0	1	0	
> 7	0.301509	0.137865	-21.556624	-52.828172	0
z 1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 1	0	1	0	1	0	1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	1	0	1	0	1	0	1	0	1	0	
> 8	0.329145	0.181227	-18.466672	-52.044786	0
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 1	0	1	0	1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	1	0	1	0	1	0	
> 9	0.340463	0.191227	-21.948410	-50.831861	0
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 1	0	1	0	1	0	1	0	1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	2	1	2	1	2	1	2	
> 10	0.312030	0.154304	-21.416400	-47.554670	0
z 1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 1	0	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	
z 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	1	2	1	2	1	2	1	2