Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2018-06-11

就学年数が長いから近視になるのか近視になるから就学年数が長いのか

読んだ。

BMJ. 2018 Jun 6;361:k2022.

観察研究では、どちらも関連があるが、この論文の結論から言えば就学年数が長いから近視になり、近視だから就学年数が長い、ということにはならない。これを近視に関する多型でメンデリアンランダマイゼーションした場合(Nat Genet. 2016 Jul;48(7):709-17.)と、就学年数に関する多型でメンデリアンランダマイゼーションした場合(Nature. 2016 May 26;533(7604):539-42.)で解析し、因果関係を推定している。

 

メンデリアンランダマイゼーションを双方向に行って因果関係を推定するのは、

例えばCRPBMI はどちらも高いことが多いが、これはCRP の多型でランダマイズしてもBMI は増加しないが、BMI に関係するFTO という多型でランダマイズするとCRP が増加する、ということでやっている。

Int J Obes (Lond). 2011 Feb;35(2):300-8.

2018-05-15

MikuHatsune2018-05-15

DNA型が一致しないときに突然変異ということで一致判定できるか

こんな記事を見かけた。

新規サイト 共通エラーページ

判例

平成29(あ)882

困惑している人が多い。

意外と知られていない知識「突然変異によるDNA型の変化はまれにある」裁判での逆転有罪の科学的根拠に驚きの声 - Togetter

 

※第二審判決には細かいallele と変異パターンが書いてあったけど頻度表が違うのでこれ以上はまあいいや…

 

最初にこのタイトルだけ聞いたときは、確かに(えっそうなん?)と思ったが、生物学的および法医学的原理を考えれば、突然変異という結果になる。

なる、というのは正しくないが、「突然変異」と「違う人物(この裁判の場合は、第三者DNAの混入)」を比較するならば、「突然変異」>>>(超えられない壁)>>>「第三者」となる、ということ。

 

前提として、DNA鑑定は、STR と呼ばれる小リピート配列領域を持って、DNA個人の結びつきを考える。ここらへんはわかりやすい、かも

DNA鑑定による法医学的個人識別の確率・統計学的背景

分子生物学的な話がわからない人はググってほしいが、ある領域を見れば、集団において「あるパターン」を持つ人は、○パーセントである、というようなことが既にわかっている。

この、ある領域をいくつも採用すれば、(基本的に各領域は独立なので)掛け合わせでどんどんそのパターンを持つ個人は減っていく、という理論である。

 

日本人の15 STR 領域の集団頻度論文では、警視庁でも使っている15STR の頻度を報告しているので、これを使ってみる。

15 STR (marker)と、そのとり得るパターン(allele)の頻度は以下のとおりである。

f:id:MikuHatsune:20180515210204p:image

 

さて、この頻度から適当に一人の人Aを作った。

     D3S1358 TH01  D21S11 D18S51 PentaE D5S818 D13S317 D7S820 D16S539 CSF1PO PentaD vWA  D8S179 TPOX FGA 
[1,] "16"    "9.3" "32.2" "15"   "5"    "9"    "8"     "11"   "10"    "12"   "9"    "17" "12"   "11" "19"
[2,] "17"    "9"   "32.2" "18"   "17"   "11"   "10"    "10"   "11"    "9"    "10"   "18" "14"   "8"  "22"

このとき、各マーカーのalelle頻度は

   D3S1358   TH01 D21S11 D18S51 PentaE D5S818 D13S317 D7S820 D16S539 CSF1PO PentaD    vWA D8S179   TPOX    FGA
16  0.3110 0.0305 0.1494 0.2195 0.1098 0.0762  0.2957 0.3323  0.2043 0.3932 0.3293 0.2713 0.1067 0.3933 0.0518
17  0.2287 0.4543 0.1494 0.0610 0.0976 0.2348  0.0945 0.2012  0.1860 0.0305 0.2317 0.2012 0.2195 0.4329 0.2043

なので、これをすべてかけあわせると、このパターンのalleleを持つ人Aの確率は

[1] 8.430659e-24

となる。

 

さてここで判例では、

(K鑑定では)精液検査により,多数の精子を認めた一方,精子以外の特異な細胞が見当たらず,また,STR型検査等により検出された15座位のSTR型とアメロゲニン型が被告人の口腔内細胞のものと一致した。

