階層ベイズモデルを用いてバトルMCの強さランキングを作成する

誰が最強のMCなのか。MCバトル好きの間では常に議論になることです。

MCバトルの勝敗は、各MCの実力だけでなく、MCのコンディションやMC同士の相性、バトルの審査方法、会場の空気など様々な要因によって決まります。実力あるMCであっても意外な人物に惨敗したり、成績が振るわない対戦形式があったりして、勝敗の結果からMCの実力を推定するのは案外難しいものです。

勝敗の結果から各プレイヤーの実力を推定する方法として階層ベイズモデリングがあります。例えば将棋の勝敗データに適用したブログ記事→階層ベイズモデルで勝敗データからプロ棋士の強さを推定する - StatModeling Memorandum。今回は、このブログの著者が書いたStan・Rベイズモデリング本（通称アヒル本）を参考にして、階層ベイズモデルを用いてバトルMCの実力を推定してみました。

勝敗データ

ネット上に散らばっている勝敗データをかき集めました。

2000年代のバトル結果はネットであまり集められなかったので、2010年～2017年のデータのみ対象にしました。また、地方で開催されているバトルを含めるとキリが無くなるので、UMB（Ultimate MC Battle）やKOK（King of Kings）の本戦などの全国大会だけを選びました。

対象にした大会は、UMB2010～2017、KOK2015～2017、罵倒Grand Championship 2014・2015・2017、戦極MC Battle 第6・8・9・10・13・16章、フリースタイルダンジョンの初代モンスター時代の1on1のみ*1、フリースタイルダンジョンの東西口迫歌合戦です（一部抜け落ちている試合あり）。全部で661試合。もし他に勝敗データを提供してくださる方がいらっしゃれば、ぜひコメントいただけると嬉しいです。

負けてしかいないMCはほとんど推定に役立たないので除外したところ、139MC・477試合が残りました。この477試合のデータを用いて、ベイズモデリングをおこないました。

こんな感じでExcelにまとめました。

モデル式

アヒル本の将棋の棋士の強さ推定モデル（p.188~190）をほぼそのまま拝借します。

Performance [g, 1] ~ Normal ( μ [Loser [g] ], σ [Loser [g] ] )

Performance [g, 2] ~  Normal ( μ [Winner [g] ], σ [Winner [g] ] )

Performance [g, 1] < Performance [g, 2] 

μ [ n ] ~ Normal ( 0, σn )

σ [ n ] ~ Gamma ( 10, 10 )

gは試合のインデックスを示し、今回は477まであります。nはMCのインデックスを示し、1~139まであります。Performance[g,1]とPerformance[g,2]は試合gにおける勝者と敗者のパフォーマンスです。各試合におけるプレイヤーmのパフォーマンスは、平均μm・σmの正規分布から生成されると考え、1・2式目のように表現しています。σは各MCの"勝負ムラ"を示し、弱情報事前分布であるガンマ分布に従うと設定しています。パフォーマンスの大きい方が勝者となる、という制約が3式目です。

詳しいことはアヒル本を読んでください。

Stan・Rコード

下のようなコードをStanで書きました。というか、アヒル本からコピーしました。

data {

  int G;

  int P;

  int<lower=1, upper=P> LW[G,2];

}
parameters {

  ordered[2] performance[G];

  vector[P] mu;

  real<lower=0> s_mu;

  vector<lower=0>[P] s_pf;

}
model {

 for(g in 1:G) {

    for (i in 1:2) {

       performance [g,i] ~  normal(mu[LW[g,i]], s_pf[LW[g,i]]);

     }

 }

 mu ~ normal (0, s_mu);

 s_pf ~gamma (10,10);

} 

 でRからStanを走らせました。

P <- max(data2[,c("Loser_id","Winner_id")])

dat <- list( G =G, P = P, LW=data2[,c("Loser_id","Winner_id")] ) 

stanmodel <- stan_model(file="model.stan")

fit <- sampling(stanmodel, data=dat, pars=c('mu','s_mu','s_pf'), seed=1234, chains =4 , iter = 2500) 

結果

推定された強さμの高いトップ20のMCは下の通りになりました（μの5~95パーセンタイル値を図示してます）。図の下に行くほど強くなっています。誰もが納得のR-指定がトップ1に輝いていますね。その後のGADORO、mol53、晋平太などの並びも、それなりに妥当ではないでしょうか。

注意すべきは、推定結果は今回扱った勝敗データに依存しているという点。当たり前ですが、扱う勝敗データを拡張すれば結果は変わります。今回は全国レベルの大きな大会だけに絞ったので、地方予選で勝ち続けていてもそれが推定結果に反映されてません。さらに、ここ2~3年の人数が大きい大会で優勝した人は、勝利した試合数が多いので有利になっています。

ちょっと意外な人がトップ20に入っているのもデータ選択の妙ですね。GASHIMAや鎮座DOPENESS、般若は、勝利した試合しか今回のデータに含まれてないのでランキング上位ですが、データ数が少ないので信頼区間が非常に大きくなっています。

