Multiple alignment

MikuHatsune2015-05-12

Multiple alignment してほしいんだけど、と言われて系統樹まで作った話。
遺伝的距離や系統樹についてはこちらが詳しいし、塩基配列の変化パターンによる重み付けについては昔やった
Multiple alignment はClustal がよく言われるが、ClustalW2を使ってみる。
まずは入力のためのFASTAファイルを作成する。

library(seqinr)
library(ape)
library(phangorn)

data(woodmouse)
w <- as.alignment(woodmouse)
write.fasta(lapply(w$seq, s2c), w$nam, file="woodmouse.fasta")

タブでDNAにして、入力ボックスにコピペしてもいいし、ファイルをアップロードしてもよい。
output option は PHYLIP にしておくと後が楽っぽい。
置換や欠損のペナルティをいくつかにするかは適宜変更する。
このあとで phylip 形式の系統樹用のデータを入手してグラフを作る。Rで作るのが面倒だったらブラウザから系統樹のデータを取ってきて、read.tree をしたら一緒になる。

ph <- read.alignment("aln-phylip.txt", format="phylip")
phyDat <- as.phyDat(ph)
tree0 <- NJ(dist.ml(phyDat))
plot(tree0)

海の人間と陸の人間の遺伝学

MikuHatsune2014-01-27

凪のあすからで、「海と陸の人間との間に生まれた子供は海中では生きていけない」という設定があるのだが、海の人間であるみおりと、陸の人間である至のこどもである美海が、急に海の中で息ができるようになったので遺伝学的に考察する。
家系図の作成にはkinshipを使っていたが、いつのまにかkinship2なるパッケージがきていたので使ってみる。
 
「海と陸の人間との間に生まれた子供は海中では生きていけない」という情報から、海で生きる表現型は劣性なのだろうと思う。ただ、元々人間は海にいて、陸で生きるものも出てきた、という設定らしいので、進化的にどうなの?とは思うがそういうことなんだろう。
で、海の人間と陸の人間のこどもである美海が、いままでは海で息ができなかったが急に息ができるようになったらしいので、海の表現型が発現(浸透率が100%ではなかった)のか、至が海のアレルを持っていてそもそも美海が劣性ホモで、表現型の発現に時間がかかったとか、そんな議論になった。
 
海の人間と陸の人間でGWASしようとか、そうなったらHWEが成立していないとか、海村が14しかこの世界に存在していないらしく、近交系になるので集団構造化が必要とか、遺伝的浮動シミュレーションはどうか、とか、そんな話で盛り上がった。

library(kinship)

id <- c("母", "灯", "みおり", "至", "あかり", "光", "美海", "晃")
Mother <- c(NA, NA, NA, NA, "母", "母", "みおり", "あかり")
Father <- c(NA, NA, NA, NA, "灯", "灯", "至", "至")
Sex <- c(2, 1, 2, 1, 2, 1, 2, 1)
affected <- c(2, 2, 2, 1, 2, 2, 2, 1)
status <- c(0, 0, 1, 0, 0, 0, 0, 0)

ped1 <- pedigree(id, Father, Mother, Sex, affected, status)
pat <- which(id == "美海")
p1 <- plot(ped1, cex=1.5, mar=c(12, 1, 6, 1))
#xy <- locator(1)
text(p1$x[pat], p1$y[pat], "→", cex=3, srt=45, xpd=TRUE, adj=c(1.7, 0.7))
title("凪のあすから 遺伝学", cex.main=3)

text1 <- "凪のあすから16話で,美海が海の中で呼吸できた。「海と陸の人間との間に\n
生まれた子供は海中では生きていけない」という設定から,海中で生きる\n
表現型は劣性形質(a)だと考えられる。しかし,美海が海の表現型であったこと\n
を考慮すると,(1)美海はAaだが,浸透率が100%ではなかったので海の表現型\n
になった,(2)至がAaで,美海はaa,つまり至の両親のどちらかは海の人間\n
だった説になる。紡の祖父の勇が元々海の人間だったことや,晃の表現型の\n
情報があればもっと説明つきそう。伴性遺伝やゲノム刷り込みとか小難しいこと\n
は知らん。"

