S.A.G.E.で分離比解析

  • S.A.G.E.のsegregを使って分離比解析に初めて挑戦するための解説文書
  • 講義90分で1-2回分相当
  • アプリケーションはこちら(S.A.G.E)
  • 分離比解析ができる
  • 分厚いマニュアルもあるが、たくさんのことができるので分離比解析に関する章となるとそれほど親切ではない。
  • パラメタファイルのサンプルが簡単に手に入らないのが、初めて使うときにはつらい
  • 下記コードでできる家系例データファイルout.txtやout2.txtと、以下に張り付けるパラメタファイルとを使って、次のコマンドをコマンドラインで打つと、動きます。
  • コマンド
segreg my.par out.txt
segreg my.par out2.txt
  • パラメタファイル
    • その1
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 = "SEGREG2"
{
title = "y"
trait = "affection"
}
    • その2
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 = "SEGREG5"
{
title = "d2"
trait = "affection"
class = "MLM"
type_suscept
{
option = "three"
}
transmission
{
option = "general"
}
}
    • その3
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 = "affection2", binary, affected = "1", unaffected = "0", missing = "."
covariate = "age", continuous, missing = "."

}

segreg, out = "SEGREG2"
{
title = "age"
trait = "affection2", type = age_onset
class = "FPMM"
fpmm
{
loci = "3"
freq = "0.5"
onset
{
type_dependent = "A"
multi_dependent = "N"
age_onset = "age"
age_exam = "age"
}
}
}
  • 詳しくは以下の通り
  • Rmdファイルです。html化、epub化できます(やり方はこちら)
  • html化、epub化が面倒くさければ、kindleで2米ドルでも(こちら)

#家系集積データの扱い

# 本文書の構成


* 遺伝性フェノタイプを持つ家系データのシミュレーション
  * 遺伝性のある2値型フェノタイプの家系データをシミュレーション作成することを通じて、遺伝因子がフェノタイプの家族集積に及ぼす影響を具体性を持って理解する
* 分離比解析の実施
  * フェノタイプの家系データを見て、そこに遺伝性があると判断するための方法としての分離比解析をフリーアプリケーションを使って実行することで、分離比解析とその中で使われている諸概念を理解する

# 遺伝性フェノタイプを持つ家系データのシミュレーション

## 家系データの基本

家系データは、親子関係の記録である。
家系データが正当であるかどうかは(少なくとも現在の生物学的には)、親は2人いて、片方が男で片方が女であることである。
従って、家系データの基本は、

* 本人ID
* その父親のID
* その母親のID
* その性別

の4項目からなる。

## 遺伝性フェノタイプデータ

遺伝性フェノタイプデータを扱うときは、上記の基本4項目に加えて、着目形質が必要である。本文書では2値型の形質を扱うので、フェノタイプは0/1の2通りである。

## その他の項目

上記の4+1=5項目が最小限の項目である。
それ以外の項目は、遺伝性を考える上で共変量として扱いうる。

具体的には

* 年齢(データ取得時の年齢と発病時年齢)
* 遺伝要因
* 環境要因

が挙げられる。分離比解析では、発病に年齢が影響すればそれを、既知の遺伝要因の関与がわかっていればそれを、環境要因の関与が疑われればそれを、共変量として考慮した上で実施できる。そうすることで、考慮した共変量の影響を除いた上での遺伝性を評価することができる。

## 家系データのシミュレーション作成

Rを用いて、遺伝性のあるフェノタイプの家系データを作成する。
* 遺伝要因・フェノタイプと無関係に家系データ(親子関係と性別)を作成し
* その上で背景としての遺伝子座位を想定し、遺伝因子以外の要因(環境要因)によるばらつきを入れて(遺伝要因と環境要因との相対的な大きさを遺伝率として与える)、フェノタイプを確定する。
* 出来上がったデータをテキストファイルで出力し、S.A.G.E.で読み込めるようにする。

### 家系データの作成


家系図を以下の要領で作成する。

* 手順10世代にn.init人を作成する
* 手順2 初期n.init人の性別はランダムに与える
* 手順3 次いで、指定人数 n.iter人を次の要領で付け加える
* 手順4 すでに登録された個人から一人xを選ぶ
  * 手順4-1 xが誰の親でもないときには、新たな配偶者を同世代に1人作成し、そのペアの間に、平均子供数mean.kidでポアソン分布に従う子供を作る
  * 手順4-2 すでに子供が居る場合で、親が居ない場合には、両親を1世代上に作成する
  * 手順4-3 すでに子供がいて親もいるときは何もしない


* 引数
  * 世代0の初期人数n.init
  * 家系拡大ステップ数n.iter
  * 平均子供数mean.kid
* 出力
  * 5列の行列
  * 作成行列は1行に1人
  * 第1列は個人ID
  * 第2列は父親ID(0は情報なし)
  * 第3列は母親ID
  * 第4列は性別(1:male,-1:female)
  * 第5列は世代番号