また、

(S鑑定では)それぞれ14座位のSTR型とアメロゲニン型が科捜研鑑定と一致したものの,1座位で,科捜研鑑定と合致する2つのSTR型に加え,これと異なる3つ目のSTR型を検出した。

とある。

要は、STR マーカーの14箇所のallele は同じだったが、1箇所(どのSTR マーカーなのか、そしてそれが何に変わっていたのかは不明)ということのようである。

 

ここで、この裁判で重要なのは、この1座位の違いは、「¥Theta_1突然変異から生じた」のか、「¥Theta_2第三者DNAが混入した」のか、が争点である。仮説を¥Thetaとした。ここで、この仮説というのは仮説なので、どちらが正しい(100%なのか)というのはわからないが、原理系に計算できるので、「¥Theta_1¥Theta_2の仮説の大小を比べましょう」ということになる。これは、仮説の比みたいなことを言って¥frac{¥Theta_1}{¥Theta_2} を考えることになる。両者の仮説が等しいならば1 だし、片方が0 もしくは1 でない限り、[0, ¥infty) の間になる。

というわけで頻度と、14STR が共通、という前提の上で、「¥Theta_1突然変異から生じた」と「¥Theta_2第三者DNAが混入した」を比較する。

いま、適当に、

[1] "D16S539"

が異なるSTR型を有していたとする。このSTR を選んでallele を取り替えた人A^{’}

     D3S1358 TH01  D21S11 D18S51 PentaE D5S818 D13S317 D7S820 D16S539 CSF1PO PentaD vWA  D8S179 TPOX FGA 
[1,] "16"    "9.3" "32.2" "15"   "5"    "9"    "8"     "11"   "15"    "12"   "9"    "17" "12"   "11" "19"
[2,] "17"    "9"   "32.2" "18"   "17"   "11"   "10"    "10"   "11"    "9"    "10"   "18" "14"   "8"  "22"

である。10 が15 になっている。このとき、AA^{’} が異なる人、つまり、「¥Theta_2第三者DNAが混入した」という確率は、14STR が一致しているA^{’}一般集団からサンプリングできる確率なので、2.21905e-22 になる。

一方、1STR が突然変異したからこうなった、と考える「¥Theta_1突然変異から生じた」は、突然変異を生じる確率は10^{-5} 程度らしい。

Genomics Proteomics Bioinformatics. 2007 Feb;5(1):7-14.

もちろん、突然変異しやすい/しにくいSTR や、変化しやすい/しにくいパターンがあるので一概には言えないが、これも仮説に取り込めばよい。

というわけで、¥frac{¥Theta_1}{¥Theta_2}=4.5¥times 10^{17} となる。

これは単純に、「¥Theta_1突然変異から生じた」のは「¥Theta_2第三者DNAが混入した」に比べて4.5¥times 10^{17} 倍起きやすい、となる。

 

いまは特定のSTR allele パターンで、変異したSTR もたった一つ検討しただけなので、たくさんシミュレーションして試してみる。

少なく見積もっても10^{10} 倍は、「¥Theta_1突然変異から生じた」と考えるほうが妥当である。

ここで、¥frac{¥Theta_1}{¥Theta_2} の最小値(の常用対数)が10 くらいなので、¥frac{¥Theta_1}{¥Theta_2}=10^{3} くらいになる突然変異確率を考えると、10^{-12} くらいになるが、これは生殖細胞系列での突然変異頻度を考えると、あまりにも小さすぎるので、分子生物学的にはおかしい、となる。

f:id:MikuHatsune:20180515214726p:image

 

科学的に考えればこうなる。ただし、検体の取り違えや、突然変異パターンの偏りは考慮していない。

司法的な考え方では、「疑わしきは被告人の利益に」だと素人的には思うので、「¥Theta_2第三者DNAが混入した」確率が0 と証明できない以上、がんばることはできると思うが、それでも、第二審で間違った判断を下しても許されるアレは、はるかに短い時間と少ない情報で億単位の賠償リスクについて判断しないといけない、とある社畜業界で働く身としてはうらやましいなぁと思いました、まる。

freq <- read.csv("freq.txt", header=TRUE, row=1)
freq <- as.matrix(head(freq, -6))
freq <- sweep(freq, 2, colSums(freq, na.rm=TRUE), "/")