次に"勝負ムラ"σの大きいMCと小さいMCのトップ10です。上の図がσの大きいMC、つまりムラがあるMCで、下の図がσの小さいMC、つまり安定しているMCです。やはり図の下側からトップ（あるいはワースト）の順になっています。

サイプレス上野が一番ムラがある、という結果でした。なぜか納得です。戦極8章でR-指定を倒して優勝した一方で、フリースタイルダンジョンなどでの連敗が反映されているのでしょうか。安定しているMCには、やはり実力のあるMC達がランクインしてますね。逆にここに入っていない晋平太や崇勲は、割と負けるけどハマったら強いというタイプでしょうか。

最後に

できれば年ごとの強さを推定したかったのですが、そのためにはデータ数が少なすぎるみたいです。アヒル本のp.231~235を参考に時系列モデルを作成してみましたが、あまり納得できる結果にはなりませんでした。計算時間もかなりかかるし…。残念。

地方予選も含めて膨大なMCバトルのデータを持っている方がいれば、ぜひ提供してくださいw。 

個人的にはR-指定と晋平太の2トップになって欲しかったなぁ。

(追記：2018.08.05)

解析のもとにしたExcelファイルアップロードしました→MC_battle.xlsx 。

だいぶ適当ですが。例えばUMB2010は自分の印象に残っているバトルしか記載してません。

StanとRでベイズ統計モデリング (Wonderful R)

*1:晋平太や呂布カルマなどチーム戦が導入されてからの1 on 1はまだ含めてない。

論文メモ: ヨコエビ Hyalella aztecaのゲノム

 「底質毒性と進化毒性学のモデル生物：Hyalella aztecaのトキシコゲノム」

Poynton HC, Hasenbein S, Benoit JB, Sepulveda MS, Poelchau MF, Hughes DST, Murali SC, Chen S, Glastad KM, Goodisman MAD, Werren J, Vineis JH, Bowen JL, Friedrich M, Jones J, Robertson HM, Feyereisen R, Mechler-Hickson A, Mathers N, Lee CE, Colbourne JK, Biales A, Johnston JS, Wellborn GA, Rosendale AJ, Cridge AG, Munoz-Torres MC, Bain PA, Manny AR, Major KM, Lambert F, Vulpe CD, Tuck P,  Blalock B, Lin YY, Smith ME, Ochoa-Acuña H, Chen MM, Childers CP, Qu J,  Dugan S, Lee SL, Chao H, Dinh H, Han Y, Doddapaneni HV, Worley KC, Muzny DM, Gibbs RA, Richards S, 2018, The Toxicogenome of Hyalella azteca: a model for sediment ecotoxicology and evolutionary toxicology, Environ Sci Technol.

ヨコエビHyalella aztecaのドラフトゲノム。Redundansでアセンブリされたゲノムデータ（Hazt_2.0）やその情報はNCBIのページにあります。全長550.9Mbpで、偽遺伝子含めた遺伝子数は20,022個。

ドラフトゲノムを作ったことでアノテーションできる遺伝子数が増えたよ、という話。ゲノムアセンブリ前は、有害物質曝露によって発現変動する遺伝子（DEGs）のアノテーション率は低く、Cdの場合は13%、Znは29%、PCB126は0%、Cyfluthrinは20%でした。リファレンスゲノムができた効果で10~30%アノテーション率が増加したそうです。それでも50%以下で、結構厳しい。Daphnia pulexのゲノム論文（2011, Science）は確か環境変化に対するDEGsはアノテーションされにくい、系統lineageに特異的な遺伝子が多い、という話でしたが、それとも一致している傾向？

「ヒトゲノムを解読した男：クレイグベンター自伝」感想

ひと月くらい前に読了。

バイタリティに富んだ人間という一言に尽きます。

彼が海軍衛生兵としてベトナムへ従軍していた時の話。腹部に傷を負った二人の男が同時期に病院に運ばれてきた。ひとりは生きて当然だと思われたが、すぐに息を引き取った。もうひとりは腸や肝臓、脾臓の一部が吹き飛ばされてすぐに死ぬと思われたが、何日もしゃべり続けて生きながらえていた。二人の対比から、ベンターは人間の精神と意志の力はどんな薬よりも効果があることを学んだ、とあります。ベンターはこの出来事以来ほとんど毎日のように、今日も生きていたいと願うようになったそうです。

これは彼の自伝を象徴するエピソードと思います。とにかく並外れた情熱と行動力が本書の節々から感じられます。ヒトゲノム計画のあとも、微生物メタゲノム*1や人工生命など新しい研究課題に取り組んでますしね。

ゲノム研究を始める前の学部生のころ（ベトナムからの帰還後）からPNASに論文を載せまくるなど、超優秀だったとは知りませんでした。

ショウジョウバエのゲノムを手に入れた際、その分析を短期間で終わらせるために世界中のショウジョウバエ研究者を一か所に集め、11日間かけて皆で解析して3本の一連の論文を発表したと言います。おかけでゲノムプロジェクトを始めてから論文発表まで1年足らずで済ませています。すごすぎ。こういう時代が動いてるというか、興奮が充満してるような感じ、ぜひ参加してみたい。

