2013-04-22

Genotypeに影響されるPhenotype分散

MikuHatsune2013-04-22

遺伝子座による形質の定量度の違いをquantitative trait loci (QTL)というらしい。
普通、ある遺伝子座がA→Gだと、収縮期血圧が10mmHg高い傾向にある、みたいな、平均値の変動を定量化することが多いが、分散の変動、つまり、普通アレルだと平均±5mmHgのところが、このアレルだと±10mmHgでばらつく、みたいな影響のことを variance QTL というらしい。
Shen X., et al. PloS Genet. 2012.では、Arabidopsis thaliana(シロイヌナズナ)のデータベースから、モリブデン輸送体MOT1のgenotypeとモリブデン濃度のばらつきを示している。
論文Figure 3より引用

　
vQTLについての手法はRönnegård L. BMC Genet 2012.にいくつかまとまっている。
論文Figure 1より引用。

　
Rではqtl, QTLRelパッケージが検索でひっかかるが、中身をじっくり読むヒマがないので、プロット的なことをしてくれるものをやってみたら、マンハッタンプロット的なものが出てきたけどよくわからないので保留。

data(miscEx)
# impute missing genotypes
gdat.imp <- genoImpute(gdatF8, gmap=gmapF8, step=Inf, gr=8, na.str=NA)
# estimate variance componentsqqPlot 27
o <- estVC(y=pdatF8$bwt, x=pdatF8$sex, v=list(AA=gmF8$AA,DD=gmF8$DD, HH=NULL, AD=NULL, MH=NULL, EE=diag(length(pdatF8$bwt))))
# genome scan
llk <- scanOne(y=pdatF8$bwt, x=pdatF8$sex, vc=o, gdat=gdat.imp)
# plotting
plot(llk, gmap=gmapF8) # gmap is needed
# plotting in another way
idx <- match(colnames(gdat.imp), gmapF8$snp)
tmp <- data.frame(chr=gmapF8$chr[idx],dist=gmapF8$dist[idx],y=llk$p)
plotit(tmp, main="Mapping Plot", xlab="Chromosome", ylab="LRT", col=as.integer(tmp$ch)%%2+2,type="p")

2013-04-15

塩基配列上での変異とアミノ酸置換その2

R 数理モデル統計遺伝医学医学 Rを使いこなす

MikuHatsune2013-04-15

塩基配列のComplexityとしてエントロピーという統計量を使うのだが、このエントロピーの定義が自分の周辺で物議を醸している。
一般的なシャノンエントロピーは
$S_1=-\sum_{i\in \{A,G,C,T\}}p_i log_2 p_i$
だが、前の論文では、HBVのquasispeciesとして、クローン数が全部で $N$ 個の配列で、配列 $s_i$ が $p_i$ の割合で存在しているとすると
$S_2=-\sum_i^n \frac{p_i log{p_i}}{log N}$
となる、と書いてある。
　
$S_1$ は、各塩基位置において、AGCTがどの頻度で出現しているかを計算している。
$S_2$ は、各(臨床)サンプルにおいて、ある配列をもつクローンがどの頻度で出現しているかを計算している。
と思う。
で、これらが一致するかという質問を受けたのだが、結論から言うと一致しない。
$S_1$ は、取りうる要素数が $C(=4)$ (塩基)のとき、 $0\leq S_1\leq log_2 C$ となる。
$S_2$ は、~~未検証~~。 $S_1$ が $0\leq S_1\leq log_2 C$ となるので、 $log N$ (エントロピーを自然対数で計算しなおせば)で割れば $0\leq S_2\leq 1$ になる。

