2005-09-14

Ｒのhaplo.statsの使い方（６）haplo.glmの実行

関連解析

必須引数
- データ,mydata(データ型はmodel.matrix:後述)
- 式(たとえば、y~x1+x2; mydataのカラムにy,x1,x2があり、x1とx2とを独立変数、yを従属変数として解析する場合。ただし、今回のジェノタイプデータの場合には、後述のような方法により、xiに複数の多型を一括して指定する(ことが可能))
- モデルファミリ(回帰にあたって使用する仮定のセット。デフォルトはgaussian。データのばらつきが正規分布に従うとの仮定を含む)
- ジェノタイプデータの欠測値の処理方法(一般線形回帰では欠測値は扱わないが、ＥＭによるハプロタイプ推定においては扱えるので、ジェノタイプデータについては、指定により扱えるようにできる)
引数の準備
- 実例に即して進めよう
- mydataの作成
  - 個人別のデータはハプロタイプを推定する一連のジェノタイプデータとそれ以外を別個に準備する
  - ハプロタイプ推定用ジェノタイプデータをテキストファイルから作る
    - まず、data.frameに取り込むryamadaの遺伝学・遺伝統計学メモ - Ｒのhaplo.statsの使い方（２）ハプロタイプ推定
    - ついで、haplo.gmlがハンドリングできるように、NA値の指定とオブジェクト型を"model.matrix"に変換する



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

geno<-setupGeno(geno,miss.val=c(9))

#不明アレルとして9を取ることに決定しているので、miss.valとして9を与えている

#ついでにラベル用の情報も取り込む

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

ハプロタイプと関係ない個人別のデータを取り込む
- 個人の並び方、人数はgenoの作成時と完全に一致させること
- 例として、第１列には、{0,1,NA}、第２列には、{Female,Male,NA}、第３列には、正の有理数の量的形質データが入った３列のデータファイルを使用している。ラベル行をもっていない。
- 欠測値{NA}については、こちらは、何の値も入っていないものとする("\t\t"のように\tと\tとで囲まれることになる)
- 行数が個人の数、列数が因子の数のタブ区切りファイル"data2.txt"から取り込むとすると



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

mydataの作成
- geno,data2という２つのデータフレームを連結して作成する



mydata<-data.frame(geno=geno,data2=data2)

式の決定
- 従属変数としてdata2$V1(data2.txtの第１列の形質)を、data2$V2,data2$V3の２形質と、genoから求めるハプロタイプとを独立変数とする



data2$V1~data2$V2+data2$V3+geno

実行
- 引数na.action="na.geno.keep"で、ハプロタイプ推定にあたっては、不明ジェノタイプコールを使うように指示している。一般線形回帰の処理では、不明データサンプルは除外されている
- 引数allele.lev=attributes(geno)$unique.allelesでは、gml関数で多型アレルをlevelとして渡している。ただし、１多型２列なので、２列に登場するアレルを併せてレベル作成の単位としている(SNPの場合には、アレルは0か1であり、あまり問題とならないが、アレル数が多い多型の場合には、さまざまなアレルが登場するためにこのような処理をしており、この処理を担当しているのは、setupGeno関数である。data2の変数の場合には、それぞれの列ごとに取りうる値がlevelとして渡っている。



> fit.gaus<-haplo.glm(data2$V1~data2$V2+data2$V3+geno,family=gaussian,na.action="na.geno.keep",allele.lev=attributes(geno)$unique.alleles,data=mydata,locus.label=label1$V1)

結果の表示は以下の通り
- data2$V2は{Female,Male}のDichotomous型データなので、それぞれを因子として結果が返っている。data2$V3は量的変数なので、１変数として結果が返っている。
- ハプロタイプはそれぞれ、geno1,geno2,...という呼び方で因子として表示され、それぞれのgenoNが現すハプロタイプについては、末尾にその構成と頻度が示されている
- ハプロタイプの寄与は、最も頻度の高いハプロタイプ"haplo.base"を１としてそれに対する値として出されている
- アレル頻度が小さいハプロタイプはgeno.rareとしてプールして処理するようになっているが、プールしても頻度が一定以上ある場合にのみ解析対象となる。そのデフォルト値は0.001である。これに満たない場合、低頻度ハプロタイプは基準ハプロタイプとプールされる

> fit.gaus

Call:
haplo.glm(formula = data2$V1 ~ data2$V2 +
data2$V3 + geno, family = gaussian, data = mydata,
na.action = "na.geno.keep", locus.label = label1$V1,
allele.lev = attributes(geno)$unique.alleles)

Coefficients:
coef se t.stat pval
(Intercept) 2.00e+00 0.06251 32.00971 0.000
data2$V2Female -1.00e+00 0.06142 -16.28090 0.000
data2$V2Male -1.00e+00 0.06394 -15.63982 0.000
data2$V3 -2.13e-03 0.01339 -0.15939 0.874
geno.1 -1.00e+00 0.01205 -82.96941 0.000
geno.4 3.99e-05 0.01740 0.00229 0.998
geno.5 7.32e-05 0.00604 0.01212 0.990
geno.rare -3.33e-01 0.03384 -9.84123 0.000

Haplotypes:
X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 hap.freq
geno.1 0 0 0 0 0 0 0 0 0 0 0.07109
geno.4 1 0 0 1 1 0 0 0 1 0 0.02868
geno.5 1 0 0 1 1 0 1 0 0 0 0.36730
geno.rare * * * * * * * * * * 0.00711
haplo.base 1 0 0 1 1 0 0 0 0 1 0.52582
>