```{r}
my.pedigree <- function(n.init=3,n.iter=20,mean.kid=2){
  trio <- matrix(0,n.init,5)
	m <- length(trio[1,])
	trio[1:n.init,1] <- 1:n.init
	trio[1:n.init,4] <- sample(c(-1,1),n.init,replace=TRUE)
	trio[,5] <- 0
	cnt <- n.init
	# 「子がいる場合」の「子の数の平均」がmean.kidなので、ポアソン分布に渡すパラメタはmean.kidより小さい
	tmp.f <- function(m,k){
		(m/(1-exp(-m))-k)^2
	}
	xmin <- optimize(tmp.f, c(0, mean.kid), tol = 0.0001, k = mean.kid)
	mean.kid. <- xmin[[1]]
	for(i in 1:n.iter){
		s <- sample(trio[,1],1)
		if(length(which(trio[,2:3]==s))==0){
			# 子がいない(相手がいない)
			trio <- rbind(trio,c(cnt+1,rep(0,m-1)))
			dp <- dpois(1:100,mean.kid.)
			dp <- dp/sum(dp)
			num.kids <- sample(1:100,1,prob=dp)
			tmp <- matrix(0,num.kids,m)
			tmp[,1] <- cnt+1+(1:num.kids)
			if(trio[s,4]==1){
				tmp[,2] <- s
				tmp[,3] <- cnt+1
			}else{
				tmp[,2] <- cnt+1
				tmp[,3] <- s
			}
			tmp[,4] <- sample(c(-1,1),num.kids,replace=TRUE)
			tmp[,5] <- trio[s,5] + 1
			trio <- rbind(trio,tmp)
			trio[cnt+1,4] <- -trio[s,4]
			trio[cnt+1,5] <- trio[s,5]
			cnt <- cnt+1+num.kids
			
			
		}else if(trio[s,2] == 0 & trio[s,3] == 0){
			# 親登録がない
			trio <- rbind(trio,c(cnt+1,rep(0,m-1)))
			trio <- rbind(trio,c(cnt+2,rep(0,m-1)))
			trio[cnt+1,4] <- 1
			trio[cnt+2,4] <- -1
			trio[c(cnt+1,cnt+2),5] <- trio[s,5]-1
			trio[s,2] <- cnt+1
			trio[s,3] <- cnt+2
			cnt <- cnt+2
		}
	}
	trio
}
```

### フェノタイプの割り付け

1個以上の遺伝子座位に3ジェノタイプを想定し、各ジェノタイプに重み(遺伝的リスク)をつける。
その上で、3ジェノタイプの一般集団における頻度を想定すると、上記の関数で作成した家系図において、親が登録されていない個人のジェノタイプは集団のジェノタイプ確率に応じて決まり、
家系図によって親が確定している個人のジェノタイプはメンデル遺伝によって確率的に定まる。

ジェノタイプが定まれば、各座位の重みが決まるので、全座位の重みを単純な和とすれば(ここではそうする)、遺伝的リスクの総和が決まる。

遺伝要因の分散と環境要因の分散との割合を遺伝率として想定することで、作成した家系における遺伝的リスクの分散から環境要因の分散を逆算し、その分散を持つ正規分布を環境リスクとして個人のリスクとして加える。

その上で、個人の総リスクに応じて、指定したリスク以上のリスクを持つ個人を有病者とする。

#### ジェノタイプ割り付け

これを実行するのに、まず、個人のジェノタイプ割り付けを行う

* 引数
  * my.pedigree()関数の出力である5列の行列
  * 一般集団のジェノタイプ確率を表す総和が1の長さ3のベクトル
* 出力
  * 長さが人数分の、ジェノタイプとして値が{0,1,2}のベクトル。{0,1,2}がリスクアレルの本数に相当する。
```{r}
make.genotype <- function(trio,prob){
  ret <- rep(0,length(trio[,1]))
	# 家系図の個人を世代の降順に
	ord <- order(trio[,5])
	# 親がいなければ、集団の確率で
	# 親がいればメンデルの法則で
	for(i in 1:length(ord)){
		if(trio[ord[i],2]==0 & trio[ord[i],3]==0){
			ret[ord[i]] <- sample(c(0,1,2),1,prob=prob)
		}else{
			parents <- trio[ord[i],2:3]
			pat <- ret[trio[ord[i],2]]/2
			mat <- ret[trio[ord[i],3]]/2
			if(pat==0.5){
				pat <- sample(c(0,1),1)
			}
			if(mat==0.5){
				mat <- sample(c(0,1),1)
			}
			ret[ord[i]] <- pat+mat
		}
	}
	ret
}
```

#### 複数のジェノタイプを作ってフェノタイプ割り付け
個々の遺伝子座位のジェノタイプ頻度ベクトルとリスクベクトルを、遺伝子座位の数だけ併せて行列とし、その上で遺伝率と有病者割合を与えて、フェノタイプを確定する。
個人のジェノタイプ割り付けには、上記のmake.genotype()関数を用いる。

* 引数
  * my.pedigree()関数の出力である5列の行列
  * 一般集団のジェノタイプ確率を表す総和が1の長さ3のベクトルを1行とし、座位数を行数とする行列
  * ジェノタイプのリスクを表す長さ3のベクトルを1行とし、座位数を行数とする行列
  * 遺伝率 h (0~1)
  * リスク閾値 prev (0~1)
* 出力
  * リスト
    * 家系図を表す行列
    * ジェノタイプの行列
    * 個人の各座位のリスクを表す行列
    * 個人の遺伝的リスクと環境リスクの総和
    * 0/1のフェノタイプ
```{r}
make.affected.pedigree <- function(trio,p.g,r.g,h,prev){
  # n.lociのジェノタイプを作る
	genotype <- matrix(0,length(trio[,1]),n.loci)
	for(i in 1:n.loci){
		genotype[,i] <- make.genotype(trio,p.g[i,])
	}