あと民間企業でもアカデミックな研究ができるんだなぁと感心しました（小並感）。彼のマネをして上手くいくかは知りませんが。

ヒトゲノムを解読した男 クレイグ・ベンター自伝

*1:クレイグベンターがメタゲノム解析の走りだったのは知りませんでした。サルガッソー海のこれ。

論文のメモ: EDA・TIEについての最近の論文

最近出たEDA（Effect-Directed Analysis）とTIE（Toxicity Identification Evaluation）に関する論文。

「底質の複雑な毒性の診断：TIEとEDA」

Li H., Zhang J., You J., 2017, Diagnosis of complex mixture toxicity in sediments: application of toxicity identification evaluation (TIE) and effect-directed analysis (EDA), Environ Pollution: in press.

総説。様々な化学物質が混在している汚染底質の中で、具体的にどの物質が悪影響の原因になっているのかを探る手法にSediment TIEとEDAがあります。その2つの手法の欠点を補い合い、より効果的に悪影響の原因物質を同定しようという話。ざっと読んでみて、Burgess et al. (2013, ETC) の総説から新しい話はほぼなさそうでした。TIEの利点はbioavailabilityを考慮してる点で、EDAの利点は有機物の細かい分画を実施してる点。それぞれの欠点は逆。↓表にするとこんな感じ。

で、EDAでおこなう有機の分画をTIEでやれば良いよね、その時は強力な有機溶媒で抽出したりせずbioavailabilityを考慮したpassive samplingなどの手法を使おうね、という話の展開。最後、おまけ程度にAOP（Adverse Outcome Pathways）の話もあります。

同じグループから同時期に似たような総説も出てるみたいです↓。 こちらは未読。

You J., Li H., 2017, Improving the accuracy of effect-directed analysis: the role of bioavailability, Environ Sci Processes Impacts. 

「パッシブドージングとin vivo毒性試験を組み合わせた底質EDA」

Qi H., Li H., Wei Y., Mehler W.T., Zeng E.Y., You J., 2017, Effect-directed analysis of toxicants in sediment with combined passive dosing and in vivo toxicity testing, Environ Sci Technol 51:6414-6421.

上の総説で引用されていたもの。総説で言うところのbioavailabilityを考慮した新しいEDAとはこの論文の内容を指すっぽいので、 実際の研究の流れを知りたければ総説よりこちらを読むべし。

河川底質を有機溶媒で抽出し、その抽出物をゲル浸透クロマト等で分画して毒性試験に使ってます。ここまでは伝統的なEDAですが、抽出物の曝露をポリジメチルシロキサン（PDMS）によるpassive dosingでおこなってることと、毒性試験はユスリカ（in vivo）でおこなってることがこの論文の新しい点だそうです。

目指している方向性は共感出来ますが、毒性原因物質であることの確認が底質含有量と毒性（致死率・バイオマーカー応答）との相関というのは少し弱い気が･･･。考察で述べられているようにどのdosing手法を用いるかによって物質ごとのavailabilityが変化するので、passive dosingを用いずに原因物質かどうかを確認するステップは大事だと思います。

というか、passive dosingがbioavailabilityを模擬?できているのかという根本的な疑問があります。このあたりは勉強したいけど、「こういう操作をしたらこういう結果が得られた」と割り切った方が良いかも。

「殺虫剤耐性に関する変異を用いて環境試料の毒性原因を同定する」

Weston D.P., Poynton H.C., Major K.M., Wellborn G.A., Lydy M.J., Moschet C., Connon R.E., 2018, Using mutations for pesticide resistance to identify the cause of toxicity in environmental samples, Environ Sci Technol 52: 859-867.

面白い！

ピレスロイド系殺虫剤に耐性のある野生集団と耐性のないlab集団の両者を用いて汚染水への曝露試験をした時に、集団間で応答の違いがあれば、その原因はピレスロイド系殺虫剤と分かる、という話。このような、変異mutationによる感受性の違いを利用して、毒性原因物質を推測する手法をbiological TIEと名付けてます。

殺虫剤への耐性がある集団をlabで非汚染下で継代飼育しても、耐性はある程度維持できるそう。ピレスロイド系殺虫剤への耐性がある集団は、DDTに対しては高い感受性を有していました。この原因は不明みたいですが、交差耐性cross-resistanceの逆が生じているんですね。この現象を考えると、biological TIEというコンセプトは実現可能性が低い気もします（というかとりあえず名付けてるだけでは･･･）。ただ方向性は面白いし、AOPとEDAの発展した先は実はこんな感じではないかと妄想させてくれます。

1塩基の変異×2がピレスロイド系殺虫剤の耐性獲得につながってるそうです（Weston et al., 2013, PNAS）。

海外ポスドク目指しての就活 ～公募アプライ編～

この時の続き。

この秋Marie-Curie FellowshipとCanon European Fellowshipに応募しましたが、12月中旬にCanonからはお祈り通知が来ました。Marie-Curieもおそらくダメだろう、ということで12月末から公募を探し始めました（遅い）。で、実際Mari-Curieはダメでした

