意訳：Finding Hidden Messages in DNA (Bioinformatics I): Week1
本日はバイオインフォマティクスの知識を深めるために、CourseraのFinding Hidden Messages in DNA (Bioinformatics I)を要約します。
見た動画はweek1の２つ。
https://class.coursera.org/hiddenmessages-007/lecture/9
https://class.coursera.org/hiddenmessages-007/lecture/11
どちらの動画もDNAの複製に関することで、問題設定としては、
Input: ゲノム配列の中の複製が関わる箇所(500塩基), Output: DNA複製開始位置(数~十数塩基)
言い換えると、環状のDNAからDNAの複製が開始される位置を把握するにはどうしたらいいかというもの(ただしゲノム配列の中の複製が関わる区間がどのあたりなのかはあらかじめわかっていると仮定する)。これに対しては我々エンジニアはなるほど、じゃあアルゴリズムを考えようかと思いがちだがこれでは問題が曖昧すぎて実際お手上げ状態。一方バイオロジストであればDNAを切っていき、DNAの複製が行われなくなることを確認すればDNAの複製が開始された位置を把握できる。エンジニアがバイオの分野に参入するには適切な問題設定や修正ということが、私達が普段直面する問題よりも遥かに重要なのかも。ではどうするかというと、そのゲノム配列には必ず複製開始位置を知らせる隠されたメッセージが含まれているので、そのメッセージを探すようにしようというもの。問題設定を以下のようにしてみる。
Input: 複製に関わるゲノム配列の一部, Output: 頻繁に出現する、長さkの文字列を抜き出す(通称k-mer)
これでエンジニアもこの問題を対処できるようになりました。エンジニア的にはこれで満足だけど、これはバイオの問題なので必ずバイオロジストの視点からこの問題設定が適切かどうか確認する必要がある。これも当たり前ではありますが、非常に重要なことっぽい。で、バイオロジストに伺ったところ、この問題設定で大丈夫とのこと。実際、DNAの複製はDNAポリメラーゼによって行われるが、その複製を開始するのはDnaAタンパクであり、DnaA boxと呼ばれる9塩基程度の短い箇所へDnaAは結合するらしい。上の問題設定で得られる解から私たちはDnaA boxの候補を見つけ出すことができたということになる。
ここで一つ注意点として、DNAは２対の鎖状になっており逆相補鎖(A-T, C-Gに対応)が存在するので、例えばATCとTAGを同列のものとして扱うこともできる。なので先程の問題のOutputが複数ある場合にその候補をさらに絞ることが出来る。
上記のようなOutputの場合他の候補よりも逆相補鎖が存在する候補のほうが可能性が高い。これで全て解決かというとそうではない。今まではめちゃ長いゲノム配列の一部のoriC(オリジンとも言われる)という、複製開始の区間の中からDnaA box候補を探す、ということをしていた。先程まではそのoriCの位置が既知であると仮定していた。なのでそもそもoriCどこよという話。ではどうするかというと、ゲノム配列を多数の窓に区分けし、先ほどのDnaA box検索をそれぞれの窓に対して行う。多数の頻出文字列が存在する窓はoriCである可能性が高いというわけだ。その窓の大きさをL、文字列の出現回数をtとする。その群れのことを(L-t)-clumpと呼ぶ。
Input: ゲノム配列, Output: (L, t)-clumpを形成する全てのk-mers
それでは実際のゲノムに対してこのアルゴリズムを適用するとどうなるか。9-mers、(500, 3)-clumpが1904個ほどヒット。これは何を意味するかというと、窓枠500の中で3回出現する塩基9個の文字列が1904個存在するということ。つらい笑
さてどうするか。次回へ続く。
DNAポリメラーゼの働き
バイオの知識を深めるために、これから毎日PDB-101: Home Pageの記事を和訳します。
本日はDNAポリメラーゼについて。以下、意訳。
生命の神秘
DNAポリメラーゼは生物的なプロセスにおいて重要な役割を果たしています。それは遺伝情報の複製を行います。細胞が分割される度、DNAポリメラーゼはDNAを複製し、親の細胞から子の細胞へと情報を受け渡し世代を超えて情報が伝えられていきます。このDNAの引き継ぎは地球に初めて誕生した頃の生命と私達の細胞が生物的につながりがあることを示しています。その点から考えると、私達のDNAはとても貴重なものであると言えるかもしれません。
驚くべき正確さ
DNAポリメラーゼは最も正確な働きをする酵素です。それは毎回正確にDNAのコピーを作成しつつも、10億に1回ミスするかしないかの精度を誇っています。これは私達の情報処理能力と比べても優れていることがわかります。なぜならそれは1000冊の小説を読んで、その中からたった一文字の誤字を見つけるようなものだからです。C-G、A-Tをマッチさせる(シトシンはグアミン、アデニンはチミンと対応づけられる)その素晴らしい能力は、それが高い正確性を兼ね備えているという特殊性を実現しています。しかしDNAポリメラーゼはそれだけではなく、DNAの複製後必ずその複製物を校正し必要とあらば誤った箇所を切り捨てることまでしてくれます。
囚人と血統
DNAはユニークですので、他人との識別に有効です。もし血痕の一つでもあれば、その血が誰のものなのか調べることができます。もちろん、血痕の中に入っているDNAの量はごくわずかのため、それを増やすためにDNAポリメラーゼが使われます。そのわずかなDNAはPCR(Polymerase Chain Reaction: DNAポリメラーゼの複製機能を利用してDNAを増やす手法)により増やされ、簡単に解析できるようになります。PCRのプロセスは以下のようです。まず、そのわずかなDNAは試験管にDNAポリメラーゼとともに入れられます。その試験管を温めると、DNAの2本鎖は2本に解かれます。DNAポリメラーゼは2本の解かれたそれぞれのDNAから、新たなDNA鎖を作成します。この新たに作成された2つのコピーDNAは同じように温められ、複製され、4本のコピーDNAになります。この反応を繰り返すことで大量かつ同一のDNAが複製されます。私達の体にある、通常のDNAポリメラーゼの場合、温めた段階で壊れて機能を失ってしまいます。しかしThermus aquaticusという、温泉に生息する細菌が持つDNAポリメラーゼの場合高温にも耐えることができるのですべてのPCRでそのDNAポリメラーゼが活用されています。ちなみに、これがPDBに登録されているその酵素のデータです。
構造
このシンプルな酵素は手のような形をしています。上記二つとも、細菌のポリメラーゼです。左がEscherichia coliで、右がThermus aquaticusのポリメラーゼです。左のE.coliポリメラーゼの下部に示された緑の領域はPDBに登録されておらず、上部がよく研究されています。面白いことに、親指とその他の指との間にはDNA鎖がちょうどフィットする間隔があり、実際その手のひらのような箇所にDNAがフィットします。写真で紫と緑に色が付いている、手の中央あたりでしょうか。金型となるDNAは紫、新たに複製されたDNAが緑色に色づけられています。この酵素は３つの活性な箇所があり、Synthesis(新規DNAが収まっている箇所), Proofreading(新規DNAを校正する箇所), そしてPrimer removal(DNA複製にRNAの断片を除去する箇所)と示されている箇所です。E.coliはProofreading機能がある一方、Thermus aquaticusはその機能がありません。
DNAポリメラーゼ
すべての生物はDNAポリメラーゼを持っています。とてもシンプルなものから全てを行うようなものまで。私達の体の中のポリメラーゼはより複雑な構造をしています。指輪のような形のものやαヘリックスを解いたようなものものなどの複数のタンパク質から成っています。
上の写真は左上がE.coli、右上がヒト、下のがウイルスのDNAポリメラーゼです。それぞれ形も大きさも異なりますが、DNAの包み方や合成反応などは似通っています。
Scalaで自己共分散と自己相関の実装
自己共分散も自己相関も時系列データ元々の分布（関数）の時間軸をずらした時にどれほど元のグラフと相関があるかについて計る指標。自己相関は自己共分散の値を使って求め、自己共分散と違い値の大きさが単位に依存しない（分散と標準偏差みたいな関係。本質的には同じだけど標準化されているかそうでないかの違い）。
自己共分散＝（（元の関数の各値ー平均）×（ずらした関数の各値ー平均））の平均。
def autocovariance(raw:Seq[Double],lag:Int):Double={
lazy val mean=meanf(raw)
lazy val lag_raw=raw.drop(lag)
if(raw.length<=lag) 0
else lag_raw.zip(raw).map{x=>(x._1-mean)*(x._2-mean)}.reduce{(a,b)=>a+b}/raw.length
}
実装ではまず初めに引数のデータの平均と、引数のlagの分だけ先頭から要素を削除したSeqを各変数に代入している。その変数がlazy（遅延評価）になっているのはその後の行を見たらわかるが、もしlagがデータ数以上だった場合自己共分散は常に0になり計算がもったいないから。んでもしlagが必要以上に大きくなかったら、時間軸をずらしたラグデータと元のデータをzipする。次にmapメソッドで対となる要素ー平均を掛け合わせたものの合計/元のデータ数をする。なぜここで同じ平均を使っているのかというと確かに自己共分散を求めるにはそれぞれのデータの平均、上の式で言えばµt,µsを使う必要があるが本来ラグデータは時間軸をずらしただけのデータであり要素数も平均も元のデータと同じ。実装の都合上やりやすいからあたかも要素を削除しているように見えるけどそれは単純に見かけの話。
次に自己相関。自己相関は自己共分散の値を利用して求めるだけ。
ρ(h)＝ r(h) / r(0)
自己相関（lag h）＝自己共分散（h）/ 自己共分散（0）
普通のRとかでは自己相関関数（acf）を実行すると20くらいまでそれぞれのlag hに応じた自己相関を作ってくれるから、それも真似てみる。
def autocorrelation(raw:Seq[Double])
:IndexedSeq[Double]={
val r0=autocovariance(raw,0)
if (raw.length<20) (0 to raw.length).map(x=>autocovariance(raw,x)/r0)
else (0 to 20).map(x=>autocovariance(raw,x)/r0)
}
ScalaTestすごい楽。
試してみたのは単純な、FunSuiteクラスを継承するタイプ。
class Datates extends FunSuite {
test("name of this test"){
assert(0==0,"if fail,message appears")
}
でもDoubleが比較に入ってくると小数点以下がちょこっと違うだけでほとんど同じ！みたいなケースが出てきて、そこは融通きかしておっけーにしちゃえよ！とかあるからそんな時はEquality[Double]クラスに対して暗黙の型変換を行い、どれくらいまで許容するかを指定してあげる。
import org.scalautils.TolerantNumerics._
implicit val doubleEquality = tolerantDoubleEquality(0.001)
これでだいたい小数点三桁以下までなら許容するはず。それと、比較する時に上では==を使っていたけど許容範囲を指定する時は===を使う。
assert(2.0009===2.0)
スピアマンの順位相関係数を実装 in Scala
この前ピアソンの相関係数を実装したので次はノンパラメトリックなスピアマンの順位相関係数を実装してみた。ピアソンの場合、外れ値の影響を大きく受けるのと、厳密にはデータが正規分布であると仮定している。つまりデータ数が小さい時は正確性に欠ける。翻ってスピアマンの方はノンパラメトリックなのでデータの分布は仮定せず、外れ値の影響を受けづらい。
- スピアマンの計算手順
1. ２つの変数の各要素に順位をつける（大きい順でも小さい順でも）
2. 対となる要素の順位の差を求める
3. それを基に相関係数を求める
4. 分析表と照らし合わせて有意かどうか判断（ここは実装してない）
相関係数の式。diは要素の順位差、nはデータ数。
単純な例を考えてみる。一日の睡眠時間と集中力が持続する時間の平均値 が相関を持っているのか考える。
- 一日の睡眠時間
8 , 4 ,5
- 集中力の時間平均
60 , 30 , 25
まずはそれぞれの変数の各要素に順位づけをする。
(8 , 3) , (4, 1) , (5, 2)
(60 , 3) , (30 , 2) , (25 , 1)
次にそれぞれの順位の差の２乗を求める
(3-3)^2 , (1-2)^2 , (2-1)^2
=0 , 1 , 1
後は相関係数の式に当てはめる。
r = (6 * 2) / (3^3-3)
= 12 / 24
= 0.5
要はこれをコードにするだけ。
def labeling(data:Stream[Double])={val placeI=data.sorted.zipWithIndex.map(a=>(a._1,a._2+1)).toMapdata.map(a=>placeI(a))}def rankzip(data1:Stream[Int],data2:Stream[Int]):Stream[(Int,Int)]=data1.zip(data2)//2つのデータをソートしてzipdef difsqured(zipped:Stream[(Int,Int)]):Stream[Double]=zipped.map(x=>math.pow(x._1-x._2,2))//対となる値の差をそれぞれ求めて２乗するdef spearman(dsdata:Stream[Double])=1-(6*dsdata.reduce{(a,b)=>a+b})/(math.pow(dsdata.length,3)-dsdata.length)//スピアマンの順位相関係数。同じ順位の場合は昇順にしているため若干本来と違う。a:4.5,b:4.5 => a:4,b:5 にしてる。
labelingメソッドが何とも効率悪そうだなあ
sbtって別にシンプルじゃないよな。って思うのは他のビルドツールを知らないからなのだろうか。
sbtについてどしっとまとめたことがなかったからちょいとメモる。sbtのダウンロードはHomebrewを使ってちょちょいと出来る。brew install sbt、てな感じで。sbtのプロジェクトを作るにはまずルートディレクトリでsbtコマンドを実行。するとprojectフォルダなどができる。それに加えてbuild.sbtファイルをルートディレクトリ下に作成し名前やバージョン、依存性（ライブラリ）等を書き込むことによってプロジェクトを管理する。以下がルートディレクトリ下のフォルダ構成（srcに対応したbinフォルダも出来上がるし成果物はtargetフォルダ内に入る）
src/main/resources/<リソースファイル>scala/<scalaソースコード>java/<javaソースコード>test/resources<テストリソース>scala/<scalaテストコード>java/<javaテストコード>project<sbtプロジェクトファイル>target<コンパイルされたクラスファイルなどsbtが自動生成したファイルの出力先>
よく使うsbtコマンド。
~compile : 保存するたびに自動コンパイルしてくれる。ごっつ便利。
clean : 生成されたファイルをtargetフォルダごと削除
reload : ビルド定義を再読み込みする(ルートディレクトリにある.sbtファイルとprojectディレクトリ下にある.sbtと.scalaファイルが対象)
run : ラン。ラン。ら〜ん。
update : 依存関係を変更したらしましょ。
test : てすてす。
package : jarファイルを作ってくれる。何これうれしい。
build.sbtとBuild.scala、どちらを使えば良いのか＝＞基本的にbuild.sbtで十分。sbtをカスタマイズしたい時はBuild.scalaの方が良い。らしい。build.scalaを書く上で大事なことがいくつか。まず、設定ではキーと値がペア。左がキー。右が値。
name := "i~am~hungry~~~~~project"
次に、それぞれの設定を区切るのは空行。<=これまじでどうにかならんのかな。なんか嫌だ。
version := "1.0" scalaVersion := "2.9.1"
ライブラリの追加。
libraryDependencies += groupID % artifactID % revision もし複数同時に追加したいのなら+=ではなく++=を使用する。sbtでは空行が大きな役割を果たすので無駄に空行は入れないようにする。たぶん。revisionはversionみたいなもの。
libraryDependencies ++= Seq(
groupID % artifactID % revision,
groupID % otherID % otherRevision
)
テストには必要だけどメインの方には必要ないみたいなライブラリはrevisionの後にに%"test"を追加する
groupID % artifactID % revision % "test"
基本的にこれで十分なんだけど全てのパッケージが１つのサーバーに置いているとは限らないのその場合はresolverとやらを追加する必要がある。Eclipseとsbt同時に使いたい場合、projectディレクトリ下にplugins.sbtを作成し、
addSbtPlugin("com.typesafe.sbteclipse" % "sbteclipse-plugin" % "2.4.0")
と書き込むことによってsbtプロジェクトをEclipseでも使うことが出来る。それでも依存性を変更した後とかはEclipseではエラーだけどsbtではエラーじゃないみたいな状況になったりして、そういう時は
clean
reload
update
