分離比解析をS.A.G.E.で - ryamadaの遺伝学・遺伝統計学メモ

分離比解析に関連して家系データの作り方や家系図描図について書いた(昨日など
そこで作った家系データをS.A.G.E.のsegregという分離比解析アプリにつなごうと思う
分離比解析のための基礎知識
- 分離比解析でやること
  - 家系で0/1のフェノタイプ情報を集める(今は0/1フェノタイプのみを対象とする)
  - 主要な遺伝子座位があるかどうかを判定する
  - そのような遺伝子座位の遺伝形式と「タイプ」別リスクの大きさを推定する
- 解析・推定の主な仕組み
  - 家系図構造の情報を使う(当然…)
  - 発病するかしないかに関してロジスティック回帰的な式のあてはめを行う
    - 共変量(生活習慣・年齢など)を取り入れることができる
  - 普通のロジスティック回帰は一般化線形回帰の一種なので、「線形回帰」だが、家系図構造があるので、そうはいかないので、最尤推定値探索を行う
  - 遺伝子座位のジェノタイプ情報を使う代わりに、個人が持っていて、親が子に伝達する「タイプ」というものを想定する
「タイプ」と遺伝形式というのが、ちょっとわかりにくいので説明する
- 「タイプ」は３種類ある。３種類あると仮定する
- ３種類の「タイプ」は"AA","AB","BB"のように書くのでSNPのジェノタイプのようだが、もう少しチュ小的なもの
- 「タイプ」は２つの顔を持つ
  - １つは0/1の疾患フェノタイプの疾患感受性を決める。発病確率であったり、年齢の関数としての発病率であったりする
  - もう１つは、疾患フェノタイプを子にどのように伝達する要因であるか、という顔である
    - たとえば「優性遺伝形式」の発病者は原因アレルのホモかヘテロであり、ホモの場合は子のフェノタイプは病気が１、非病気が０となる。「優性遺伝形式の原因アレルのホモ」という「タイプ」は(1,0)の確率を子に起こす、メイティングの相手によらずそのような確率で子のフェノタイプが決まる、という性質を持つ。このように、個人は「タイプ」を持っており、それはメイティング相手の「タイプ」が決まると、子のフェノタイプがどちらにどれくらいの確率でなるかを決めるものである、というように考える
    - ただし、ここで、子のフェノタイプの確率が親の「タイプ」で確定してしまうことにすると、ペネトランスとか、「タイプ別の発病確率」を推定対象とする必要がなくなる。であるから、ここでいう「タイプ」が、「子に伝わる影響」を「発病」とは分離して、ピュアなものにすることにしよう、として持ち出した考え方が、"transmission"である。ここで、親が子に伝える要素は連続量で考えてもいけなくはないのだが、(おそらく、離散的な考えをする伝統の影響か？）、親の「タイプ」が、子の「タイプ」に影響を与える、というときに、次の様に考えることにした
    - "AA"というタイプが"A"を伝達する確率、"AB"というタイプが"A"を伝達する確率、"BB"というタイプが"A"を伝達する確率、というのを変数にする。A,Bを実在としてのアレルとすると、"AA"が"A"を伝える確率は1であり、"AB"のそれは0.5であり"BB"のそれは0なのだが、「抽象的な概念」として導入したから、こだわらなくてもよいことにして、推定対象にしましょうということにした
推定するのは結局何か
- タイプ"AA","AB","BB"の疾患感受性・発病率(共変量と組み合わせて個人の発病率がモデル式で定められる。そのモデル式を構成する変数が推定対象。ロジスティック回帰の１成分としてタイプが用いられるので $\beta_{AA},\beta_{AB},\beta_{BB}$ )
- "AA",AB","BB"のタイプが子に"A"を伝達する確率( $\tau_{AA},\tau_{AB},\tau_{BB}$ )
- 家系図上の個人の「タイプ」は親がいれば親の「タイプ」と $\tau_{..}$ とで決まるが、家系図の中には親を持たない個人もいて、その個人の「タイプ」の割合は「集団」における「タイプ」の割合の影響を受けるので、「集団」における「タイプ」割合も推定対象
「推定するだけではない」が、では「何をするか」
- いくつかのモデルに合わせて、上記推定対象の最尤推定を行った後、どのモデルを採用するのが妥当であるかを、モデルの自由度とモデルに得られた最大尤度とから、採択するべきモデルを決める(AICで)
- モデルは大きく分けて次の二つ
  - 遺伝因子がないとするモデル
  - 遺伝因子があるとするモデル
- 遺伝因子があるとするモデルは次の２セットのパラメタに設定する制約の有無・制約の具合で複数に分けられる
  - $\beta_{AA},\beta_{AB},\beta_{BB}$ の大小関係の制約
  - $\tau_{AA},\tau_{AB},\tau_{BB}$ の大小関係の制約
S.A.G.E.アプリの導入と使い方
- フリーソフトS.A.G.E.はこちらからダウンロード
- 基本的には、指示に従ってインストールするだけ
- 入力ファイルは２つだけ
  - (1)家系情報(親子関係と性別)と個人のフェノタイプ情報、その他の共変量情報のファイル
  - (2)パラメタファイル
- (1)家系情報＋個人情報ファイルは、後述のように、「Rで家系データ作成」からテキストファイルに書きだすので、それを使えばよい。適切なデリミタで区切られたテキストファイルであって、ヘッダー行にカラムの内容を指定してあればそれでよい
  - 少し注意するのは、「欠測値」が何であるかの指定が必要だること
- (2)パラメタファイルは、２つの部分からなる
  - (i) 家系情報＋個人情報ファイルの作り(なんという名前のカラムに何の情報があるか、欠測値は何か、など)を説明する部分と
  - (ii) S.A.G.E.の中の特定の解析(今回はsegreg)の実行条件を指定する部分と
- 実は、この「パラメタファイル」が厄介で、なかなか「さらさら」と書けない。GUIを使うと「あーら、簡単にできますね！」というのを目指しているらしいのだが、それもなかなか思うようにいかない
- 従って、「これがパラメタファイルのデフォルト」というのをテキストファイルで作っておいて実行するのがよさそう
- 慣れてきたら、手作業で書き換えるのもGUIを使って作成するのも見通しがよくなって、大丈夫なはず

pedigree
{
delimiters = "\t"
delimiter_mode = "single"

individual_missing_value = "0"
sex_code, male = "1", female = "-1", missing = ".", trait

pedigree_id = "family"
individual_id = "ind"
parent_id = "father"
parent_id = "mother"
sex_field = "sex"

trait = "affection", binary, affected = "1", unaffected = "0", missing = "."

}

segreg, out = "SEGREG14"
{
title = "b"
trait = "affection"
class = "MLM"
type_suscept
{
option = "three_add"
suscept = "A*"
suscept = "AA"
suscept = "AB"
suscept = "BB"
}
transmission
{
option = "general"
tau = "AA", val = "1"
tau = "AB", val = "0"
tau = "BB"
tau = "BB", val = "0"
tau = "A*", val = "1"
}
}

Rで作った家系データをS.A.G.E.のsegregに解析させる
- う

### 多家系作成条件
# n.loci個のリスク座位
n.loci <- 10
# それらのアレル頻度とHWEを仮定したジェノタイプ頻度
p.a <- runif(n.loci)*0.1
p.a[1] <- 0.4
p.g <- cbind(p.a^2,2*p.a*(1-p.a),(1-p.a)^2)

# 各座位のジェノタイプ別リスク
r.a <- c(1000,rep(1,n.loci-1))
r.g <- cbind(r.a*2+1,r.a+1,rep(1,n.loci))

herit <- 0.95

prev <- 0.4
n.init <- 1
n.iter <- 50
mean.kid <-2 

k<- 5

n.family <- 10
#jpeg("test.jpeg",width=10000,height=10000)
#par(mfcol=c(5,5))
F <- list()
for(i in 1:n.family){
	my.ped <- my.pedigree(n.init=n.init,n.iter=n.iter,mean.kid=mean.kid)
	aff.pedigree <- make.affected.pedigree(my.ped,p.g,r.g,herit,prev)
	#plot.affected.pedigree(aff.pedigree)
	el <- el.from.trio(aff.pedigree$pedigree)
	# エッジリストからグラフオブジェクトを作る
	g <- graph.edgelist(el)

	smpl.aff.ped <- sample.aff.sub(g,aff.pedigree,k)
	F[[i]] <- smpl.aff.ped
	#plot.affected.pedigree(smpl.aff.ped,10)
}
#par(mfcol=c(1,1))
#dev.off()

genotype.change <- function(g.mat){
	ret <- g.mat
	ret[which(ret==0)] <- "0/0"
	ret[which(ret==1)] <- "0/1"
	ret[which(ret==2)] <- "1/1"
	ret
}

make.multi.family.table <- function(F){
	family.id <- rep(1,length(F[[1]]$pedigree[,1]))
	pedigree <- F[[1]]$pedigree
	genotype <- genotype.change(F[[1]]$genotype)
	gen.risk.mat <- F[[1]]$gen.risk.mat
	gen.risk <- F[[1]]$gen.risk
	affected <- F[[1]]$affected
	
	if(length(F)>1){
		for(i in 2:length(F)){
			family.id <- c(family.id,rep(i,length(F[[i]]$pedigree[,1])))
			pedigree <- rbind(pedigree,F[[i]]$pedigree)
			genotype <- rbind(genotype,genotype.change(F[[i]]$genotype))
			gen.risk.mat <- rbind(gen.risk.mat,F[[i]]$gen.risk.mat)
			gen.risk <- c(gen.risk,F[[i]]$gen.risk)
			affected <- c(affected,F[[i]]$affected)
		}
	}
	tmp <- cbind(family.id,pedigree,affected,genotype,gen.risk.mat,gen.risk)
	dimnames(tmp)[[2]] <- c("family","ind","father","mother","sex","generation","affection",paste("genotype",1:length(genotype[1,]),sep="_"),paste("locus_risk",1:length(genotype[1,]),sep="_"),"genetic_risk")
	tmp
}

tbl <- make.multi.family.table(F)
#write.table(tbl,"out.txt",sep="\t",quote=FALSE,na=".")
write.matrix(tbl, file = "out2.txt", sep = "\t")

S.A.G.E.のsegregの出力ファイルの読み