text(par()$usr[1]-0.07, par()$usr[3]+0.8, text1, xpd=TRUE, pos=4, cex=1.5)

idx <- which(affected == 2 & id != "美海")
text(p1$x[idx], p1$y[idx], "aa", cex=2, pos=3, xpd=TRUE)
text(p1$x[which(id == "至")], p1$y[which(id == "至")], "(1) AA\n(2) Aa", cex=2, pos=3, xpd=TRUE)
text(p1$x[pat], p1$y[pat], "(1) Aa\n(2) aa", cex=2, pos=3, xpd=TRUE)

Multifactor Dimensionality Reduction (MDR)

MikuHatsune2013-09-02

Multifactor Dimensionality Reduction(MDR)という手法がある。意味のある組み合わせを考えたいが、組み合わせ爆発が起こるのでうまくbin化しながら考えるようなイメージだと思った。これを遺伝医学に応用した話。
Am J Reprod Immunol. 2013 Jul 31.
RではMDR, mbmdrなどでできる。元々はBioinformatics. 2003 Feb 12;19(3):376-82.


STEP1
CVするために n-fold CVとして今回は1/10にした。
STEP2
N個のlocusがfactorとしてある。
STEP3
遺伝子座Lについて、L_iL_j を取ってきて、3*3のカテゴリーでcase と controlの人数を考える。
STEP4
リスクが1以上なら high, 1以下なら low, 割り付けが存在しなかったら emptyのラベルに貼り替える。
STEP5
L_iL_jの組み合わせの時L_{i,j}のモデルのときの誤差をひたすら計算しておく。
STEP6
検証データで誤差を確認する。

library(mbmdr)
data(simSNP)
fit <- mbmdr(y=simSNP$Y, data=simSNP[, 3:12], order=2, family=binomial(link=logit), printStep1=TRUE) # print the results of step one

SNP1とSNP2を取ってきて計算した時の結果を示す。SNP1の遺伝子型がAA/Aa/aaのときとSNP2の遺伝子型がBB/Bb/bbのときの9パターンにおいて上のステップの通り計算し、Wald testのp値を返している。

Model:  2 1 
 SNP2 SNP1 cases controls    beta   p.value category
    0    0     0        9 -2.9904 0.0510640        L
    1    0    49       19  1.1286 0.0001103        H
    2    0     0       15 -3.5117 0.0182311        L
    0    1    50       19  1.1554 0.0000727        H
    1    1     0       57 -5.0794 0.0004079        L
    2    1    43       30  0.4396 0.0938642        H
    0    2     0       14 -3.4397 0.0211741        L
    1    2    58       26  1.0056 0.0001219        H
    2    2     0       11 -3.1919 0.0348002        L

permutation testでp値が小さかった順にSNPの組み合わせを返す。

print(fit)
  SNP1 SNP2 NH  betaH     WH        PH NL   betaL     WL        PL     MIN.P
  SNP2 SNP1  4 6.1135 18.285 1.902e-05  5 -6.1135 18.285 1.902e-05 1.902e-05
  SNP8 SNP2  1 0.5195  2.771 9.598e-02  2 -0.7656  7.972 4.750e-03 4.750e-03
  SNP8 SNP7  1 0.6271  3.592 5.806e-02  1 -0.8908  6.901 8.614e-03 8.614e-03
 SNP10 SNP4  2 0.6802  6.182 1.291e-02  1 -0.7414  3.099 7.836e-02 1.291e-02
 SNP10 SNP9  0     NA     NA        NA  1 -1.1196  6.156 1.309e-02 1.309e-02
 SNP10 SNP1  0     NA     NA        NA  1 -0.7777  5.435 1.974e-02 1.974e-02
  SNP6 SNP1  0     NA     NA        NA  1 -0.5302  4.930 2.640e-02 2.640e-02
  SNP8 SNP6  1 0.7379  4.845 2.773e-02  0      NA     NA        NA 2.773e-02
  SNP8 SNP4  1 0.8045  2.940 8.640e-02  1 -0.6804  4.518 3.355e-02 3.355e-02
  SNP6 SNP2  0     NA     NA        NA  1 -0.6804  4.518 3.355e-02 3.355e-02
  SNP7 SNP3  1 1.0338  4.457 3.475e-02  0      NA     NA        NA 3.475e-02
  SNP7 SNP2  0     NA     NA        NA  1 -0.5757  3.609 5.747e-02 5.747e-02
  SNP6 SNP4  1 0.8744  3.540 5.989e-02  0      NA     NA        NA 5.989e-02
  SNP8 SNP1  1 0.5572  3.503 6.126e-02  0      NA     NA        NA 6.126e-02
 SNP10 SNP2  1 0.5390  3.132 7.677e-02  0      NA     NA        NA 7.677e-02
  SNP9 SNP6  1 0.8045  2.940 8.640e-02  1 -0.5849  3.089 7.883e-02 7.883e-02
  SNP4 SNP2  1 0.5214  3.037 8.137e-02  0      NA     NA        NA 8.137e-02
  SNP5 SNP1  1 0.4901  2.867 9.040e-02  0      NA     NA        NA 9.040e-02

