2015-03-04

sQTL　GLiMMPS

GLiMMPS sQTL

論文
- 背景モデルの説明PDF
アプリのサイト
概要
- NGSでデータを取ると、エクソンが抜けるか抜けないかについて、「抜けていないことが確実なリード(対象エクソンとその前後どちらかのエクソンとにまたがったリード)本数と、エクソンをスキップしたことが確実なリード(エクソンをまたいで連結した２エクソンにまたがるリード)との本数がわかる
- 前者をIJ、後者をSJとすると、y=IJ/2,n=IJ/2+SJが、「抜けていないトランスクリプト」「抜けていないのと抜けているのを併せたトランスクリプト」の本数になる
- これを使って、「抜けていないトランスクリプトが全体に占める割合」について考えたい
- もっとも単純には、 $\psi=y/n$ という算術割合を使うのもよい。たくさん読めていればそれでもかなりよい
- 問題は、あまり読めていない場合も扱わなくてはいけないこと
- 考え方として三段階
  - (1)y/nを使う。今、興味があるのは、(周辺SNPの)ジェノタイプ $g$ とy/nの関係なので、それを線形回帰する( $y/n = \beta_0 + \beta_1 \times g + \epsilon$
  - (2)リード数が多い人と少ない人とのy/nを平等に扱っているので、気持ち悪い。それに対処するべく、 $\psi$ がジェノタイプの関数 $\log{\frac{\psi}{1-\psi}}=\beta_0 + \beta_1 \times g$ であるとして、このジェノタイプごとの $\psi$ を使って、リード総本数n(これはポアソン分布だったり、それよりばらつきの大きい負の二項分布だったりすると考えるのがふつう)に対して、二項分布で決まるのがyと考える
  - (3)実際のデータは、それほどジェノタイプが同じなら、が似ている、というわけではないので、個人別の影響を加えてのように、iさんのがジェノタイプで決まる値に正規乱数があって、しかも、この正規乱数は、iさん個人で決まるバラツキを持つようにとする、というもの
    - 最後の「個人差」を考慮するので、generalized linear mixed model (GLMM)
まずは、データをシミュレーション作成することで、何を考慮しているかを確認する

# (1) 線形回帰のモデルに相当する
beta0 <- 0.3 # 切片分
beta1 <- 0.1 # ジェノタイプの寄与分
# ジェノタイプを0,1,2で作る
g <- c(rep(0,490),rep(1,420),rep(2,90))
# 総人数
N <- length(g)
# y/nという「割合」を対象にする
# 線形の式に乱雑項を入れる
y.n.lm <- rep(0,N)
for(i in 1:N){
	y.n.lm[i] <- beta0 + beta1*g[i]+rnorm(1,0,0.1)
}

par(mfcol=c(1,3))
hist(y.n.lm)
for(i in 0:2){
	hist(y.n.lm[which(g==i)],col=2+i,add=TRUE)
}

# psiがジェノタイプで決まるとするモデル
# 回帰するときは一般化線形回帰glm
library(gtools) # inv.logit()関数を使うために入れる
# nをポアソン発生させ、そのうえでpsiに応じてyを決める
y <- n <- psi <- rep(0,N)
n <- rpois(N,30)
for(i in 1:N){
	psi[i] <- inv.logit(beta0+beta1*g[i])
	y[i] <- rbinom(1,n[i],psi[i])
}
y.n.glm <- y/n
hist(y.n.glm)
for(i in 0:2){
	hist(y.n.glm[which(g==i)],col=2+i,add=TRUE)
}

# GLiMMPSモデル
# 回帰するときは、GLMMする
y <- n <- psi <- rep(0,N)
n <- rpois(N,30)
# 個人ごとに正規分布に従うSDを取らせる
sigmai <- abs(rnorm(N,0,0.1))
for(i in 1:N){
	tmpu <- rnorm(1,0,sigmai[i])
	psi[i] <- inv.logit(beta0+beta1*g[i]+tmpu^2)
	y[i] <- rbinom(1,n[i],psi[i])
}
y.n.glimmps <- y/n
hist(y.n.glimmps)
for(i in 0:2){
	hist(y.n.glimmps[which(g==i)],col=2+i,add=TRUE)
}
par(mfcol=c(1,1))
mean(y.n.lm)
mean(y.n.glm)
mean(y.n.glimmps)
var(y.n.lm)
var(y.n.glm)
var(y.n.glimmps)

GLiMMPSはRとpythonのパッケージだが、その中で、どうやっているのかをみると回帰の方がわかりやすい
- アプリのサイトからとってくれば、Rscriptsディレクトリにソースがある
- (1)の方法。以下のようにy/n=psi自体をジェノタイプと線形回帰で評価している

lm(psi ~ SNP)

- (2)の方法。個人要素を入れないで、psiをジェノタイプの関数としてinv.logit()で対応付けるときには、以下のように、ジェノタイプでpsi対応値にして(そのあとで観測本数と突き合わせる)

glm(response ~ SNP, family = binomial)

- (3)の方法。個人要素を入れるには、GLMMをする
  - (1|obs)が個人要素を考慮する項。obsは全員に異なる値を与えて使っており、glmer()関数のヘルプ記事にもある通り"Random-effects terms are distinguished by vertical bars ("|") separating expressions for design matrices from grouping factors."

library(lme4) # glmer()関数を使うため
glmer(response ~ SNP +(1|obs), family=binomial))