ここでは、種絶滅に伴ってどれだけ遺伝的機能多様性が失われるかという生態学的問題（生物多様性ー生態系機能関係）を例に、簡単なモンテカルロシミュレーションについて学びます。

元ネタは、

Takeshi Miki, Taichi Yokokawa, Kazuaki Matsui (2014) Biodiversity and multifunctionality in a microbial community: a novel theoretical approach to quantify functional redundancy. Proceedings of The Royal Society B(OPEN ACCESS!!) vol.281 no.1776 20132498

その日本語解説

定量生態学のすすめ：生物多様性減少の影響を定量的に予測する

[1] 生物多様性と生態系機能の関係について学ぶ

まずは、上の日本語解説の１～６のセクションをよく読みましょう。

[2] 全ゲノム情報を用いた多機能性評価

上の日本語解説の７のセクションの図３のようなイメージのグラフをシミュレーションによって描くのがゴールです。

特に第二段落目「具体的にどのように～群集全体での機能は低下にしくいだろうと予測できるわけです」に出てくる、キーワード「オーソログ遺伝子」、「細菌群集を構成する各細菌種間で，オーソログ遺伝子の共有の度合い」というところの意味をよくかんがえてください。

[2-1] オーソログ遺伝子リストの読みこみ

まずは次のファイルをダウンロードして解析用のフォルダに保存してほしい。

matrix_part.csv

このファイルをエクセル等の表計算ソフトで開く（もしくは次のコードでRに取り込む）と０、１の数字が並んでいるがそれの意味は以下の図で説明しよう。

gene_data <- read.csv("./matrix_part.csv", header=F)

View(gene_data)

[2-2] シミュレーションで描きたい図（ゴール）

このcsvファイルのデータを使って、全52種類が存在する（仮想的な）細菌群集から次々に種が絶滅していったときに群集全体からはどの程度オーソログ遺伝子が失われるかをRを使ってシミュレーションしたい。そのゴールとなる図は日本語解説の図３のようなものである（使っているデータが同じではないので全く同じ図になるわけではないことに注意）。

種数20のところで説明したステップを、種数について、１～５２で繰り返せば、(SR, MF)の組がたくさん得られるので、それを散布図にすれば、図３のようになるはずだろう。

[2-3] Rを用いて乱数を用いた計算（＝シミュレーション）をするには？

１）乱数発生の出発点を決める（何千、何万の乱数を繰り返し発生させる時でも一度だけ実行する）ために、set.seed()関数を使う。引数は任意の自然数を用いるとよい。これを乱数の初期化という。

set.seed(1234)

Rに限らず多くのコンピュータ言語で使用される乱数は、正確には「擬似乱数」であり、隣り合う値がほとんど相関しておらず、どの値も同じくらいの頻度で登場する、非常に周期の長い周期関数である。周期関数なので、出発点が同じならば常に同じ乱数のセットが発生する。したがって、蘭州発生の出発点を固定しておけば、シミレーションの結果に再現性が保たれることになる。ちなみに現在主流の高品質の擬似乱数発生器はメルセンヌ・ツイスター法という「メルセンヌ素数」に関連する超絶に長い周期（2^19937-1→10の6,000乗）を実現可能である。擬似乱数についてちょっとだけ聞いたことがある人は、周期関数であることを知っているので、たくさんの乱数を発生させる時には周期が一周してしまうことを心配して、初期化を繰り返してしまうミスを犯す。初期化は一回だけと覚えておこう。

２）次に、(1,2,...,52)から任意の数を重複なく選ぶ方法であるが、これについて純粋に[0,1)の一様乱数を複数発生させることによって達成可能ではあるがそのアルゴリズムをゼロから考えるだけでも日が暮れてしまうので（それはそれで興味深いのだが）、今回はRに標準で用意されているsample()という関数を使ってみよう。

１～５２までが格納されているベクトルを作って、そこから重複なく２つの数字をランダムに選ぶ方法は次のコードで実現可能である。

species_ID <- c(1:52)

sample(species_ID, 2, replace=F)

３）したがって、sample()の結果を別のベクトルに格納しておけば、

sample_ID <- sample(species_ID, 2, replace=F)

このsample_IDの各要素iを使って、ランダムに選ばれた種の遺伝子リストにアクセスできる。

i <- 2

gene_data[, sample_ID[i]

４）最終的には、1～52種が含まれるランダムな群集をそれぞれの種数について20個ずつ作って、上の散布図と同じスタイルのグラフを作ってほしい（ただし、注意してほしいのは、全部で52種類しかいないので52種類が含まれる群集は１つしか作れない）。この作業を自動的に完遂するには、for loopを何度も使う必要があることをイメージしてほしい。ただその前に、コアとなる計算方法について考える必要がある。

５）一つの種からなる群集に含まれる遺伝子数はどのようにカウントすればよいだろうか？ 次のステップとして２つの種からなる群集に含まれる遺伝子数はどのようにカウントすればよいだろうか？

• 一つの種からなる群集であれば、各遺伝子の有無は０、１で表されているので、sum()関数をうまく使って、注目する種の遺伝子リストに対応するベクトルのすべての要素の総和をとれば、存在する遺伝子数になるだろう。
• では、二つの種からなる群集の時はどうか？ 簡単のため、種Aの遺伝子リストを(1,1,0,0,0)、種Bの遺伝子リストを(1,0,1,1,0)であるとすると、群集全体での遺伝子の有無は(1,1,1,1,0)となるので遺伝子数は４と計算されないといけない。実は、「比較演算子」とsum(), ベクトルの足し算(加法)を使えば、結構簡単にできます。それが思いつかなければ、地道にif文などを使って、各遺伝子の有り無しを1～10000までチェックしていく必要があるだろう。

６）最終的にシミュレーションの結果は、二列（一列目が種数、2列目が遺伝子数）のデータフレームに保存されるのが望ましい。そうすれば、散布図も容易に書けるし、日本語解説の図3の下数行後にあるべき乗関数での近似も、両辺自然対数（底がe）変換したデータについて線形回帰をすれば求まるはずである。

というわけで今週の課題は、以下のことを実行してそのRスクリプトを提出してください。（わざわざできたグラフをレポートにまとめる必要はありません）

• 1～52種が含まれるランダムな群集をそれぞれの種数について20個ずつ作って各群集の遺伝子数をカウントし、図３のようなグラフを作成する
• 両側対数の線形回帰によって、（MF：遺伝子数）= c x (SR: 種数)^a の形式のべき乗近似をして散布図に重ねて表示する