論文では、決定木を使って疾患リスクを推定している。

library("rpart")
classification_tree <- rpart(Y ~ . -X, data = simSNP)

#finding the best optimal size of tree
cp <- classification_tree$cptable[which.min(classification_tree$cptable[, "xerror"]), "CP"]# select the complexity parameter associated with minimum error
Regtree_optimal <- prune(classification_tree, cp=cp)
plot(Regtree_optimal, uniform=T, main="pruned regression Tree for trait")
text(Regtree_optimal, use.n=TRUE, xpd=TRUE)

ベイズ的思考な遺伝相談

医師国家試験より
97I8

32歳の女性。遺伝相談のために来院した父が常染色体優性遺伝疾患に罹患している。
この疾患の浸透率(penetrance)は50%。近親婚はなく、本人は発症していない。
家系図を次に示す。第1子が発症する確率はどれか。

常染色体優性遺伝疾患とは、普通ヒトには染色体が父親由来と母親由来の2本あって、そのうち少なくとも1本に遺伝素因があれば病気になる疾患のことで、
浸透率(penetrance)とは、じゃあその遺伝素因があれば100%その疾患になるの?というとそうでもなく、発症しないままのことがあり、その発症確率のこと。
ベイズ的には、「相談者が発症していない」という情報がカギ。ちなみに相談者は家系図の3番の人。
この手の問題は条件付き確率で考えるのが重要。いま、病気であることをD^+、遺伝素因を持っていること、すなわり常染色体優性遺伝のアレルを持っていることをA^+とする。
設問の条件は、「相談者が発症していない」ということで、問題の「第1子が発症する確率」は、
大親(1)のアレルが親(3)に遺伝して親(3)がアレルを持っていて、なおかつそれが子(5)に遺伝して、なおかつ発症しなければならない
とならないと発症しない。
ここで大親(1)から親(3)にアレルが伝わらなかったから親(3)が発症していないのか、それとも大親(1)から親(3)に遺伝したけれども浸透率のせいで発症していないか、そこが問題だ。
ベイズ的に真面目に解くと、親(3)について、病気が発症していないことがわかっていて、それで病気になるアレルを有する条件付き確率P_3(A^+|D^-)は、浸透率を\alphaとすると
P_3(A^+|D^-)=\frac{P_3(A^+\cap D^-)}{P_3(D^-)}=\frac{0.5(1-\alpha)}{1-0.5\alpha}=\frac{1-\alpha}{2-\alpha}
(大親(1)からアレルが遺伝して、かつ発症しなかった)/(発症しないのは、アレルが遺伝して発症した場合の余事象)となる。
これが親(3)がアレルを持つ確率である。ここから、子(5)にアレルが遺伝した場合に、確率\alphaで発症するから、\alpha=0.5であるから\frac{1}{3}*\frac{1}{2}*\frac{1}{2}=\frac{1}{12}となる。
 
R的シミュレーションはこちら。

