自分のパソコンに統計解析ソフトRをインストールしておいてください。
方法は、「生物統計解析基礎」の項目１および２を参照してください。
0. はじめに
生物の多様性を研究し、それを示すには、「なにを比較するか？」を考えなければなりません。
その比較対象が、生物進化にともなった場合、生物系統分類に大きな影響を与えます。
生物進化とともに継承されてきた遺伝情報の比較は、生物の多様性解析では、最も重要なものであることは言うまでもありません。
例えば、カール・ウーズによるリボソームRNAやそのコード遺伝子の比較による生物系統の根本的な見直しなどの成果を挙げてきました。
その背景には、比較対象の生物に、その遺伝情報が存在している必要があります。リボソームは全生物に存在しています。
DNAのシークエンスを決めるために時間や労力がかかった時代は過ぎ去り、現在では全ゲノム情報を取得することが容易にできます。
すなわち、遺伝情報の比較は、ゲノムレベルでの比較となりました。
ここでは、遺伝情報の総体であるゲノムを比較するため、特定の遺伝子を持っていないということも情報として扱いますし、遺伝子という枠組みにとらわれずに塩基配列そのものを比較することも扱います。
ゲノムの比較には、大きく分類して、遺伝子アノテーションを必要とするもの（下記、項目５など）と必要としないもの（項目４など）に分けられます。
1. 解析の場の準備
C: においてgenomeというフォルダを作成し、そこにダウンロードしたすべてのファイルを入れてください。
Rを立ち上げ、getwd( )と打ち込みworking directoryを確認してください。
working directoryをC:/genomeに変更してください（ディレクトリの変更はsetwd("C:/genome")とすることもできます）。
それが面倒な方は、作業場所をデスクトップにして、ファイルをそのままデスクトップに置いて実行してください。
その場合、Rを立ち上げて、「ファイル」をクリックして、「ディレクトリの変更」をクリックし、デスクトップを選択し、OKとしてください。
2. DNAの塩基組成（GC含量）
DNAは二重らせん構造をとっていますが、二つの鎖はアデニン（A）とチミン（T）およびグアニン（G）とシトシン（C）がそれぞれ水素結合してニ重鎖となっています。
よって、AとTの数は等しく、GとCの数も等しくなっています。
このことはワトソンとクリックによるDNAのニ重らせん構造の発表より以前にシャルガフ（Erwin Chargaff）によって示されています。
ゲノムDNAを比較する際（特にバクテリアのゲノムDNAの場合）、全塩基に占めるGとCの割合が使われることがあり、GC含量と呼ばれています。
もちろん、GC含量とAT含量の和は1（100%）となっています。
それでは実際にゲノムのGC含量を調べてみましょう。
NCBIのトップページにアクセスしてください。
ここではバクテリアにおけるクロモソームとプラスミドのDNAのGC含量を比較します。
Pseudomonas resinovoransのクロモソーム情報は
http://www.ncbi.nlm.nih.gov/nuccore/NC_021499.1
ここにおいて、FASTAの塩基配列情報をダウンロードしてください。
次に、そのファイルをメモ帳などで開き、第一行の説明部分を削除し、塩基配列だけにして、AおよびTを0 [数字の0とスペース]、CおよびGを1 [数字の1とスペース]に置換してください。
わからないひとは、生物統計解析基礎の9. 行列とデータフレームを参照してください。
この0と1からなるファイルをchr.txtとしてフォルダgenomeに保存してください。
Pseudomonas resinovoransのプラスミド情報は
http://www.ncbi.nlm.nih.gov/nuccore/NC_021506.1
クロモソームのときと同様にして、pla.txtとしてフォルダgenomeに保存してください。
chr.txtおよびpla.txtのファイルには、0と1がスペースを空けて並んでいるものとなっていますので、それらを確認してください。
それではRを使ってそれぞれのファイルを読み込みGC含量をみてみましょう。
chr <- scan("chr.txt")
pla <- scan("pla.txt")
と打ち込んでください。それぞれのファイルがchrおよびplaとして読み込まれます。
GおよびCの数を知るためには、総和を求めればよいので、
sum(chr)
sum(pla)
で知ることができ、また、全塩基数は
length(chr)
length(pla)
により知ることができます。
よって、GC含量は
sum(chr)/length(chr)
sum(pla)/length(pla)
とすれば、求めることができます。
クロモソームのGC含量よりも若干プラスミドの方が低くなっていることがわかります。
この現象は上記のバクテリアに限ったものではなく、ほとんどのプラスミドとその宿主クロモソームにおいて成り立っています。
詳細を知りたい方は下記の論文などを参照してください。
http://www.hindawi.com/journals/ijeb/2012/342482/
上記と関連いたしますが、外来性のDNA領域はクロモソームにおいても周りよりも低いGC含量となる傾向があります。
細胞に水平伝播したものが低いGC含量のDNAであり、それがクロモソームへ挿入されたことを考えると理屈に合っています。
では、プラスミドのGC含量は領域における偏りはないのでしょうか？それらを確かめるためには、GC含量をある幅をもってスキャンし、その偏りを見る必要があります。
そこで、先のpla.txtにおけるGC含量の分布をスキャンしてみましょう（chr.txtも同様にできますが、少し時間がかかるため、ここではプラスミドの方で行います）。
例えば、塩基ポジション1から500までのGC含量、2から501までのGC含量、・・・と500塩基の幅でGC含量をスキャンします。
pla500 <- c()
for(i in 1:(length(pla) - 499)){
pla500[i] <- sum(pla[i:(i+499)])
}
この結果、ほしいデータがpla500に収納されます。
GC含量の偏りの様子を図で見る場合には
plot(pla500/500, type="l", xlab="Position", ylab="GC content")
としてください。
このプラスミドのサイズは198965塩基対です。
連続する500塩基が残りの領域に比べ有意に違ったGC含量となっているかどうかをχ2（カイ2乗）検定で調べてみましょう。
500塩基中i塩基がGまたはCの場合、
a <- c()
for(i in 1:499){
a[i]<-chisq.test(matrix(c(i, 500-i, 112171-i, 198965-112171-500+i), 2))$p.value
}
b <- data.frame(1:499, a)
range(b[b[,2] >= 0.05, 1])
これによって、iが260塩基から304塩基の範囲でp値が5%以上（信頼区間）になることがわかります。
segments(0, 260/500, 200000, 260/500, col="red")
segments(0, 304/500, 200000, 304/500, col="red")
すなわち、0から259塩基および305塩基から500塩基の場合には帰無仮説を棄却することになります。
課題１３ 上記は500塩基中の範囲でしたが、有意水準５％のとき、帰無仮説が棄却されない信頼区間のGCの数を、1000塩基中の範囲の場合には、どのような結果になるか示してください。
このことがわかれば、統計的にほかの領域とは有意にGC含量が異なっている領域を抽出できます。
プラスミドとクロモソームでは、そのような領域の出現頻度は違っているでしょうか？
重要なところですので、補足いたします。
それぞれの生物によって、そのゲノムDNAのGC含量は違っています。
また、ある幅の領域におけるGC含量をスキャンして、そのばらつきを見ても、生物種によって違っています。
その状況の中で、「外来のDNA領域は、内在のDNA領域よりもGC含量が低くなっている」ということで、どの領域が候補として抽出できますか？ということが問われています。
この問いへの答えはもちろん一つではありません。
それをいろいろと考えることがバイオインフォマティクスであると考えています。
3. GC skew
1996年、Lobryによってリーディング鎖とラギング鎖におけるGとC間およびAとT間の偏りが報告されました（Mol Biol Evol, 13, 660-665）。すなわち、リーディング鎖ではGおよびTに富み、ラギング鎖ではCとAに富むというバイアスです。この偏りを示す指標をGC skewと呼び、
（Cの数－Gの数）/CとGの数の和
で表します。よって、この数値は複製開始点を挟んで正負が逆になります。
上記のプラスミドについて、GC skewの様子を調べましょう。
GC含量のプロファイルを調べる際にはGもCも1としましたが、GC skewの場合にはそれでは都合がわるいですので、Cのみ1（あとは0）およびGのみ1のファイルを作成します。
それらをplaC.txtおよびplaG.txtとします。
スキャンする幅を1000塩基とすると、
plaC <- scan("plaC.txt")
plaG <- scan("plaG.txt")
plaCG <- c()
for(i in 1:(length(plaC)-999)){
plaCG[i] <- (sum(plaC[i:(i+999)])-sum(plaG[i:(i+999)]))/(sum(plaC[i:(i+999)])+sum(plaG[i:(i+999)]))
} #時間がかかります。クロモソームの場合はさらに時間がかかります。
plot(plaCG)
として、GC skewの様子を見ることができます。
4. ゲノムシグネチャ
数塩基配列のゲノムにおける出現頻度をゲノムシグネチャ（genome signature）といいます。
例えば、2塩基配列の場合には4×4=16通りがありますが、AAの相補鎖はTTですので、AAの出現頻度とTTの出現頻度は同じです。
そこで、2重鎖において両鎖の5'末端側から16通りのすべてのパターンの出現回数をカウントします。
この作業をゲノム塩基配列に対してスキャンすることによって行うと時間がかかりますので、行列にして処理することが速く解析する方法です。
すなわち、4塩基を1, 2, 3, 4など適当な数で置換して、そのベクトルを1塩基ずらして並べ、2列の行列にします（3塩基配列の場合には、さらに1塩基ずらしたものをもう1列追加して3列からなる行列にします）。
環状DNAの場合には、最後の部分に最初の部分を付け足した形で行列を完成させます。
その後、この行列に対して、各塩基配列に対応する組み合わせの行をカウントすれば終了です。
行数のカウントについては生物統計解析基礎を参照してください。
例えば、4種のバクテリアBacillus subtilis、Escherichia coli、Mycoplasma genitalium、Streptomyces griseusの2塩基配列のゲノムシグネチャは次の表になります。
このデータをRに入力します。
ec <- c(832142,389101,473369,686598,561841,386030,413156,473369,551175,508090,386030,389101,436052,551175,561841,832142)
bs <- c(677352,512270,473998,619638,647389,541810,693340,473938,534535,767862,541810,512270,423922,534535,647389,677352)
mm <- c(165916,60071,66470,103798,73024,33037,11290,66470,55959,34754,33037,60071,101356,55959,73024,165916)
sg <- c(275030,947123,815257,336016,827183,2010163,2519900,815257,1110967,2104250,2010163,947123,160246,1110967,827183,275030)
それぞれの生物種における塩基数が異なるため、単位ベクトルにして比較します（単位ベクトルは長さが１のサイズになっています）。大腸菌の場合には、
ec/sqrt(sum(ec^2))
になります。よって、先ほどの１６行４列のデータにおいてそれぞれを単位ベクトルに補正したものは下記となります。
gs <- data.frame(ec/sqrt(sum(ec^2)), bs/sqrt(sum(bs^2)), mm/sqrt(sum(mm^2)), sg/sqrt(sum(sg^2)))
課題１０ 上記においてそれぞれを単位ベクトル化した理由を述べてください（単位ベクトル化せずに比較するとはどういうことかを考えてください）。
gsは行が2塩基のパターン、列に生物種となっています。こちらをそのまま解析したものと、転置行列を行ったのちに解析したものを比べてください。
plot(hclust(dist(gs)))
plot(hclust(dist(t(gs))))
ヒートマップは
heatmap(t(gs))
階層的クラスタリングについては生物統計解析基礎の「16. クラスタリング」を参照してください。
ゲノムシグネチャに基づく比較は、遺伝子の枠組みを必要としないことより、メタゲノムのようなヘテロなデータ解析に使用することができます。
5. 遺伝子アノテーションに基づく系統解析
同じ祖先を持ち、機能が同じ遺伝子はオルソログと呼ばれています。
オルソログは遺伝子産物の構造に類似性があることより、類似配列検索において高く保存されているものをオルソログと便宜的にする場合があります。
これらオルソログ遺伝子を比較することはそれほど難しいものではありませんが、すべての遺伝子を搭載しているゲノムを比較することは一筋縄ではいきません。
ゲノムDNAを比較する方法は大きく分けると２通りあります。
１つはマルチプルアライメントに基づく比較であり、もう１つはマルチプルアライメントによらない比較です。
前者の方法にはオルソログ遺伝子のゲノムにおける有無のパターンのマルチプルアライメントに基づく方法と、共通するオルソログ遺伝子を抽出し、それらの配列のマルチプルアライメントを取り、すべてのオルソログのマルチプルアライメントを連結させたものを作成し、それに基づく方法に分けることができます。
後者の方法には数塩基の塩基配列の出現頻度（ゲノムシグニチャ）の類似度に基づく比較があります。
5.1. 類似配列の有無パターンに基づく比較
オルソログを特定の生物種から抽出し、それらをデータベース化して公開しているサイトがあります。
Microbial Genome Database for Comparative Analysis
このサイトでは、様々な生物分類群におけるオルソログを知ることができます。
では、自分でオルソログを抽出するには、どのようにすればよいでしょうか？
オルソログに対応することばは、パラログです。
オルソログは、進化的起源が共通で、機能が同じですが、パラログは、オルソログ集団と進化的起源が共通で、機能が違っています。
塩基配列からの推定遺伝子に基づく解析では、類似度を比較できますが、機能を知ることはできません。
よって、例えば、異なるバクテリアのゲノムに推定されたすべての遺伝子産物（アミノ酸配列）間での総当たりの類似度を算出します。
その方法は、類似配列検索の「2. 複数の配列に対する類似配列を一度に検索する」を参考にしてください。
データベースとqueryを変えても、変化ないトップヒットのペアはオルソログとしてもよいという感じになっています。
機能が推定である以上、オルソログの定義に従うと「感じ」としか言えません。
例えば、類似配列検索では、Deinococcus grandisのRodZと表現していますが、機能が明確になるまで、RodZ類似タンパク質となります。
私たちは、この遺伝子破壊株を取得し、細胞形態が桿状から球状に変化したことから、RodZホモログとしました。
ある機能を持つ遺伝子が存在しているとき、それと同じ機能を持つ遺伝子が異なる生物にあるかどうかということです。
例えば、４種の異なる生物種において、異なる機能の１０の遺伝子の分布の様子が下記のようであったとします。
このデータをもとにして、Molecular Evolutionary Genetics Analysis (MEGA)を用いて解析を行うこともできますが、ここでは4.でRで行った方法でやってみましょう。
「あり」を１、「なし」を０として、それぞれの情報をベクトルにします。
d1 <- c(1, 0, 1, 1, 0, 0, 1, 1, 1, 1)
d2 <- c(1, 1, 1, 1, 1, 1, 0, 0, 1, 1)
d3 <- c(1, 0, 0, 1, 1, 0, 1, 1, 1, 1)
d4 <- c(1, 1, 1, 0, 1, 1, 1, 1, 1, 1)
4.の解析と同様にして、ヒートマップを作成できます。
課題１１ 上記の方法では、遺伝子の有無というパターンになっていますが、実際のゲノムにおいては、特定の機能を持つ遺伝子が複数存在していることがあります。そのような時、上記の「遺伝子」は「遺伝子群（あるいは遺伝情報）」となります。その上で「０」「１」に数値化した場合、問題点がないとは言えません。どのような問題があるか答えてください。
5.2. オルソログ遺伝子産物の連結マルチプルアライメントに基づく比較
分子系統解析にあるように、類似配列を多重に並べて整列させ、マルチプルアライメントを得ることができます。
オルソログは異なる生物種に存在している類似配列ですので、それが可能となります。
例えば、上記の遺伝子の有無の表において、遺伝子１はすべての生物種にありますので、その産物であるアミノ酸配列を多重化できます。
遺伝子１の系統関係を示す際には、生物１に１０のオルソログがあり、生物４には１つしかなくても問題ありません。
しかし、すべてのオルソログを連結させる場合、具体的には、マルチプルアライメントを連結させる際、生物種において数が違っているとできません。
マルチプルアライメントを連結させるためには、各生物種からは１つのオルソログで行う必要があります。
よって、遺伝子の有無表において、オルソログ１ついう情報の遺伝子に基づいて、全オルソログ連結データを作成することできます。
そのデータをMolecular Evolutionary Genetics Analysis (MEGA)を用いて解析すればできます。