# 15 STR の頻度
cols <-jet.colors(100)
x <- seq(nrow(freq))
y <- seq(ncol(freq))
par(mar=c(5, 6, 2, 2))
image(x, y, freq, col=cols, xaxt="n", yaxt="n", xlab="", ylab="")
axis(1, at=x, labels=rownames(freq), las=2)
axis(2, at=y, labels=colnames(freq), las=2)
mtext("Allele", 1, line=3)
mtext("STR marker", 2, line=5)
abline(v=head(x, -1)+0.5, lty=3)
abline(h=head(y, -1)+0.5, lty=3)


freq[is.na(freq)] <- 0
allele <- rownames(freq)
marker <- colnames(freq)
# 各STR と、そのアレルの組み合わせを作る
dip <- dip1 <- mapply(function(z) sample(allele, size=2, prob=freq[,z], replace=TRUE), marker)

# どのSTR が変異していたか適当に選ぶ
mut <- sample(marker, size=1)

# その変異していたSTR がどのallele パターンに変異していたか選ぶ
change <- sample(setdiff(allele[freq[,mut]>0], dip[,mut]), size=1)
# 戻す
dip1[sample(2, size=1), mut] <- change

# 各STRとallele の頻度
genotype <- mapply(function(z) freq[dip1[, z], z], setdiff(marker, mut))

# diplotype として取りうる頻度
apply(genotype, 2, prod)

# 全部積にしてしまえば、その個人の確率
subject <- prod(genotype)
# 変異する確率
rate <- 10e-5

# シミュレーション
sim <- function(){
  dip <- dip1 <- mapply(function(z) sample(allele, size=2, prob=freq[,z], replace=TRUE), marker)
  mut <- sample(marker, size=1)
  change <- sample(setdiff(allele[freq[,mut]>0], dip[,mut]), size=1)
  dip1[sample(2, size=1), mut] <- change
  genotype <- mapply(function(z) freq[dip1[, z], z], setdiff(marker, mut))
  subject <- prod(genotype)
  rate <- 10e-5
  return(rate/subject)
}