library(kinship)
id <- 1:5
Mother <- c(0, 0, 2, 0, 3)
Father <- c(0, 0, 1, 0, 4)
Sex <- c(1, 2, 2, 1, 3)
affected <- c(2, 1, 1, 1, 1)
ped <- pedigree(id=id, dadid=Father, momid=Mother, sex=Sex, affected=affected)
plot(ped)
niter <- 100000
p1 <- c(0.5, 0.5) # 遺伝確率。遺伝する vs 遺伝しない
p2 <- c(0.5, 0.5) # 浸透率。発症する vs 発症しない
res <- matrix(0, niter, 4)
colnames(res) <- c("F1allele", "F1phenotype", "F2allele", "F2phenotype")
res <- as.data.frame(res)
res$F1allele <- sample(1:0, size=niter, prob=p1, replace=TRUE)
res$F1phenotype <- unlist(mapply(function(x) sample(c(res$F1allele[x], 0), size=1, prob=p2), seq(niter)))
res$F2allele <- unlist(mapply(function(x) sample(c(res$F1allele[x], 0), size=1, prob=p1), seq(niter)))
res$F2phenotype <- unlist(mapply(function(x) sample(c(res$F2allele[x], 0), size=1, prob=p2), seq(niter)))
# 親(3)が発症していない、という状況だけを考えればよい。
colMeans(subset(res, F1phenotype == 0))
   F1allele F1phenotype    F2allele F2phenotype 
 0.33334223  0.00000000  0.16767906  0.08432014 

親(3)がアレルを持つ確率は\frac{1}{3}、子(5)が発症する確率は\frac{1}{12}となった。

Scan statistics

MikuHatsune2013-07-04

Scan statisiticsというのを聞いたのだが、空間統計学の一種らしい。

Spatial scan statistics are used to determine hotspots in spatial data, and are widely used in epidemiology and biosurveillance. In recent years, there has been much effort invested in designing effcient algorithms for finding such "high discrepancy" regions, with methods ranging from fast heuristics for special cases, to general grid-based methods, and to efficient approximation algorithms with provable guarantees on performance and quality.

空間統計学は一般的には疫学によく用いられるらしい。コレラと井戸の話はよく出てくるが、この scan statistics を遺伝学的に使うのはどういうことかというと、数直線状に伸ばしたDNA鎖に対して、例えばレアバリアントとなるものが乗っていて、その hotspot 具合を ある幅をもった window をチマチマ動かしてその window 上に hotspot がいくつあるかをカウントするらしい。
時系列があればburst modelに似てないこともないような気がするけど、数直線的(空間的)な密度で考えると window になるんだろう。
 
Rでは DCluster, SpatialEpiといったパッケージが scan statistics に使えるらしい。

x1 <- runif(5,  min=-5, max=-1)
x2 <- runif(50, min=-1, max=1) # hotspot
x3 <- runif(5,  min=1 , max=5)
x0 <- sort(c(x1, x2, x3))

Wlen <- 1 # window の幅
Wstart <- seq(min(x0) - Wlen, max(x0) + Wlen, by=0.5) # window の開始地点

plot(x0, rep(1, length(x0)), pch=16, cex=0.5, xlim=range(Wstart), ylim=c(0, 1), xlab="position", ylab="")
mapply(function(i) segments(Wstart[i], 1 - 0.06*i, Wstart[i] + Wlen, lwd=3), head(seq(Wstart), -1L)) # 各window
i <- 8 # ある window の例
segments(Wstart[i], 1 - 0.06*i, y1=1, lty=2)
segments(Wstart[i] + Wlen, 1 - 0.06*i, y1=1, lty=2)

# 各 window に含まれる spot の数
spots <- mapply(function(i) sum(x0 >= Wstart[i] & x0 <= Wstart[i] + Wlen), head(seq(Wstart), -1L))
mapply(function(i) text(Wstart[i], 1 - 0.06*i, spots[i], adj=c(1.3, NA)), head(seq(Wstart), -1L))

plot(head(Wstart, -1L) + Wlen/2, spots, type="s", lwd=3, xlab="position", ylab="# of spots")


Ka/Ks (dN/dS) の計算

MikuHatsune2013-06-21

ある塩基配列の変異パターンを解析しているのだが、それでKa/Ks (dN/dS)比というモデルを勉強した。
しかしこれには計算法がたくさんあるらしい。
サンプルデータは http://evolution.genetics.washington.edu/book/primates.dna をコピペして使う。

