2005-09-13

Ｒのhaplo.statsの使い方（２）ハプロタイプ推定

関連解析

推定アルゴリズム
- ＥＭベース
  - ＥＭとの相違は、
    - 推定の過程で、設定閾値に満たない頻度のハプロタイプを除去する
    - 推定の初めから、全多型の作るハプロタイプ空間を相手にせずに、解析多型数を増やしながら進める(SNPHAP by David Clayton http://www-gene.cimr.cam.ac.uk/clayton/software/snphap.txt:SNPHAPと同様)
使用法
- 入力は２方法
  - (1)個人別のdiploidジェノタイプデータと、多型名データを入力する
  - (2)同一ジェノタイプの個人について、人数をカウントし、diploidジェノタイプデータとその人数、および、多型名データをを入力する
- さまざまな入力形式に対応できるように作られているが、それが、かえって、使用法のオフィシャル説明文をわかりにくくしているので、思い切って、入力フォーマットを次のように固定したテキストファイルで所有しているものとする
- 入力方法(1)と(2)とでは、圧倒的に(2)の方が早い
- 入力方法(1)の場合
  - ジェノタイプファイル
    - 行数は人数、列数は多型の個数x2
    - １多型２列
    - タブ区切り
    - ２アレルは、0か1、不明は9とする
    - ホモの場合は"0\t0(\t)" もしくは "1\t1(\t)" ヘテロの場合は"0\t1\t" 不明の場合は"9\t9(\t)"
    - サンプルファイル(４人、４多型)"genotype.txt"



1	1	0	0	1	1	9	9

0	1	0	1	0	0	1	1

1	1	0	1	0	0	0	0

9	9	0	1	1	1	0	1

多型名データファイル
- Haploviewの.infoファイルと同形式とする
- 行数は多型数、列数は２
- 第１カラムは多型名、第２カラムは多型位置
- タブ区切り
- サンプルファイル(４多型)"snp.info"



SNP34	47063263

SNP56	47063331

SNP62	47063675

SNP78	47067844

実行
- コマンドは以下の５つ



library(haplo.stats)

geno<-read.table("genotype.txt",FALSE,"\t")

label1<-read.table("snp.info",FALSE,"\t")

save.em.keep<-haplo.em(geno,locus.label=label1$V1,miss.val=9)

print.haplo.em(save.em.keep)

第１コマンドでhaplo.statsパッケージの読み込み
第２コマンドでgenotypeデータをgenoオブジェクトに取り込み
第３コマンドで、多型名データをlabel1オブジェクトに取り込み
第４コマンドで推定実行(ジェノタイプデータをgenoで引渡し、多型名とlabel1オブジェクトの第１列(label1$V1)で渡している。不明値は9なので、miss.val=9と指定している)
第５コマンドで結果の表示そのサンプル(10SNP,251人の結果。最上段はハプロタイプ名(この場合は1,2,...10)、２段目からハプロタイプIDとその多型のアレルと推定頻度。最下段は推定にあたって出た値(詳細、省略)



> print.haplo.em(save.em.keep)

============================================================================================================ 

                                                 Haplotypes                                                  

============================================================================================================ 
   1 2 3 4 5 6 7 8 9 10 hap.freq

1  0 0 0 0 0 0 0 0 0  0  0.06834

2  0 0 0 0 0 1 0 0 0  0  0.00244

3  1 0 0 1 1 0 0 0 0  1  0.51803

4  1 0 0 1 1 0 0 0 1  0  0.02763

5  1 0 0 1 1 0 1 0 0  0  0.37622

6  1 0 0 1 1 0 1 1 0  0  0.00228

7  1 0 1 1 1 0 1 0 0  0  0.00258

8  1 0 1 1 1 0 1 1 0  0  0.00002

9  1 1 0 1 1 0 1 0 0  0  0.00111

10 1 1 0 1 1 0 1 1 0  0  0.00002

11 1 1 1 1 1 0 1 0 0  0  0.00132

12 1 1 1 1 1 0 1 1 0  0  0.00002

============================================================================================================ 

                                                  Details                                                    

============================================================================================================ 
lnlike =  -369.8037 

lr stat for no LD =  312.1477 , df =  1 , p-val =  0

入力方法(2)の場合
- ジェノタイプファイル
  - 行数はジェノタイプパターン数、列数は多型の個数x2
  - １多型２列
  - タブ区切り
  - ２アレルは、0か1、不明は9とする
  - ホモの場合は"0\t0(\t)" もしくは "1\t1(\t)" ヘテロの場合は"0\t1\t" 不明の場合は"9\t9(\t)"
  - サンプルファイル(全８人が４パターンのジェノタイプに分かれたものとする、４多型、)"genotype_w.txt"



1	1	0	0	1	1	9	9

0	1	0	1	0	0	1	1

1	1	0	1	0	0	0	0

9	9	0	1	1	1	0	1

ジェノタイプパターン別人数ファイル"weight.txt"(第1パターンに３人、第2、３パターンに２人、第4パターンに１人の計、3+2+2+1=8人)

多型名データファイルは、(1)と同じ
- Haploviewの.infoファイルと同形式とする
- 行数は多型数、列数は２
- 第１カラムは多型名、第２カラムは多型位置
- タブ区切り
- サンプルファイル(４多型)"snp.info"



SNP34	47063263

SNP56	47063331

SNP62	47063675

SNP78	47067844

実行
- コマンドは以下の６つ



library(haplo.stats)

geno<-read.table("genotype_w.txt",FALSE,"\t")

weight1<-read.table("weight.txt",FALSE,"\t")

label1<-read.table("snp.info",FALSE,"\t")

save.em.keep<-haplo.em(geno,locus.label=label1$V1,weight=weight1,miss.val=9)

print.haplo.em(save.em.keep)

第１コマンドでhaplo.statsパッケージの読み込み
第２コマンドでgenotypeパターン別人数のデータをweightオブジェクトに取り込み
第３コマンドでgenotypeデータをgenoオブジェクトに取り込み
第４コマンドで、多型名データをlabel1オブジェクトに取り込み
第５コマンドで推定実行(ジェノタイプデータをgenoで引渡し、多型名とlabel1オブジェクトの第１列(label1$V1)で渡している。不明値は9なので、miss.val=9と指定している。パターン別人数をweight=weight1として渡している)
第６コマンドで結果の表示そのサンプル(10SNP,41パターンに分かれる251人の結果。最上段はハプロタイプ名(この場合は1,2,...10)、２段目からハプロタイプIDとその多型のアレルと推定頻度。最下段は推定にあたって出た値(詳細、省略)
(1)でのサンプル出力と元データは同一。推定ハプロタイプ頻度がわずかに異なることと、推定にあたっての尤度などが異なることに注意。



> print.haplo.em(save.em.keep)

============================================================================================================ 

                                                 Haplotypes                                                  

============================================================================================================ 
   1 2 3 4 5 6 7 8 9 10 hap.freq

1  0 0 0 0 0 0 0 0 0  0  0.06834

2  0 0 0 0 0 1 0 0 0  0  0.00244

3  1 0 0 1 1 0 0 0 0  1  0.51803

4  1 0 0 1 1 0 0 0 1  0  0.02763

5  1 0 0 1 1 0 1 0 0  0  0.37622

6  1 0 0 1 1 0 1 1 0  0  0.00230

7  1 0 1 1 1 0 1 0 0  0  0.00258

8  1 0 1 1 1 0 1 1 0  0  0.00001

9  1 1 0 1 1 0 1 0 0  0  0.00103

10 1 1 0 1 1 0 1 1 0  0  0.00001

11 1 1 1 1 1 0 1 0 0  0  0.00140

12 1 1 1 1 1 0 1 1 0  0  0.00001

============================================================================================================ 

                                                  Details                                                    

============================================================================================================ 
lnlike =  -369.8035 

lr stat for no LD =  -491.8254 , df =  1 , p-val =  1