（→Marie-Curie fellowshipの結果が返ってきたよ）。

出した公募10個ほども全部ダメで、4月から国内の某ポスドク職として採用していただけることになりました。採用されるまでには本当に多くの方のお世話になり、ご迷惑をおかけしました。感謝しても感謝しきれません。

結局国内に就職するので、失敗の記録ではありますが、どこかの誰かの参考になるかもしれないので、海外ポスドク目指しての就活をここに記録しておきます。ただ就活の過程で、海外のポスドクになるには（海外学振などお金をこちらが持って行く場合を除けば）受け入れ先とコネクションがないと厳しいとか、日本で確固たる地位を確立してから海外へ出る方が帰国時の就職先の心配をしなくても良いのではないか、というお話を複数の方々から受けたので、そのことも念頭に置いて読んでもらえればと思います。

公募探しに利用した主なポータルサイトは以下。

・Research gate

・Nature jobs

・EuroScienceJobs

・jobs.ac.uk

個人的に良さげな公募が多く見つかったのは、Research gateとNature jobsでした。Research gateは自分の"スキル"に適した公募を勧めてくれるのが良いけど検索しにくい。Nature jobsは絶対数が多くて更新も頻繁っぽい。

生態毒性というニッチな分野でのポスドク情報は、ドイツのThe University of Koblenz Landauの研究室が運営しているEcotox blogというページがとても参考になりました。自分は利用しなかったけど、水関係の公募はJosh's Water Jobsというサイトに出てました（後輩のYさんに教えてもらった）。

大してコネクションの無いまま出した公募ですが、数件はSkype面接まで進みました。例えばこんなスケジュールでした。

・2018.01.04 公募締切

・2018.01.31 面接への招待メール

・2018.02.07 面接の際に使用するプレゼンスライドを送ってくれとのメール

・2018.02.09 面接本番

2/7のメールがなぜか迷惑メールに振り分けられていて、面接にスライドが必要なことを面接の10分前に知らされるという最悪のアクシデントが起きたこともあり、この面接は散々なことに…。それまでのメールは普通に受け取れていたのに、なぜ！？ ありあわせの資料を急遽つなぎ合わせて送信し、「だいぶシンプルな資料だね」みたいな感想を言われ、それでも45分くらいの面接を最後までしてもらえましたが…。

この面接で知ったのは、実験手法のかなり細かい部分まで聞かれるということです。もちろん「このポストでどういうことがやりたいか？」とか「何が求められていると思いますか？」みたいな良くある質問もありましたが、中には「YYYの正式名称は何？」とか「XXXをしたいときに用いる統計手法は？」とかクイズみたいな問いまで。

なぜそこまで詳しく聞くのか当時は分かりませんでしたが、おそらく海外ではテクニシャンがいたり分業が盛んだったりで、学生やポスドクが全ての作業を担うことは少ないのでしょう。だからカバーレターに「XXXができます！」と書いてても、実際は全くできないというケースが多いみたいですね。もしくは、自分の受け答えがあまりにお粗末だったために初歩的な質問が出されたのかも…。

Marie-Curie fellowshipの結果が返ってきたよ

結果は惨敗。サクラ散る。

Marie-Curie Individual fellowshipはEUがお金を出しているフェローシップです。獲得すると、日本の大学教授並みの給料がもらえ、プラス旅費と家族手当、研究費が手に入り、同業者から一目置かれる（であろう）高ステータスなフェローシップです。平たく言うと学振のEU版ですね。

去年の夏から申請書の準備をし、9月中旬に応募して、1月29日に結果が来ました。Marie-Curieの申請書はExcellence (weight 50%)・Impact (weight 30%)・Implementation (weight 20%) の3つのパートで構成されてますが、それぞれのパートの点数と長所・短所が記述形式で返ってきました。下のような感じです。

・Excellence Score: 3.60 (Threshold: 0/5.00 , Weight: 50.00%)

Strengths:
- The state-of–the–art in the project field is comprehensively presented and major gaps in the existing knowledge are specifically identified.
- The research objectives are clearly and logically presented.

(略)

Weaknesses:
- The appropriateness of the research approach is not convincingly demonstrated. The proposed methodology, especially in its experimental
part, is not presented in sufficient detail.

- The proposal does not comprehensively describe what specific measures will be taken to integrate the researcher within the host’s research
team.

科研費でも不採択の場合に審査員からの評価を教えてもらえますが、Marie-Curieは記述式の評価がもらえる分、より具体的に改善ポイントのわかる点が良いですね。申請書作りにも審査にも多大な時間が使われているのだから、科研費でもこれくらいのフィードバックがぜひ欲しいところです。

総合点は、100%中68%でした。70%が足切りなので惜しいところだったのかと思いきや、実際にpassするには85~90%は必要みたいです。なので全然惜しくない…。