library(phangorn)
library(seqinr)
primates <- read.phyDat("http://evolution.genetics.washington.edu/book/primates.dna", format="phylip", type="DNA")
tree0 <- NJ(dist.ml(primates))
plot(tree0)


 
Miyata & Yasunaga. J Mol Evol. 1980 Sep;16(1):23-36.
Nei & Gojobori. Mol Biol Evol. 1986 Sep;3(5):418-26.
Sqdif Plotで両者計算できる。計算結果もコピペで取ってこれるが、複数の pairwise alignment をやったときには一括ダウンロードできるようにしてほしい。
 
Li. J Mol Evol. 1993 Jan;36(1):96-9.
seqinrの kaks 関数に実装されているアルゴリズム
解説
縮退 (degenerate)というコドンの特徴があるが、これを
L_0:すべてのありうる変化が非同義
L_2:変化のうち 1/3 が同義
(L_3:イソロイシンのコドン3番目の場合。L_2とみなしてしまう)
L_4:すべてのありうる変化が同義
A_iL_i;\hspace{3}i=\{0,2,4\}に対応するTransition (purine/pyrimidine が変わらない変異)の数
B_iL_i;\hspace{3}i=\{0,2,4\}に対応するTransversion (purine/pyrimidine が変わる変異)の数
と定義して、
Ks=\frac{L_2A_2+L_4A_4}{L_2+L_4}+B_4
Ka=A_0+\frac{L_0B_0+L_2B_2}{L_0+L_2}
を計算する。

# コドンのため3の倍数の塩基長になっていないとエラー
k0 <- kaks(as.alignment(length(primates), seq=apply(as.character(primates)[,1:231], 1, paste, collapse="")))

 
最尤法 Mol Biol Evol. 1998 Dec;15(12):1600-11.
i番目のコドンとj番目のコドンについて、q_{i,j}
q_{i,j}=0:if i and j di er in more than one position
q_{i,j}=\pi_j:for synonymous transversion
q_{i,j}=\pi_j\kappa:for synonymous transition
q_{i,j}=\pi_j\omega:for nonsynonymous transversion
q_{i,j}=\pi_j\omega\kappa:for nonsynonymous transition
とする。ここで、\omega、は Ka/Ks,\kappaは Transition と Transversion の比、\pi_jはコドンjの平衡頻度である。\omegaが1かうんたらかは前も言ったとおり
phangornパッケージの pml 関数でできる。optim.pml で type を codon_n の何に指定するかで、\omega\kappaを1に固定するか自由に推定するか決められる。
サンプルはこちらを参考にしたが、サンプルがないわコピペ仕様になってないわで苦労した。

dat <- phyDat(as.character(primates), type="CODON")
fit0 <- pml(tree0, dat)
fit00 <- optim.pml(fit0, model="codon0", control=pml.control(trace=0))
fit1 <- optim.pml(fit0, model="codon1", control=pml.control(trace=0))
fit2 <- optim.pml(fit0, model="codon2", control=pml.control(trace=0))
fit3 <- optim.pml(fit0, model="codon3", control=pml.control(trace=0))
anova(fit00, fit2, fit3, fit1)
# 同じ結果
Likelihood Ratio Test Table
  Log lik. Df Df change Diff log lik. Pr(>|Chi|)    
1  -2905.6 25                                       
2  -2385.7 26         1       1039.65     <2e-16 ***
3  -2292.6 26         0        186.27     <2e-16 ***
4  -2291.4 27         1          2.39     0.1218    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

pairwise で比較して Ka/Ks を求めたいのだが、なぜかできない。
というわけでゴリ押ししようと思ったら、pairwise では optEdge=FALSE にしないと unroot がどうたらこうたら言われてエラーになる。
NJ法のところもエラーになってRがオシャカになるので、他の方法で距離行列を求める。