iter <- 1000
res <- replicate(iter, sim())
xl <- c(-1, max(log10(res)))
p0 <- mapply(function(p) (1-p)/p, 10^seq(0, -10)[-1])
hist(log10(res), col=7, xlim=xl, main="", xlab="log likelihood ratio")
Allele,D3S1358,TH01,D21S11,D18S51,PentaE,D5S818,D13S317,D7S820,D16S539,CSF1PO,PentaD,vWA,D8S179,TPOX,FGA
5,,,,,0.1098,,,,,,,,,,
6,,0.2225,,,,,,,,,0.0031,,,0.0031,
7,,0.2439,,,,0.0061,0.0091,,,0.0031,0.0031,,,,
8,,0.0396,,,0.0092,0.0031,0.2957,0.1281,,0.0031,0.0213,,,0.4329,
9,,0.4543,,,0.0031,0.0762,0.122,0.0427,0.3445,0.0305,0.3293,,0.0031,0.1189,
9.3,,0.0305,,,,,,,,,,,,,
10,,0.0092,,,0.0396,0.2073,0.0945,0.2012,0.2043,0.2531,0.2317,,0.1372,0.0152,
11,,,,0.0031,0.1372,0.2348,0.2073,0.3323,0.186,0.189,0.186,,0.1281,0.3933,
12,0.0031,,,0.0396,0.1463,0.2409,0.1982,0.247,0.1829,0.3933,0.1189,,0.1067,0.0335,
13,,,,0.2134,0.0274,0.1646,0.064,0.0457,0.0671,0.1128,0.0793,,0.2226,,
14,0.0335,,,0.1951,0.0457,0.0671,0.0091,0.0031,0.0122,0.0122,0.0274,0.2378,0.2195,0.0031,
15,0.3476,,,0.2195,0.125,,,,0.0031,0.0031,,0.0183,0.1067,,
16,0.311,,,0.1098,0.0579,,,,,,,0.2043,0.0701,,
17,0.2287,,,0.0732,0.0976,,,,,,,0.2713,0.0061,,
18,0.0671,,,0.061,0.0793,,,,,,,0.2012,,,0.0396
19,0.0091,,,0.0305,0.0518,,,,,,,0.0488,,,0.0518
20,,,,0.0244,0.0335,,,,,,,0.0152,,,0.1128
21,,,,0.0183,0.0092,,,,,,,0.0031,,,0.1494
22,,,,0.0031,0.0183,,,,,,,,,,0.2043
22.2,,,,,,,,,,,,,,,0.0031
23,,,,0.0031,0.0061,,,,,,,,,,0.1921
23.2,,,,,,,,,,,,,,,0.0061
24,,,,0.0061,0.0031,,,,,,,,,,0.1159
24.2,,,,,,,,,,,,,,,0.0061
25,,,,,,,,,,,,,,,0.061
26,,,,,,,,,,,,,,,0.0427
27,,,,,,,,,,,,,,,0.0091
28,,,0.0366,,,,,,,,,,,,0.0061
28.2,,,0.0092,,,,,,,,,,,,
29,,,0.2287,,,,,,,,,,,,
29.2,,,0.0031,,,,,,,,,,,,
30,,,0.3476,,,,,,,,,,,,
31,,,0.1006,,,,,,,,,,,,
31.2,,,0.064,,,,,,,,,,,,
32,,,0.0183,,,,,,,,,,,,
32.2,,,0.1494,,,,,,,,,,,,
33.2,,,0.0396,,,,,,,,,,,,
34.2,,,0.0031,,,,,,,,,,,,
ET,0.5637,0.071,0.5047,0.805,0.189,0.9677,0.7927,0.0623,0.0687,0.386,0.143,0.8173,0.5177,0.4243,0.3837
H-obs,0.6951,0.7805,0.7988,0.8476,0.8659,0.8293,0.7683,0.689,0.7805,0.7134,0.7683,0.8232,0.8415,0.6402,0.8659
H-exp,0.7301,0.682,0.787,0.8436,0.9042,0.8064,0.8022,0.7678,0.7647,0.7317,0.7817,0.7847,0.8407,0.6424,0.8629
PD,0.876,0.8273,0.9209,0.9521,0.9782,0.9331,0.9328,0.911,0.8993,0.8842,0.9195,0.9112,0.9506,0.8099,0.9629
PIC,0.6753,0.6295,0.7596,0.8252,0.8632,0.778,0.7745,0.7317,0.729,0.6892,0.7501,0.7503,0.8462,0.5738,0.8481
MEC,0.4788,0.43,0.598,0.6903,0.8063,0.6148,0.6124,0.5528,0.5503,0.5006,0.5797,0.5745,0.6778,0.3714,0.724

2018-04-17

遺伝統計学と疾患ゲノムデータ解析

書いた。

遺伝統計学と疾患ゲノムデータ解析

COI:著者謹呈。編集者とはズブズブの関係。

 

序文を引用すると、

学問の学びの最初の一歩は,背景となる基礎理論の丁寧な理解から始まる。第1章では,遺伝統計学の基礎理論について,日本を代表する専門家の先生方にご執筆いただいた。

ゲノム配列解読技術の発達により,大容量かつ多彩なゲノムデータが得られるようになった。エピゲノム・メタゲノムといった多彩なオミクス情報の展開も著しい。これらの最新技術の原理を理解し,その特性を活かしたデータ解析を行うことが重要である。第2章では,最新のデータ解析技術を駆使して活躍されている若手研究者にご執筆いただいた。

ゲノム情報の適切な社会実装を実現するためには,コホート研究,データシェアリング,個別化医療,ゲノム創薬など,ゲノム情報と社会との接点に寄り添った活動が重要である。第3章では,本邦における社会実装の最前線での取り組みをご紹介いただいた。

一般論として,既存のデータ解析手法を適用することより,新たなデータ解析手法を生み出すほうが難しい。多くの研究者に認知され広く使われる手法を生み出すことは,なおさら困難である。しかし国際的なイニシアティブを獲得するためには,新規解析手法の構築においても存在感を発揮することが必須である。第4章では,ゲノム・エピゲノム解析とがんゲノム解析に分けて,独自の解析手法の開発者を対象にご執筆いただいた。

ゲノムデータ解析の現場の声を伝える目的で,新しい試みとして,本章とは別にコラムを設けた。

とある。この前にも、

近年では,疾患ゲノムデータを多彩な生物学医学データと分野横断的に統合し,疾患病態の解明や新規ゲノム創薬,個別化医療へと活用する際に有用な学問としても注目を集めている。しかし,遺伝統計学を担う人材は本邦では特に少なく,学問としての知名度も十分とは言い難い状況にある。