評価を見てみると、アウトリーチやスケジュールの緻密さに重きを置いているのが分かります。一応簡単なアウトリーチの予定は書いてたのですが、"The quality of communication of the project results to different target audiences is not sufficiently demonstrated."とのこと。小中高への課外授業とかマジで書かないとプラスにならなさそう…。

気持ちを切り替えて、今は雇われポスドクへの道を模索中です。

相互情報量ベースのネットワーク分析 using R minet

先日WGCNA（Weighted Gene Coexpression Network Analysis）を試しに使ってみました（この記事 → 論文のメモ: 生態毒性研究へのWGCNAの適用）。しかし、遺伝子発現のデータは別に直線関係ばかりではないので（最近読んだ論文のメモ： omicsデータの用量応答モデル）、ピアソンの相関係数に基づいて計算することに違和感を覚えました。

そこで、非線形の関係も扱える相互情報量（mutual information）に基づくネットワーク分析を勉強してみました。

参考にした情報は以下の論文＋Wikipediaの"相互情報量"と"情報量"。

ARACNEをRで実装した論文。

Meyer PE, Lafitte F, Bontempi G, 2008, minet: AR/Bioconductor package for inferring large transcriptional networks using mutual information, BMC Bioinfor 9(1): 461.

相互情報量とピアソン相関係数のどちらが適切かは時と場合によるという話。

法隆大輔, 林武司, 2013, 遺伝子発現プロファイル類似度としてのピアソン相関係数と相互情報量の比較, 計量生物学 33(2):125-143.

まず相互情報量を求める

相互情報量の定義はwiki参照。数式で書くとこんな感じらしい。

(wikiより)

連続値の場合の周辺確率と同時確率の求め方は、色々提案されているようですが（たぶん）、一番分かりやすいのはデータを分割して離散化する方法でしょう。

下のデータで遊んでみる。G1とG2が遺伝子発現量を、Sが試料番号を示してます。

> data

　    S1 S2  S3 S4 S5   S6 S7  S8 S9
G1   62 41     5 30 15    52 45    1   39
G2 167 42 136 50 98  114 90 167  38

ピアソンの相関係数は-0.20です。

一方、相互情報量の場合。

まずデータを分割して離散化します。区切り方は、区切りの幅を同一にする"equalwidth"か区切り内のデータ数を同一にする"equalfreq"か、など選択します。

> d <- discretize ( t(data) , nbins=3 ) #デフォルトだとnbinsはサンプルサイズの平方根

> d

  G1 G2
1  3   3
2  2   1
3  1   3
4  2   1
5  1   2
6  3   2
7  3   2
8  1   3
9  2   1

## こういう感じに分割される

# このデータの場合equalfreqもequalwidthも同じ

次に相互情報量を下記の式から求めます。Hはエントロピー。

(wikiより)

上の図をもとに集計表を作成し、X（=G1）とY（=G2）のエントロピーを求めます。

• H(X) = -  { 1/3*log(1/3) + 1/3*log(1/3) + 1/3*log(1/3) }  = 1.584963
• H(Y) = -  { 1/3*log(1/3) + 1/3*log(1/3) + 1/3*log(1/3) }  = 1.584963
• H(X,Y) = - { 1/3*log(1/3) + 2*1/9*log(1/9) + 2*2/9*log(2/9) } = 2.19716
• I (X; Y) = 1.584963 + 1.584963  - 2.19716  =  0.972766 

一応、Rのinfotheoライブラリーのmutinformation関数で検算しておきます。求め方は"empirical (=maximum likelihood)" か"shrink entropy"か、"schurmann-grassberger"など選べます。

empirical methodで求めた相互情報量は0.973で、上記の手計算と一致していますね。

> mutinformation(d, method="emp") ## 底はe

          G1      G2
G1 1.099 0.674
G2 0.674 1.099

> natstobits ( mutinformation(d, method="emp") )  ## 底2のビットに変換

                  G1             G2
G1 1.5849625 0.9727653
G2 0.9727653 1.5849625

ちなみにRのminetパッケージ（infotheoを利用している）だとbuild.mim関数で一気に求められます。

> build.mim( t(data), estimator="mi.empirical", disc="equalfreq")
                   G1             G2
G1 0.0000000 0.6742695
G2 0.6742695 0.0000000

相互情報量をもとにネットワーク分析をしてみる

ぶっちゃけRのminetライブラリにぶち込むだけですが、とりあえずやってみます。

次はデータを拡張して、下のものを使用します。G1~G10が遺伝子で、各サンプルS1~S6における発現量のデータです。サンプル6個で相互情報量を求めるのは厳しいかも…。

ピアソン相関係数の行列はこんな感じ。

相互情報量を、build.mim関数で求めたのが下。

> mim.eqwi <- build.mim( t(data) , estimator ="mi.empirical", disc="equalwidth")

このように類似度行列が得られたら、次は図示するための準備。WGCNAの場合と同様に、閾値で区切って0/1データに変換するなどの処理を行います。とりあえずminetライブラリに実装されているARACNEを用います。

> net.eqwi <- aracne(mim.eqwi)

#赤字の場所がARACNEによってカットされた部分。indirectな関係を除外できるということですが…？