p0 <- as.character(primates)[, 1:231]
res <- diag(0, nrow(p0))
dimnames(res) <- list(rownames(p0), rownames(p0))
cmb <- combn(nrow(p0), 2)
for(i in seq(ncol(cmb))){
	tmp <- p0[cmb[, i], ]
	tree0 <- as.phylo(hclust(dist.ml(as.phyDat(tmp)))) # NJ だと壊れる
	dat <- phyDat(p0, type="CODON")
	fit0 <- pml(tree0, dat)
	# pairwise だと optEdge=FALSE にしないといけないっぽい
	pml1 <- optim.pml(fit0, model="codon1", optEdge=FALSE) # k and omega are not fixed to 1
	res[cmb[2, i], cmb[1, i]] <- pml1$dnds
}

系統樹を作成するときにもいろいろなアルゴリズムが使える。

# Maximau likelifood
# pairwise だと optEdge=FALSE にしないといけないっぽい
pml1 <- optim.pml(fit0, model="codon1", optEdge=FALSE) # k and omega are not fixed to 1
pml1$dnds

# Jukes-Cantor (starting tree from NJ)
fitJC <- pml(tree0, dat)

# optimize edge length parameter
fitJC <- optim.pml(fitJC, model="codon1", optEdge=FALSE)

# JC + Gamma + I - model
fitJC_GI <- update(fitJC, k=4, inv=0.2)

# optimize shape parameter + proportion of invariant sites
fitJC_GI <- optim.pml(fitJC_GI, optGamma=TRUE, optInv=TRUE, optEdge=FALSE)

# GTR + Gamma + I - model
fitGTR <- optim.pml(fitJC_GI, optBf=TRUE, optQ=TRUE, model="codon1", optEdge=FALSE)

# Pamilo-Bianchi-Li
k0 <- kaks(as.alignment(2, seq=apply(p0, 1, paste, collapse="")))
k0$ka / k0$ks

 
20171225追記
dnds の時系列の解析や検定ができないのかと聞かれたので調べものメモ
https://www.ncbi.nlm.nih.gov/pubmed/23650353
https://www.ncbi.nlm.nih.gov/pubmed/25388108
https://www.ncbi.nlm.nih.gov/pubmed/24748652
https://www.ncbi.nlm.nih.gov/pubmed/12509511
https://en.wikipedia.org/wiki/McDonald%E2%80%93Kreitman_test

SNP-set (Sequence) Kernel Association Test (SKAT)

MikuHatsune2013-06-17

SKATという手法があるらしい。
kernel methodというなんとも中二病的ネーミングの手法を用いて、レアバリアントがどうたらこうたら言っていたがちょっと追いつけなかった。
SKATをやってみる。
昔やったロジスティック回帰をやると

library(SKAT)
data(SKAT.example)
attach(SKAT.example)

g0 <- glm(y.b ~ y.c, data=SKAT.example, family="binomial")
p <- function(x){
	1/(1+exp(-sum(g0$coefficients * c(1, x))))
}
s0 <- seq(-4, 4, length=1000)
plot(SKAT.example$y.c, SKAT.example$y.b, xlim=range(s0))
points(s0, mapply(p, s0), type="l", col=2, lwd=3)

閾値を決めて表現型が決められているようだ。

# Compute the P-value of SKAT with default Beta(1,25) Weights
# - without covariates
# continuous trait
obj <- SKAT_Null_Model(y.c ~ 1, out_type="C")
SKATc0 <- SKAT(Z, obj)
# dichotomous trait
obj <- SKAT_Null_Model(y.b ~ 1, out_type="D")
SKATd1 <- SKAT(Z, obj)

# Compute the P-value of SKAT with default Beta(1,25) Weights
# - with covariates
# continuous trait
obj <- SKAT_Null_Model(y.c ~ X, out_type="C")
SKATc1 <- SKAT(Z, obj)
obj.b <- SKAT_Null_Model(y.b ~ X, out_type="D")
SKATd1 <- SKAT(Z, obj.b)

# Compute the P-value of SKAT with default Beta(1,25) Weights
# - Optimal Test
SKAT2 <- SKAT(Z, obj, method="optimal.adj")
SKAT3 <- SKAT(Z, obj.b, method="optimal.adj")

# Compute the P-value of SKAT with Beta(1,30) Weights
SKAT3 <- SKAT(Z, obj, weights.beta=c(1,30))