遺伝統計をやっている、と言っている人は本当に少なく、それ故何を学べばいいのか本当に分かりづらい、という状況である。なのでこの本は遺伝統計学に精通した専門家を集めて執筆されているのだが、専門家たちが自身の研究・業績について執筆したため、相当な知識がないと読みこなせない、というのがネック。

 

個人的には、遺伝現象から離れて、バイオインフォマティクスやアルゴリズム屋が1章のデータの捉え方、や、4章のデータ解析手法のところ、を読むと、結構面白いと感じるのではないか、と思う。

2017-11-20

MikuHatsune2017-11-20

logit の代わりにtanh を使ってはいけないのですか?

講義の最中にこんな質問があって、「logit は0-1 の範囲にあって、確率として扱いやすいから」という回答だったが、あとで

スケールが変わっただけで本質的には同じということがアナウンスされていた。

 

logit は

f(x)=¥frac{e^x}{1+e^x}

tanh は

g(x)=¥frac{e^x-e^{-x}}{e^x+e^{-x}}

となり、これらの関係は

g(x)=2f(2x)-1

となる。f(x)=g(x) となるのは、x_0=¥log(1+¥sqrt{2}) のときである。

 

logit が確率として頻用されるのは[0, 1] で収まるからだし、tanh が[-1, 1] で収まるのはdeep learning の活性化関数で頻用される。

また、こういう関数たちは連続なので微分が可能だから扱いやすいのだろう、と数学素人的に思う。

x <- seq(-5, 5, length=3000)
y <- cbind(logistic=gtools::inv.logit(x), tanh=tanh(x))
x0 <- log(1+sqrt(2))
cols <- c("blue", "red")
matplot(x, y, type="l", col=cols, lwd=5, lty=1)
legend("bottomright", legend=colnames(y), col=cols, lwd=5, cex=2)
points(x0, tanh(x0), pch=16, cex=2)

f:id:MikuHatsune:20171120191439p:image

2017-11-18

awk を使って行の要素がすべて0 の行を除外する

RNAseqをしているが、行に5つくらいのデータがあって、それらすべてが0 の行を取り除きたいのだが、いい方法はないかと聞かれた。

列和であれば、$2 とかしてすぐに結果が出せるのはぐぐってよく出る。

しかし、行和をどうこうするのはなかなか出てこない。

 

というわけで、行がいくつあるかの変数 NF と、行についてfor を回して0 かどうかの論理判断をして、すべて0 だった場合は除外、そうでなければprint して残す、というif 文でゴリ押しする。

 

R のワンライナーで適当にinput.txt を作成する

R -q -e 'nr<-8;nc<-5;a<-matrix(runif(nr*nc,1,10),nr);a[3,]<-0;a[sample(seq(a), 6)]<-0;write.table(a,"input.txt",quote=F,col.names=F)'
cat input.txt
1 7.884 0 7.431 7.233 7.049
2 9.302 5.414 8.19 2.165 1.838
3 0 0 0 0 0
4 6.113 8.062 3.974 6.062 6.251
5 7.478 5.364 0 2.168 3.363
6 6.258 3.689 8.887 3.095 1.854
7 3.396 1.824 0 1.51 1.189
8 3.728 1.779 0 0 8.062

ここで、1列目には遺伝子名がはいっているらしく、行がすべて0 なのは3行目、あとはダミーで適当に0 のセルを作っている。

 

ということで、for は2列目から参照して、最後の列 NF まで回す。

その最中に、sum という箱を用意しておいて、真偽判定で真なら0、偽なら1 が入るように加算して、最終的にsum が0 かどうかを判定してprint するかどうか決める。

awk '{sum=0; for(i=2;i<=NF;i++){sum+=($i!=0)}; if(sum!=0) print $0}' input.txt
1 7.884 0 7.431 7.233 7.049
2 9.302 5.414 8.19 2.165 1.838
4 6.113 8.062 3.974 6.062 6.251
5 7.478 5.364 0 2.168 3.363
6 6.258 3.689 8.887 3.095 1.854
7 3.396 1.824 0 1.51 1.189
8 3.728 1.779 0 0 8.062

3行目だけ削除できた。