生命情報解析第4回シグナル配列の統計解析(3) 慶應義塾大学先端生命科学研究所確率分布と有意性 (1) サイコロを10回ふったときの6の目が出る回数とその確率との関係 0.35 0.30 P = 0.00243815649926 0.25 確率 0.20 0.15 0.10 0.05 0.00 0 1 2 3 4 5 6 7 8 9 10 6の目が出る回数棒グラフの右側部分の面積の合計が確率、すなわち有意性を表す確率分布と有意性 (2) 確率有意性確率変数が取る実数検定対象の値 • • • • 確率分布をはっきりさせる検定対象の値から右側の面積を求める “こんなにも大きな値”が出る確率が求まる有意性の指標として使う Z Scoreの特徴 0 Z Scoreの分布 1.96 1.20E-01 • 平均が0、分散(データの散らばり)が1になる • 元の分布が正規分布なら、そのZ Score は標準正規分布となる • どんな正規分布でも、Z Scoreに直せば同じ土俵(?)で確率計算ができる • Z Scoreが1.96を超える確率は0.025 1.00E-01 確率 8.00E-02 6.00E-02 4.00E-02 2.00E-02 0.00E+00 -4.47 -3.13 -1.79 -0.45 0.89 2.24 3.58 4.92 6.26 7.60 8.94 10.29 11.63 12.97 14.31 15.65 16.99 18.34 19.68 21.02 22.36 Z Score Z Scoreの計算サイコロを100回振って、90回”6”の目が出るときのZ Scoreは N obs − Np 90 − 100 ⋅1 / 6 Z Score = ≅ 21.17 = Np (1 − p ) 90 ⋅1 / 6 ⋅ 5 / 6 塩基の方も… • ここでは簡単のため、１塩基の偏りだけを考える • ゲノム全体の塩基組成を考えて、塩基iが対象となる場所において観測される確率はpiとする • 今、N本の配列のうち、 Ni個について、対象となる位置に塩基iが観測された • この条件では通常、Niは正規分布に従う頻出塩基の統計的有意性 Z Score = N Ni Bi N i − NBi NBi (1 − Bi ) :解析する配列数 :観測された塩基iの数 :ゲノム中における塩基iの割合 Z Scoreは標準正規分布に従う Z Score > 1.96なら、P < 0.05 複数の塩基の有意性を同時に検定するには? • ゲノム全体の塩基組成をBa=0.3, Bc=0.3, Bg=0.3, Bt=0.1として、与えられた位置における塩基iの個数Niが Na = 40, Nc =40, Ng = 10, Nt = 10のとき、偏りは有意か？ • Z Scoreを４つも計算すると… – ４つも値が出て、取り扱いが煩雑になる – 偶然に高い値を示すものが出やすくなる • Χ2値を使う Χ 2値 • n個の互いに独立なZ Score：Z1, Z2, Z3, …, Zn があるとき、 Χn2値 = Z12+Z22+Z32+…+Zn2 • Χn2値は自由度nのΧn2分布に従う Χ 2分布 1.4 1.2 1 自由度 1 2 3 4 5 確率 0.8 0.6 0.4 0.2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 2 Χ 値 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 Χ 2分布に従う値を求める • しかし、Na, Nc, Ng, Ntは互いに独立ではないため、これらのZ Scoreを足しても自由度４のΧ2分布には従わない。そこで… 2 P B ( − ) i χ 32値 = N ∑ i Bi i = a,c,g, t 但しNは解析対象の配列数 • 塩基がそれぞれBiの頻度で出現するとき、上記Χ2 値は自由度３のΧ2分布に従う • Χ2値>12.84ならP < 0.005 Χ 2計算の例 • ゲノム全体の塩基組成をBa=0.3, Bc=0.3, Bg=0.3, Bt=0.1として、与えられた位置における塩基iの個数Niが Na = 40, Nc =40, Ng = 10, Nt = 10のとき、偏りは有意か？ (0.4 − 0.3) 2 (0.4 − 0.3) 2 (0.1 − 0.3) 2 (0.1 − 0.1) 2 + + + 0.3 0. 3 0.3 0.1 ≅ 0.033 + 0.033 + 0.133 + 0 ≅ 0.2 • この例では偏りが有意とは言えない。 14000 大腸菌開始コドン周辺の塩基のχ2乗値 12000 Χ 自乗値 10000 8000 6000 4000 2000 0 -100 -50 0 開始コドンからの相対位置 50 100 ２つの数式の関係(1) 2 P B ( − ) i χ n −12 値 = N ∑ i Bi i =1 n χn 値 = Z + Z + Z + L + Z 2 2 1 2 2 2 3 • 4塩基を2種類に分類して考える – プリン(A,G)、ピリミジン(C,T) • ゲノム中のプリン、ピリミジンの頻度をそれぞれBpur, Bpyrとする • 対象となる位置で観測されたプリン、ピリミジンの頻度をそれぞれPpur, Ppyrとする。但し、Ppur + Ppyr = 1 2 n ２つの数式の関係(2) ⎛ ( Ppur − Bpur ) 2 ( Ppyr − Bpyr ) 2 ⎞ ( NPpur − NBpur ) 2 2 ⎟= N⎜ = Z pur + ⎟ NB (1 − B ) ⎜ B B pur pyr pur pur ⎠ ⎝ • ２種類の塩基の数をもとに計算したΧ2値は自由度1 のΧ2分布に従う • 4種類の塩基の数をもとに計算したΧ2値は自由度3 のΧ2分布に従う • 自由度は自由に動ける変数の数を意味する演習問題 Ba=0.3,Bc=0.2,Bg=0.2,Bt=0.3として、与えられた位置における塩基iの個数 Niが (1) Na = 50, Nc = 30, Ng = 10, Nt = 10 (2) Na = 500, Nc = 300, Ng = 100, Nt = 100 のときの増加情報量、χ2値を求めよ。 log23 ≒1.585、log25 ≒2.322 演習問題解答 • (1), (2)ともに • (1) 1 5 3 3 1 1 1 1 増加情報量 = log 2 + log 2 + log 2 + log 2 2 3 10 2 10 2 10 3 ≅ 0.368 + 0.175 − 0.1 − 0.159 ≅ 0.285 ⎛ ( 12 − 103 ) 2 ( 103 − 15 ) 2 ( 101 − 15 ) 2 ( 101 − 103 ) 2 ⎞ ⎟ + + + χ = (50 + 30 + 10 + 10)⎜⎜ ⎟ 3 3 1 1 5 5 10 ⎝ 10 ⎠ ≅ 100 × (0.13 + 0.05 + 0.05 + 0.133) ≅ 36.7 2 • (2) 3 2 3 3 2 ⎞ 1 1 2 1 1 2 1 ⎛ ( − ) ( − ) ( − ) ( − ) χ 2 = (500 + 300 + 100 + 100)⎜⎜ 2 310 + 10 1 5 + 10 1 5 + 10 3 10 ⎟⎟ 5 5 10 ⎝ 10 ⎠ ≅ 1000 × (0.13 + 0.05 + 0.05 + 0.133) ≅ 367 特定のシグナル配列の存在頻度 • 様々な塩基配列の偏りを調べるのではなく、特定のシグナル配列の存在頻度を調べたい (ex. SD配列 “AGG”) • 最も単純なのは、頻度＝あるシグナル配列が観測される配列数÷解析対象の配列数大腸菌開始コドン周辺の”AGG”の頻度 0.14 0.12 "AGG"の頻度 0.1 0.08 0.06 0.04 0.02 0 -100 -50 0 開始コドンからの相対位置 50 100 頻出塩基配列パターンの統計的有意性 Z Score = N obs − Np Np (1 − p ) Nobs :パターン観測数 N:解析する配列数 p:ゲノム中におけるパターンの割合パターンの出現頻度がpのとき、Z Scoreは標準正規分布に従う Z Score > 1.96なら、P < 0.05 大腸菌開始コドン周辺の”AGG”のZ-Score 80 70 60 Z Score 50 40 30 20 10 0 -100 -50 0 -10 開始コドンからの相対位置 50 100 翻訳開始シグナル抽出結果 Escherichia coli 16S rRNA 3- terminal: gcggttggatcacctcctta3 Expected SD Sequence: 5taaggaggtgatccaaccgc Pat. agga ggag aagg gagg gaga Z-Sc. 94.97 82.94 58.15 53.08 42.23 Pos. -11 -10 -11 -11 -9 シグナル配列出現の評価 • 塩基の偏り – 偏りの程度 … エントロピー、増加情報量 – 偏りの有意性 … Χ2値 • 配列パターン – 出現の程度 … 頻度 – 出現の有意性 … Z Score 分子レベルの生命現象の根幹～セントラルドグマ～ ATG TAA DNA 転写 RNA UAA AUG 翻訳タンパク質機能 RNAレベルで機能する分子 • tRNA • rRNA • Other non-coding RNA • Translational regulation by mRNA tRNA tRNA ACC ACGAGUACA UGCUCAUGUUGG rRNA リボソーム Methionine fMet-tRNAf 16S rRNA AUUCCUCC mRNA AUG AGGAGG 開始コドン Shine-Dalgarno sequence 16S rRNAの3‘末端はShine-Dalgarno配列と対合する翻訳での遺伝子の発現制御 Fe Ferritin gene 5’ 3’ 二次構造による終止コドンの読み飛ばし二次構造 mRNA UAA 通常の長さのタンパク質 Steneberg, P. 2001 リードスルーによってできた長いタンパク質 Function of readthrough product is stronger 1034 hdc AUG 2981 4274 UAA UAA hdc gene is expressed in tracheoles in larvae of D. melanogaster Possibility of Regulation by readthrough? Long product Short product Branching of lumens are inhibited strongly Branching of lumens are inhibited weakly Steneberg and Samakovlis, 2001 cDNA配列を用いた転写産物の収集ＡＴＧ遺伝子 TAA ＤＮＡ転写 mRNA ＡUG UAA 逆転写 cDNA ＡTG UAA マウスcDNA配列の網羅的収集コード領域を持たないcDNA? ゲノム cDNA Numata et al. 2003 さらにコード領域を持たない多数のcDNA? タンパク質をコードしないcDNA配列が多くある cDNA ゲノムゲノムの62.5%をカバー多くのRNAは翻訳されなくても機能を持つ？非翻訳RNAが多量に存在？多数の非翻訳RNAの存在が予想されているものの、ほとんどは機能未知 RNAの二次構造予測 • 一本鎖RNAはDNAに比べ、自由な構造を取ることが可能 • RNAが機能する上で立体構造が重要になってくる • 二次構造は、どの塩基とどの塩基が結合しているかを表す • 一次配列から二次構造を予測しよう！ tRNAの二次構造予測の例 GenBank tRNA配列 http://www.genome.ad.jp/dbget-bin/www_bget?gb:ECOCPTGG http://www.genome.ad.jp/dbget-bin/www_bget?gb:ECOPHER Zukerのmfold http://www.bioinfo.rpi.edu/applications/mfold/old/rna/form1.cgi