最後に、相互情報量と相関係数の結果を比較。snaライブラリーで図示。

>  mim.cor <- build.mim( t(data) , estimator ="pearson")  # 相関係数使用, 離散化はされない

> net.cor <- aracne(mim.cor)

> library( sna ) 

> par (mfrow=c(2,1) )

> gplot(net.eqwi, usearrows=FALSE, displaylabels=TRUE)

> gplot(net.cor, usearrows=FALSE, displaylabels=TRUE)

#左が相互情報量ベースで右がPearson相関係数ベース

論文のメモ: omics technologyとAOP

AOP（Adverse Outcome Pathways）の構築に、網羅的な生物応答の解析技術（omics）を用いて何ができるかというお話し。omicsの中でも特にtranscriptomicsについて。

「化学物質リスク評価へAOPを適用するうえでomicsの果たす役割」

Brockmeier EK, Hodges G, Hutchinson TH, Butler E, Hecker M, Tollefsen KE, Garcia-Reyero N, Kille P, Becker D, Chipman K, Colbourne J, Collette TW, Cossins A, Cronin M, Graystock P, Gutsell S, Knapen D, Katsiadaki I, Lange A, Marshall S, Owen SF, Perkins EJ, Plaistow S, Schroeder A, Taylor D, Viant M, Ankley G, Falciani F, 2017, The role of omics in the application of adverse outcome pathways for chemical risk assessment, Toxicol Sci 158: 252-262.

2014年リバプール大学でのワークショップをもとに書かれた総説。

現状ではomicsは生態リスクの評価をおこなうに十分な証拠を提供できていない（"omics datasets cannot provide sufficient evidence to characterize risk within ERA"）とか、omicsが既存の手法を完全い代替することはない（we ... do not fee that they will completely replace all approaches used in traditional risk assessments）とか、omicsの限界点に触れていて好感触。他に挙げられているomicsの課題は、標準法がないこととインフォマティクスの手法が確立してないこと。この辺りはいつもの議論。

omicsを使うことによって 、既知のパスウェイを確認するだけではなく、新規パスウェイを発見することもできる例としてAntczakら（2015）が引用されてます。

「トランスクリプトーム解析によるエマメクチン安息香酸塩の急性毒性メカニズム洞察」

Song Y, Rundberget JT, Evenseth LM, Xie L, Gomes T, Høgåsen T, Iguchi T, Tollefsen KE, 2016, Whole-organism transcriptomic analysis provides mechanistic insight into the acute toxicity of emamectin benzoate in Daphnia magna, Environ Sci Technol 50(21):11994-12003.

AOP構築に至るまでにomics解析をどう活用できるかの事例として読みました。この論文で提唱されたメカニズムはAOP wikiにも仮登録されてます。

既往文献から導き出せるAOP仮説をマイクロアレイ解析で検証してみたという印象。あくまでomicsは確認作業っぽい。

エクジステロイド受容体（EcR）に関するin vitro試験の結果が、単純なomicsだけではAOP構築できないことを示している気がします。エマメクチン安息香酸（EMB）に曝露させるとDaphniaのEcR遺伝子は活性化されるが、in vitro試験ではそのような反応は見られない。これは、EcRシグナル経路がEMB曝露に対する応答の下流にあるから。

因果関係を踏まえてAOPを構築するためには、key event間の関係をこういう風に一つ一つ検証していかないとダメなんですね。当たり前のことでしょうが…。自分にはin vitro試験の必要性を感じさせる事例論文として面白かったです。レポーターアッセイなどin vitroでなくともノックアウト生物を簡単に作れれば良いのかもしれない。

「生理学を理解するツールとしての機能ゲノミクスの多層式データ統合」

Davidsen PK, Turan N, Egginton S, Falciani F, 2015, Multilevel functional genomics data integration as a tool for understanding physiology: a network biology perspective, J Applied Physiol 120(3):297-309.

ネットワーク分析に関して上のBrockmeierら (2017) で引用されていた総説。FF氏が文章を書いたところはなんとなく分かってしまう。

前半の総説部分は、データからネットワークを推定しようという話。非線形のデータにも使える相互情報量に基づいたARACNEでの推定が良いよ、という話など。ネットワーク推定に生物replicatesは>50~100必要とのこと。

「ゼブラフィッシュ胚を用いた環境毒性評価のための縮小トランスクリプトーム」

Wang P, Xia P, Yang J, Wang Z, Peng Y, Shi W, Zhang X, 2017, A reduced transcriptome approach to assess environmental toxicants using zebrafish embryo test, Environ Sci Technol 52: 821-830.

ちゃんと読んでないしAOPともそれほど関係ありませんが、ampliseqを生態毒性の研究に用いた例としてメモ。omicsは標準試験法がないとのことでしたが、マイクロアレイに代わってampliseqが標準法になるかも。遺伝子ごとに用量応答反応のモデル（直線・シグモイド・U字型）を当てはめ、どの濃度範囲で応答するかを調べて、そのEffect Concentrationと生物学的機能とを絡めて議論しているのが面白いです。このETCの論文が引用されてます。k-meansクラスター法で同様の議論をしている論文は多いですが、この論文のアプローチはなんとなく新鮮でした。