library(MCMCpack)
p1 <- rdirichlet(10000, rep(1, 4))
plot(sort(rowSums(-log(p1, 2) * p1)), ylab="Entropy")

　
例えば、10塩基長、15クローンのデータがあったとして、このデータのシャノンエントロピー( $S_2=0.10$ )となる。
各塩基でのエントロピー $S_1$ は、総和が2.65で平均が0.265なので、どうすれば一致するんだろうか。たぶんしない。

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [2,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [3,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "C"  "G"  
 [4,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [5,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "T"  "A"  "G"  
 [6,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "G"  "G"  
 [7,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [8,] "C"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [9,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[10,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[11,] "T"  "A"  "G"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[12,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[13,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[14,] "T"  "A"  "T"  "C"  "G"  "T"  "T"  "A"  "T"  "G"  
[15,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "C"  "G"

library(seqinr)
library(bio3d)
bp <- 10
clone <- 15
niter <- 1000
base <- c("A", "G", "C", "T")
par(mfrow=c(1, 2))
for(m0 in c(5, 50)){
	mut_bases <- rpois(1, m0)
	res <- matrix(0, niter, 2)
	for(i0 in seq(niter)){
		s1 <- sample(base, size=bp, replace=TRUE)
		m1 <- mapply(rep, s1, clone)
		mut <- sample(seq(length(m1)), mut_bases)
		mutbase <- sample(base, size=length(mut), replace=TRUE)
		m2 <- matrix(replace(c(m1), mut, mutbase), nr=clone)
		
		#AGCTベースのエントロピー
		Sn_p1 <- mapply(function(x) table(factor(m2[ ,x], levels = base)), seq(bp)) / clone
		Sn1 <- -colSums(Sn_p1 * log(Sn_p1, 2), na.rm=TRUE)
		names(Sn1) <- seq(bp)
		#barplot(Sn1, ylim=c(0, log(length(base), 2)), xlab="base", main="Shannon entropy")
		
		#クローンベースのエントロピー
		seq0 <- table(apply(m2, 1, paste, sep="", collapse=""))
		Sn2 <- -(seq0 / clone) * log(seq0 / clone) / log(clone)
		
		res[i0, ] <- c(mean(Sn1), sum(Sn2))
	}
	xl <- yl <- range(0, res)
	plot(res, xlim=xl, ylim=yl, cex=0.7, xlab="Base", ylab="Clone")
	title(paste("Correration of Shannnon entropy\nMutation in", round(m0/(bp*clone) * 100, 1), "% in sequence"))
	abline(lm(res[, 2] ~ res[, 1]))
}

あるひとつの配列中で、変異する塩基が増えると、クローンのシャノンエントロピー( $S_2$ )はある定数に近づくようである。

　
bio3dパッケージにあるentropy関数は、塩基エントロピー( $S_1$ )で計算しているようである。
このなかのサンプルデータとして、HIVのプロテアーゼのどっかの領域にひっかかる。
protein BLASTより、データの1行目

PQITLWQRPLVTIKIGGQLKEALLDTGADDTVLEEMSLPGRWKPKMIGGI--GGFIKVRQYDQILIEICGHKAIGTVLVGPTPVNIIGRNLLTQIGCTLNF

を検索すると出てくる。

これのサンプルコードを実行すると、いろいろ計算して描いてくれるが、各アミノ酸残基でのエントロピーを計算している。

# Read HIV protease alignment 
aln <- read.fasta(system.file("examples/hivp_xray.fa",package="bio3d"))
# Entropy and consensus
h   <- entropy(aln)
con <- consensus(aln)
names(h$H)=con$seq
print(h$H)
H <- h$H.norm

col <- mono.colors(100)
aa  <- rev(rownames(h$freq))
oldpar <- par(no.readonly=TRUE)
layout(matrix(c(1, 2), 2, 1, byrow = TRUE), widths = 7, heights = c(2, 8), respect = FALSE)

# Plot 1: entropy
par(mar = c(0, 4, 2, 2))
barplot(H, border="white", ylab = "Entropy", space=0, xlim=c(3.7, 97.3),yaxt="n" )
axis(side=2, at=c(0.2, 0.4, 0.6, 0.8))
axis(side=3, at=(seq(0, length(con$seq), by=5) - 0.5), labels=seq(0, length(con$seq), by=5))
box()

# Plot2: residue frequencies 
par(mar = c(5, 4, 0, 2))
image(x=1:ncol(con$freq),
   y=1:nrow(con$freq),
   z=as.matrix(rev(as.data.frame(t(con$freq)))),
   col=col, yaxt="n", xaxt="n",
   xlab="Alignment Position", ylab="Residue Type")
axis(side=1, at=seq(0,length(con$seq), by=5))
axis(side=2, at=c(1:22), labels=aa, cex.axis=0.7)
axis(side=3, at=c(1:length(con$seq)), labels =con$seq, cex.axis=0.4, tick=0)
axis(side=4, at=c(1:22), labels=aa, cex.axis=0.7)
grid(length(con$seq), length(aa))
box()
for(i in 1:length(con$seq)) {
	text(i, which(aa==con$seq[i]),con$seq[i],col="white")
}
abline(h=c(3.5, 4.5, 5.5, 3.5, 7.5, 9.5, 12.5, 14.5, 16.5, 19.5), col="gray")
par(oldpar)

2013-04-04

塩基配列上での変異とアミノ酸置換

R 遺伝医学数理モデル Rpackage Rを使いこなす

MikuHatsune2013-04-04

J Hepatol. 2009 May;50(5):895-905
HBVの薬剤応答具合と逆転写酵素領域の配列の変化をモデル化して考えた論文。
治療開始後48週間後では、反応が良かった群で治療開始前の塩基配列と比べて複雑度と多様度が減少している、ということらしい。
その数理モデルをシミュレーションしてみる。
　
Complexity
シャノンエントロピー
クローン数が全部で $N$ 個の配列で、配列 $s_i$ が $p_i$ の割合で存在しているとすると
$S_n=-\sum_i^n \frac{p_i log{p_i}}{log N}$
と計算される。
　
Diversity
Hamming距離
ベクトルのうち、異なる要素があればそこの距離を1にしてしまうやり方。
配列にしてしまえば、異なる文字であれば1にして、和をとる。
vwrパッケージでできる。

library(vwr)
slen <- 20
set.seed(831)
s0 <- paste(sample(c("A", "G", "C", "T"), size=slen, replace=TRUE), sep="", collapse="")
s1 <- paste(sample(c("A", "G", "C", "T"), size=slen, replace=TRUE), sep="", collapse="")
s0; s1; hamming.distance(s0, s1)

[1] "CGTTACGATGTGATTAACCC"
[1] "GTCGAGGTAATAATTGGGTA"
GTCGAGGTAATAATTGGGTA
                  14

　
$\frac{Ka}{Ks}$ 比というモデルがある。
これは、塩基配列中のコドンでの塩基変異が、アミノ酸置換を起こす(non-synonymous, Ka or dN)か起こさないか(synonymous, Ks or dS)という塩基配列上でのポジションごとにいろいろ計算する、というもの。
ググっても解説してくれているサイトはあまりない気がする。検索するとけっこう上位に出てくる。
Rではseqinrパッケージで計算できる。

　
いま、適当な長さの配列を発生させたとする。
この配列の長さは、タンパクに翻訳されるように3の倍数にしてある。
この配列に、配列の長さのうち%の塩基がこれまた適当に塩基変異するとする。ただし、G→Gも塩基変化としてしまうことにする。
塩基変化する%をいろいろ変えてみて、KaとKsがどう変化するかをシミュレーションした。

library(seqinr)
# seqinr のexample
s <- read.alignment(file = system.file("sequences/test.phylip", package = "seqinr"), format = "phylip")
s0 <- kaks(s)
aa <- lapply(lapply(lapply(s$seq, s2c), translate), paste, sep="", collapse="") # s2c で文字列をバラバラにして、translate で翻訳する。
aa_aln <- as.alignment(5, s$nam, aa) # alignment でalignment オブジェクトにする。

AAs <- 100 # アミノ酸数
#HS <- s2c(s$seq[[3]]) 
HS <- sample(c("a", "c", "g", "t"), size=3*AAs, replace=TRUE) # 配列を適当に作成。
n1 <- 100 # 変異を起こす配列数
mutation_base <- length(HS) # 変異を起こす塩基数
niter <- 1 # シミュレーションしようと思っていたけど作成する配列数を増やせばいいことに気づいた。
res <- array(0, c(mutation_base, n1 - 1, 2)) #ka と Ks
ent <- rep(0, mutation_base) # シャノンエントロピー
for(m0 in seq(mutation_base)){
	mHS <- vector("list", n1 - 1)
	for(i in seq(mHS)){
		#m0 <- 2 #sample(length(HS), 1)
		mHS[[i]] <- replace(HS, sample(seq(HS), m0), sample(c("a", "g", "c", "t"), m0, replace=TRUE))
	}
	mHS <- lapply(c(list(HS), mHS), paste, sep="", collapse="")
	t0 <- table(unlist(mHS))
	ent[m0] <- -sum((t0/sum(t0))*log(t0/sum(t0))/sum(t0))
	s1 <- kaks(as.alignment(length(mHS), seq=mHS))
	res[m0, , 1] <- tail(as.matrix(s1$ka)[1, ], -1L)
	res[m0, , 2] <- tail(as.matrix(s1$ks)[1, ], -1L)
}
par(mfcol=c(1, 2))
for(i in seq(2)){
	boxplot(t(res[ , , i]), xaxt="n", xlab="proportion of mutated bases", ylab="value", ylim=c(0, 4), cex=0.5)
	title(c("Ka", "Ks")[i])
	axis(1, at=seq(0, 3*AAs, length=11), labels=seq(0, 1, length=11))
}

KaとKsはそれぞれ規格化されているので、適当な塩基変異の場合、 $\frac{Ka}{Ks}$ は1になる。
$\frac{Ka}{Ks}>1$ の場合、タンパク質の変化を伴うような変異、つまりnon-synonymousな変異が積極的に起こっていると考えられ、生き残りをかけているような正の選択がなされていると考える。これは、未知の細菌やウイルスといった抗原に対する免疫反応を多様化させるといったところで役に立つ。
$\frac{Ka}{Ks}<1$ の場合はその逆で、タンパク質変化を伴わないsynonymousな変異で保守的なんだろう。

図は配列中の変異する塩基数が多くなると、Ka, Ksともに大きくなる。最小値は0、最大値は10のようである。
最小値は定義上0以上のはずだがなぜか-0.3とかが混じっていた。謎である。
比をとれば1に近い感じになった。実際の塩基配列はGCリピートなどで偏っているだろうし、変異する塩基ポジションもホットスポットのような領域があるはずなので、 $\frac{Ka}{Ks}$ が1に近いことはあまりないのではないだろうかと思う。

2012-07-28

伊藤誠で作成しようと思ったら沢越止のほうが鬼畜だと気づいた

R 遺伝医学 Rを使いこなす

昔、家系図を作成したのだが、某業界で有名な一家の家系図をRで作成したいと思った。
ここやここを参考にデータ入力し、作図した。
もっと色づけとかやりたかったけど保留。
こんな家系の家族歴問診は激ムズだろう。

library(kinship)
school_days_cast <- seq(59)
Mother <- c(6, 9, 11, 0, 0, 0, 0, 0, 5, 5,
					0, 0, 0, 7, 8, 10, 10, 10, 11, 11,
					12, 12, 13, 0, 30, 14, 14, 14, 14, 9,
					0, 9, 17, 18, 25, 21, 21, 21, 24, 23,
					23, 0, 30, 30, 30, 42, 31, 32, 45, 49,
					49, 0, 52, 52, 50, 50, 55, 43, 30)
Father <- c(4, 1, 1, 0, 0, 0, 0, 0, 4, 4,
					0, 0, 0, 1, 1, 2, 2, 1, 1, 1, 
					1, 2, 1, 0, 1, 2, 2, 2, 1, 1,
					0, 1, 2, 2, 20, 20, 20, 1, 3, 3,
					3, 0, 1, 1, 1, 2, 2, 3, 1, 1,
					1, 0, 1, 1, 1, 1, 20, 3, 2)
Sex <- c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
				2, 2, 2, 2, 1, 2, 2, 2, 2, 1,
				2, 1, 2, 2, 2, 1, 1, 2, 1, 2,
				2, 2, 1, 2, 3, 1, 1, 3, 3, 2,
				2, 2, 2, 2, 2, 2, 1, 2, 2, 2,
				3, 2, 2, 1, 2, 3, 3, 2, 1)
Affected <- rep(1, length(school_days_cast))
Affected[1:3] <- 2 #沢越止、伊能歩、間瞬
Status <- rep(0, length(school_days_cast))
Status[c(48, 54, 58)] <- 1 #伊藤誠、西園寺世界、清澄刹那
cbind(school_days_cast, Mother, Father)

Makoto <- pedigree(
	id = school_days_cast,
	momid = Mother,
	dadid = Father,
	sex = Sex,
	affected = Affected,
	status = Status
	)
plot(Makoto)
plot(Makoto, col=2)

2012-01-16

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

R 遺伝医学

家系データをシミュレーションするセミナー
遺伝学をシミュレーションする上で必要な用語について

2011-07-27

第五回

医学 2011研究科横断統計遺伝学遺伝医学

今回の話
パラメーターで表すこと
モデルを立てる
パラメータ：本当に知りたいことを変数Xにおいてみる
Xを知るために：実験しましょう。観測するのは別のこと。
X←Y（数えられるもの）：何かしらの関係

振り返ってみて
塩基
X：種間の差　Y：塩基配列の違い
YがなぜXに関連づけられるか?：変異の個数
配列が同一だと、種は同じ
変異が、自然選択の結果起きたのだろう：変異の個数、変異率（数/時間）が種が分岐した瞬間だろう
2つの異なる遺伝子
種間を隔てている時間

siの実験
発現量と反応
X：ストレス応答　Y：生細胞率
ストレス応答があがると、生細胞数あがる
ストレス応答が下がると、生細胞数あがる

配列の違いと変異率、ストレス応答と生存率の関係は、少なくとも単調増加（減少）だろうと思われる。（簡単なモデル）
面倒くさい：山とか多峰性とか

t.testについて
si(-)とsi(+)で生存率を比較
どう差があると思っている?
si(+) だと、生存率80%くらいだとしよう。
80%付近に多くなってほしい：正規分布

発現量を落とすと、生存率80%になる。いろいろになると思うけど、その周辺に分布するとします（モデル化）。
対照は、生存率75%でその付近に分布。
（正規分布は中央と分散を仮定している。）
というモデル。

両者を合わせて、中央77.5%、分散‥の一山と仮定するより、
中央80%、分散‥、中央75%、分散‥の二山あると仮定したほうがいいんじゃない、というモデル。
これを判断してくれるのが、t.test
ちょっと考える：生存率が100% or 0%に近いとき
正規分布を考えるモデルが違うのではないか?
死に絶えるのがほとんどで、ほんのちょっとだけ生きる、とか

一番おきやすいことは?
一山にした生存率78%か
二つ山がある80%と75%か?

モデル
前回のSNPと優性/劣性遺伝モデル
優性と劣性との間に、相加モデルを入れる。
A/A 2
A/G 1
G/G 0
として
優性：A/A and A/G 1 G/G 0
劣性
遺伝型に点数を与える
A/A A/G G/G
優性 2 2 0
相加 2 1 0
劣性 2 0 0
何が違う?A/G列：ヘテロの点数の与え方
2は扱いにくいので1に変換
A/A A/G G/G
優性 1 1 0
相加 1 0.5 0
劣性 1 0 0
多型のアレルの強さ
ヘテロのときの強さを、2つのホモタイプのどの辺におくか。
相加平均にする：相加モデル
偏らせる：優性/劣性モデル
A/Gの力がG/Gの半分か?：実のところ、そういうことはあまりない
そうなれば、0.8とか‥‥

ヘテロの力がx(0

2011-07-20

第四回　多重検定とp値補正

R 2011研究科横断統計遺伝学遺伝医学医学

前々回
DNA配列比較
コーディング配列とRNA配列の評価は違う（かも）
　
前回
細胞に介入　タンパクの発現量を変える
表現型に違いが出るか
介入したかしないか、と表現型、を比較しているのではなくて
発現量、と表現型、の比較になっている→量的関係の回帰直線
　
横軸が、だんだん変化するタンパク発現量とかなら、回帰直線とかできる
横軸が、si(-)、si.1、si.2　（タンパク発現量は定量していない）とかいった固有名詞ならどうする？
りんご、みかん、ばなな、みたいな
りvsみ
みvsば
ばvsり　でそれぞれp値が出る
k群でペアワイズ検定： $_kC_2=\frac{k(k-1)}{2}$ 通り
　
たくさんやらなければならない（やりたい）ときはどうする？
　
たくさんやるときの弊害‥の前に
検定するなら
方針1
すべてのペアでやる
変な小細工をするとややこしい。シンプルに、ルールは少なく。
　
何回か検定を繰り返すと、閾値を超えるp値が出る確率が大きくなる：素のp値を使ってコミュニケーションしてはいけない
ずるい
$\textrm{p<0.05}$ になるのは、珍しい、と思っている→「帰無仮説」を信じたままだと珍しい
　
「帰無仮説」を信じるとは?
[tex:\textrm{0