(追記 2018.03.23)

Leung KM, 2017, Joining the dots between omics and environmental management, Integr Environ Assess Manag 14: 169-173.

IEAMのBrief Communication。あんまり真新しい話はない。

論文のメモ: トンネル洗浄水による遺伝子発現変動

「交通関係の汚染に曝露されたブラウントラウトの肝臓の遺伝子発現プロファイル」

Meland S, Farmen E, Heier LS, Rosseland BO, Salbu B, Song Y, Tollefsen KE, 2011, Hepatic gene expression profile in brown trout (Salmo trutta) exposed to traffic related contaminants, Sci Total Environ, 409(8), 1430-1443.

久しぶりに路面排水関係の論文。

トンネルの洗浄水に曝露させたブラウントラウトの肝臓の遺伝子発現をマイクロアレイで調べた論文。何が分かったかとなると微妙…。化学分析の結果、というか実験を始めるまでの想定で大体のことは分かってる気がする。一応adverse outcomeにつながるまでの想定pathwayを描いてます。

論文のメモ: 生態毒性研究へのWGCNAの適用

重み付き遺伝子共発現ネットワーク解析（WGCNA; weighted gene co-expression network analysis）のはなし。

「遺伝子共発現ネットワークによるDaphniaの繁殖影響予測」

Asselman J, Pfrender ME, Lopez JA, Shaw JR, De Schamphelaere KAC, 2017, Gene co-expression networks drive and predict reproductive effects in Daphnia in response to environmental disturbances, Environmental Sci Technol, in press.

殺虫剤とシアノバクテリアのbinary mixtures (全144サンプル) に曝露させたミジンコのマイクロアレイデータと産仔数データを、WGCNAで統合して分析してます。目的は、遺伝子発現データからmixture exposureの産仔数を予測したい、というもの。産仔数と相関のない遺伝子モジュールを省いていき、予測精度を上げていく流れ。最終的には一般化加法モデル (GAM) で予測。

Mixureだと産仔数の予測精度が低下するなど結果自体の新鮮さはあまりないかもですが、モデルの交互作用項で遺伝子発現レベルの複合影響を考慮するという手法など面白い。

「オオミジンコの初期転写反応は甲殻類特異的な遺伝子ネットワークに位置する」

Orsini L, Brown JB, Shams Solari O, Li D, He S, Podicheti R, Stoiber MH, Spanier KI, Gilbert D, Jansen M, Rusch D, Pfrender ME, Colbourne JK, Frilander MJ, Kvist J, Decaestecker E, De Schamphelaere KAC, Meester LD, 2017, Early transcriptional response pathways in Daphnia magna are coordinated in networks of crustacean specific genes, Mol Ecol, in press.

オオミジンコDaphnia magnaを12種のストレス（シアノバクテリア, 金属, 殺虫剤, 低pHなど）にそれぞれ曝露させ、RNA-seqした論文。オオミジンコは3つのgentotypeを使用してます。

この論文の結果の一つは、2011年ScienceのD. pulexゲノム論文の確認。 すなわち、環境ストレスによって発現変動した遺伝子は甲殻類特異的だということ。さらにgenotypeにも特異的だそうです。この論文では甲殻類crustaceansと言っているけど実際はD. pulexとD. magnaしか相同性検索に使用してないので、環境ストレス応答遺伝子が系統によって異なるのはDaphniaに限った話じゃないかと思ったのですが、C. elegansでも報告されている話らしい（Zhou et al., 2015, Genomics）。

ストレス曝露によって発現変動した遺伝子の解析（DEG解析）では、キチン・表皮の考察話が多い。これは正直知見が多くて、議論しやすいから議論されてるだけな気もしますが。Zebrafishのメタ解析の論文と同じく、DEG解析ではあまり多くのストレス応答遺伝子を見つけられなかったので（全処理で変動したのは1,396genes; およそ全3万遺伝子中）、追加でWGCNA解析をしてます。WGCNAでモジュールを調べることで、アノテーションのない遺伝子の機能を推測できるかも、という話ですが、いまいち結論が見えにくい…。データが膨大過ぎて解析が追いついていない感。

以下、WGCNAの勉強。

超絶初心者の自分用メモです。ご容赦ください。

共起発現ネットワークについては、英語版wikiが分かりやすい。発現量の相関係数sijを計算し、その相関係数がある閾値を超えたら1、閾値以下なら0というルールに従って隣接行列に変換する。

ただWGCNAでは、隣接行列は2値（0 or 1）ではなく、連続値（相関係数sijの累乗sij^a）を用いるそうです。ここが"weighted"たる所以だそうな。どうやって指数aを求めるかは元論文のZhang and Horvath (2005) を見ましょう（適当）。RのWGCNAパッケージだとpickSoftThreshold関数で求められるっぽい。

モジュール（=関連のある遺伝子の集まり）の決め方は普通の階層的クラスタリングでもOK。課題は、どこでクラスターを切るか。

以下、RのWGCNAパッケージを使った解析の練習。下のパッケージ開発の論文とその詳しいチュートリアル（ここ）とを見ながらやってみました。

Langfelder P, Horvath S, 2008, WGCNA: an R package for weighted correlation network analysis, BMC Bioinformatics 9(1): 559.

チュートリアルのデータは複雑なので、まずもっと簡単なデータとして先ほどの遺伝子共起発現ネットワークのwikiページのものを使います。

data <- matrix( c(43.26, 166.6, 12.53, 28.77, 114.7, 119.1, 118.9, 3.76, 32.73, 17.46, 40.89, 41.87, 39.55, 191.92, 79.7, 80.57, 156.69, 2.48, 11.99, 56.11, 5.05, 136.65, 42.09, 236.56, 99.76, 114.59, 186.59, 136.78, 118.8, 21.41), nrow=10, ncol=3)
rownames(data)<- c("G1","G2","G3","G4","G5","G6","G7","G8","G9","G10")

adj <- abs( cor(t(data) , method="pearson") )

隣接行列adjは次のようになる。

G1 G2 G3 G4 G5 G6 G7 G8 G9 G10
G1 1.000 0.234 0.610 0.706 0.029 0.353 0.860 0.998 0.971 0.366
G2 0.234 1.000 0.627 0.523 0.979 0.992 0.295 0.295 0.458 0.990
G3 0.610 0.627 1.000 0.992 0.774 0.525 0.929 0.559 0.405 0.513
G4 0.706 0.523 0.992 1.000 0.687 0.413 0.969 0.660 0.518 0.400
G5 0.029 0.979 0.774 0.687 1.000 0.945 0.485 0.092 0.265 0.941
G6 0.353 0.992 0.525 0.413 0.945 1.000 0.174 0.412 0.565 1.000
G7 0.860 0.295 0.929 0.969 0.485 0.174 1.000 0.826 0.714 0.160
G8 0.998 0.295 0.559 0.660 0.092 0.412 0.826 1.000 0.985 0.425
G9 0.971 0.458 0.405 0.518 0.265 0.565 0.714 0.985 1.000 0.577
G10 0.366 0.990 0.513 0.400 0.941 1.000 0.160 0.425 0.577 1.000

pairs(t(data),cex=3,cex.axis=2,pch=21,bg="grey")

始めは、重みづけなしの共起発現ネットワークで処理する。0.8を閾値とします。

# wiki通り0.8を閾値とする

adj2<- ifelse ( adj > 0.8, 1,0)

# sna packageで描画

library(sna)

gplot(adj2,usearrows=FALSE, displaylabels=TRUE)

試しに閾値の取り方を変えてみると、ネットワークのつながり方も大きく変化します。 

次にWGCNAパッケージを使用。

sft = pickSoftThreshold(data, dataIsExpr = TRUE)

# Sclale-free topology fit indexを描画

plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,col="red")

20乗（power=20）してもScale-free topology fit indexが0.4程度と低いまま。FAQにその対処法が書かれてますが、今回は相関係数の絶対値を取っているunsigned networkなので、FAQに従いpowerを9とします。

このあたりから、かなり手探り状態…。

adj_w <- adjacency ( t(data),power=9 )

# snaで描画してみる, 辺の太さは隣接行列の大きさに従う

gplot(adj_w, usearrows=FALSE, displaylabels=TRUE, edge.lwd=adj_w*2)

太い辺edgeに注目したら、上のunweightedの場合のネットワークとほぼ同じですね。

次にクラスタリング。

TOM = TOMsimilarity(adj_w)
dissTOM = 1-TOM

geneTree = hclust(as.dist(dissTOM), method = "average")
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",hang =-1)

# 枝の剪定

dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize=2)
dynamicColors = labels2colors(dynamicMods)

dynamicColors
[1] "blue" "turquoise" "brown" "brown" "turquoise" "turquoise" "brown" "blue"
[9] "blue" "turquoise"

最後にeigengene。Eigengeneとは、各モジュールにおける第一主成分のことらしいです。モジュールを代表するような発現だと考えて良いみたい。

MEList = moduleEigengenes( t(data), colors = dynamicColors)
MEs = MEList$eigengenes
MEDiss = 1-cor(MEs)
METree = hclust(as.dist(MEDiss), method = "average")

plot(METree, main = "Clustering of module eigengenes", hang=-1)

#blueはGene1,8,9で、brownはGene3,4,7。

疲れたので、とりあえずここまで。

あとは、遺伝子発現以外の情報（生態毒性だと致死率・繁殖阻害・成長阻害とか）と各モジュールとの関係を見るなどの解析をして、面白いモジュールを見つけるのが王道でしょうか。

ちなみに同じ発現データの相関を取ってMDSプロットしてみた場合。

library(MASS)
p<- isoMDS( 1-abs(cor(t(data))), k=2, maxit=100)
plot(p$points, type="n",cex.axis=1.5,cex.lab=1.5)
text(p$points, rownames(p$points), col="red",cex=1.5)