Rによる主成分分析

主成分分析（Principal Component Analysis）は、多数のデータを俯瞰的に見るためのツールです。今回はフリーで使える統計解析ソフトのRを使って、主成分分析を行う方法を説明します。

データの準備

図1

ここでは例として有名な アヤメという植物のデータで説明します。まずExcelなどを用いて図1のようにデータをまとめます。図のとおり1行に1個体、各列にそれぞれの形態データという形でまとめてください。このデータには、sepal.length（がく片の長さ）、sepal.width（がく片の幅）、petal.length（花弁の長さ）、petal.width（花弁の幅）の4つの量的変数と、Species（種、setosa, versicolor, virginicaの3種類）という質的変数が含まれています。まとめ終わったら、「ファイル」>「名前を付けて保存」を選択し、csv形式で保存してください。 ちなみに、ファイル名は半角英数字と「_」だけで構成することを強くおすすめします。

図2

　csvとは、Comma-Separated Valuesの略で、Excelのそれぞれの列を「,」で区切った形となります。保存したファイルをメモ帳などのテキストエディタで開くと、図2のように表示されます。この形式はRにデータを読み込ませる上で重要なので覚えておいてください。

データの読み込み

データが準備できたら、Rを起動します。起動できたら以下のコマンドでデータを読み込ませます。

x <- read.table("C:/Users/shinka/Documents/iris.csv", header=T, sep=",", na.strings="NA", 
dec=".", strip.white=T)

それぞれのコマンドの意味は、以下の通りです。

• x <- read.table():　()内の表を、x（任意の文字）に代入します。
• "C:/Users/shinka/Documents/iris.csv":　読み込むファイルの場所と名前を指定します。ダブルコーテーション(")で囲ってください。
• header=T:　ファイルの先頭行をラベルとして扱います。今回の例では先頭に変数名（sepal.lengthなど）を書いているので、T(=True)です。
• sep=",":　列の区切りに「,」を使っていることを示します。ダブルコーテーションで囲ってください。
• na.strings="NA":　空欄がある場合に「NA」という文字で埋めます。これがないと空欄に対して、システムが勝手な解釈で値を振ってしまうそうです。
  ダブルコーテーションで囲ってください。
• dec=".":　小数点記号を指定します。ダブルコーテーションで囲ってください。
• strip.white=T:　データに含まれる空白文字を自動で削除します。

上のコマンドを打ち終えたらEnterを押してください。もしエラーがでた場合はコマンドのどこかが間違っています。確認してください。 読み込みができたら、先に指定した変数名（ここでは「x」）を入力してEnterを押してみましょう。読み込ませた表が表示されるはずです。さらに、以下のコマンドで表の要約を見ることができます（図3）。

summary(x)

図3

主成分分析の実行

メインの主成分分析に移ります。まず、今回主成分分析に使うprcompというコマンドは質的変数（今回のデータではSpeciesが該当する）に対しては使用できません。そこで、「Species」の列を削ります。具体的には、「Species」以外の列を新しい変数に代入します。

y <- x[,-5] 

上のコマンドはxに入っている表の5列目以外をyという新しい変数に代入することを意味しています。[]内のコンマを忘れないでください。もし3列目を削りたければ、「-3」とします。
　実際にprcompを使います。ここでは、主成分分析の結果をpcaという変数に代入します。summary(pca)で結果の要約が見れます（図4）。

pca <- prcomp(y, scale=T)
summary(pca)

図4

図の「Proportion of Variance」は寄与率と呼ばれるものです。簡単にいえば、「一つの主成分軸で、データ全体の何割を説明できるか」を示しています。図の例では、PC1という主成分軸で全体の70%以上が説明可能です。また、下の「Cumulative Proportion」は累積寄与率で、寄与率を第一主成分から順に足していったものです。PC2の寄与率は22.85%なので、第二主成分までで、データの95.81%が説明できます。
　最後に主成分分析の結果を図示します。図示にはggbiplotを用います。これはラボのパソコンによってはインストールする必要があります。インストールは以下のコマンドで可能ですが、わからなければ先生か誰かに聞いてください。

install.packages("devtools")
devtools::install_github("vqv/ggbiplot")

ggbiplotを以下のコマンドで呼び出してから実行します。ここで変数xに入れているSpeciesを利用することで、種ごとに色分けすることができて便利です（groups=x$Speciesで指定する。x$Speciesは変数xのSpeciesというラベルを意味する）。

library(ggbiplot)
ggbiplot(pca, obs.scale=1, var.scale=1, groups=x$Species)

図5

結果は右のようになります（図5）。アヤメの種ごとにはっきりと分かれていることがみてとれます。矢印は各形態の方向性を示しています。これを見るとpetal.lengthとpetal.widthはともに第一主成分軸に沿っており、両者は「ほぼ同じ意味をもつ形態である」といえます。また、sepal.lengthも少し逸れてはいるものの、petalとよく似た方向性であり、主にこの3つが第一主成分軸を構成していることがわかります。一方で、sepal.widthはこの3つとは全く別の方向へ伸びており、これにより第二主成分軸が形成されています。以上より、ざっくりとですが、アヤメの分類では、sepal.widthとその他の形態というまとまりでデータを見ることができそうであることがわかりました。

以上が主成分分析の一連の方法となります。主成分分析について勉強したい方は以下のブログがわかりやすいと思います。

https://logics-of-blue.com/principal-components-analysis/

BEAST2による分岐年代推定②

Tracerによる収束の判定

図1

BEASTによる解析が終了すると、logファイル（*.log）とtree（*.trees）ファイルができています。解析結果が収束しているかを確認するためには、logファイルをTracerに読み込ませる必要があります。Tracerを起動したら「+」をクリックするか、「File」>「Import Trace File」でlogファイルを選択します。
　
　図1がTracer の画面です（表示されるパラメータは解析の条件によって異なります）。この解析では200,000,000回の試行を行っています（図1の「States」）。また、MCMCでは、試行の最初の方は数値が安定しないため結果から省きます。これが「Burn-In」で、この例では全試行の10%にあたる20,000,000回が指定されています。
　

図2

続いて各パラメータが収束しているか見ていきます。例として「birthRate」を見てみます。左側のパネルで「birthRate」をクリックし、右側のパネルで「Trace」タブを表示してください（図2）。表示されたグラフの横軸は試行回数で、縦軸は各試行での推定値です。グラフの左端の白くなっている部分はBurn-Inによって省かれたところです。このグラフでは、Burn-Inによって削られて以降、推定値が一定の範囲に収まっているのがわかるかと思います。次に左側のパネルで「ESS」（有効サンプルサイズ）を見ると、「birthRate」では1714と大きな値となっています。この値はパラメータの収束に関する基準で、200を超えるようにすべきとされています。以上より、「birthRate」は十分に収束していると判断できます。

図3

　悪い例を見てみます。これは「prior」 のパラメータですが、グラフを見ると試行により推定値が一定の範囲に収まっていないことがわかると思います（図3）。ESSの値も61と低く、収束しているとは判断できません。このような場合、試行回数をさらに増やしたほうが良いでしょう。試行回数を増やしてもESSが改善されない場合には、パラメータの初期値などを変更したほうが良いかもしれません。


Treeファイルの要約

全てのパラメータの収束が十分であると判断されたら、treeファイルの要約を行います。BEASTの解析で得られたtreeファイルには各試行の結果が含まれています。すなわち、先の例でいえば、200,000,000 / 1000（どのくらいの頻度でファイルに記録するかという値。xmlファイル作成時に変更可能。）= 200,000本分の系統樹と分岐年代推定値が記録されています。当然一つ一つ見ていくわけにはいかないので、treeファイルの内容を要約します。これにはTreeAnnotatorを用います（BEASTと同じフォルダに入っています）。

図4

TreeAnnotatorを起動します（図4）。「Input Tree File」でtreeファイルを選択します。「Output Tree File」では出力するファイル名を入力します。「Burnin percentage」では何%削るかを指定します。これはTracerにおいて収束が確認された際のBurninに合わせます。200,000,000回の試行、20,000,000回のBurninで収束しているならば、10と入力してください。「Target tree type」はどの樹形で出力するかを指定します。基本的にはそのままで大丈夫でしょう。「Node heights」は「Mean heights」 にしておけば良いと思います。入力ができたら「Run」をクリックします。要約が完了したら出力ファイルをFigTreeに読み込ませて結果を確認してください。以上がBEASTによる分岐年代推定の基本的な方法となります。

2021/12/14 追記

「prior」が極端に低い場合の対応について

ほとんどのパラメータは収束しているのに「prior」と「posterior」のESSが低いままのことがしばしばあります。「posterior」は「prior」と「likelihood」に関係してくるので、実質「prior」のみの問題と考えて良いでしょう。グラフを見ても収束の兆しが見えず、試行回数を増やしても改善されない場合、推定するパラメータ数が多すぎることに起因している可能性があります。対処方法としては、とにかく推定するパラメータ数を減らすことです。例えば、Tree priorとしてBirth-Death modelを使っているならYule modelへ変更したり、インプット配列のパーティション数が多いのならばPartitionFinderなどを用いてパーティション数を減らす、また、塩基置換モデルとしてGTRなどの複雑なモデルを使用している場合、よりシンプルなHKYモデルにするなどすることで収束しやすくなります。

BEAST2による分岐年代推定①

Bayesian Evolutionary Analysis by Sampling Trees (BEAST) 2 はベイズ解析を行うプログラムです。系統解析や分岐年代推定、有効集団サイズの推定などさまざまな機能がありますが、ここでは分岐年代推定の手順を説明します。なお、以下の説明は全てバージョン2.5.2で行っています。他のバージョンだと異なる場合があります。

インプットファイルの作成

まず最初に解析に使用する配列ファイル（NEXUSファイル、FASTAファイル）を準備します。このファイルをBEAUTiという同梱ソフトに読み込ませ、xml形式のインプットファイルを作成します。

配列ファイルが準備できたらBEAUTiを起動します（図1）。

図1

配列ファイルを読み込むには左上の「File」 > 「Load」の順にクリックするか、左下の「+」マークをクリックします。クリックするとファイルの選択画面が表示されるので、解析に用いるファイルを選択します。配列ファイルがFASTA形式の場合、アミノ酸配列か塩基配列かを聞かれますので、適当な方を選択します（今回の例では塩基です）。ファイルが複数ある場合、上記の作業を繰り返して全てのファイルを読み込んで下さい。読み込み完了後、配列が複数ある場合は、全ての配列を選択した上で、「Link Clock Models」, 「Link Trees」をクリックします。

図2

続いて、「Site Model」のタブをクリックします（図2）。ここでは塩基置換モデルの設定を行います。図2はGTR+G+Iモデルでの設定例です。他のモデルを使用する場合は、適宜変更を加えて下さい。まず、「Gamma Category Count」を「4」にします。すると「Shape」という項目が現れますので、そこの「estimate」チェックボックスにチェックを入れます。「Subst Model」の項目で「GTR」を選択し、「Frequencies」を「Empirical」に設定します。最後に、「Substitution Rate」の「estimate」チェックボックスにチェックを入れます。複数の配列がある場合、これを全てに対して行います。

図3

次に、「Clock Model」タブに移動します（図3）。ここでの選択肢は4つです。「Strict Clock」は、系統樹上の全ての枝で塩基置換速度が一定であることを仮定しています。「Relaxed Clock Exponential」および「Relaxed Clock Log Normal」では、系統樹上の枝ごとの塩基置換速度はそれぞれ独立に指数分布および対数正規分布に従うとそれぞれ仮定します。「Random Local Clock」は、ベイズ法により系統樹上で塩基置換速度が一定である部分系統樹を提案し、部分系統樹ごとに塩基置換速度を推定するようです（間違っていたら教えて下さい）。そのため、塩基置換速度の変動数としては、「Strict Clock（変動は全くなし）」より大きく、「Relaxed Clock（全ての枝で変動）」より小さくなります。自分は一度使ったことがあったのですが、全く収束しなかったので諦めた覚えがあります。
　どれが最適なモデルかは正直わかりませんが、近縁種間の分岐年代推定でなければ、たいてい 「Relaxed Clock Log Normal」が使われることが多いと思います。

図4-1

「Priors」タブに移動します（図4-1）。ここでは主に較正点（Calibration point）の設定を行います。一番下にある「+」マークをクリックします。次に表示された画面では「MRCA prior」を選択します。一番上の「Taxon set label」には較正点の名前（任意）を入力します。続いて較正点を入れるノードに含まれる全てのタクサを右のボックスに移動させます。例えば、図4-2のような系統樹があり、ノード1を較正点としたい場合はC、D、Eのタクサを全て右のボックスへ移してください。タクサの移動が完了したら「OK」をクリックします。

図4-2

　「OK」をクリックすると、先程設定した較正点の名前がついた新しい項目が追加されます。ここで較正点として使用する分布を選択します。分布を選択したら左端の「▶」をクリックして分布の数値（平均、分散、上限、下限など）を入力します。また、較正点があるノードの単系統性が確実であるならば、「monophyly」にチェックを入れておきます。これにより、系統推定の間このノードが固定されます（系統推定を行わない場合は必要ないと思います）。較正点が複数ある場合、全てに対してこの作業を繰り返します。

　最後に「MCMC」 タブに移動します。「Chain Length」では、MCMC（マルコフ連鎖モンテカルロ法）の試行回数を指定します。BEASTを使っている論文では、少なくとも1000万回は行われているかと思います。足りない場合は後から追加することができるので、とりあえず1000万回から始めてみれば良いと思います。次に、「trace log」「tree log」のファイル名を入力します。これらが終わったら「File」>「Save as」を選択し、xml形式で保存してください。

初期系統樹の指定方法

図5

ここでは初期系統樹の指定方法を説明します。必要ない場合は次のステップに飛んでください。最近は初期系統樹の指定もBEAUTiを使ってできるようになりました。画面の上の方にある「View」をクリックし、「Show Starting tree panel」を選択してください。「Starting tree」というタブが新しくできるので移動します（図5）。移動したら「Initial Tree」で「Newick Tree」を選択し、その下の「Newick」ボックスに系統樹を貼り付けてください。ここで、「Adjust Tip Heights」と「Adjust Tree Node Heights」にチェックが入っていることを確認してください。万が一、入ってなかったら付けておきます。

図6

系統樹の指定は以上ですが、 次にBEAST解析の間、系統推定を行わせないようにする方法の説明をします。「View」>「Show Operators panel」を選択して、そのタブに移動してください。ここで、「Subtree Slide」と「Exchange: Tree.t ~ Narrow exchange」、「Exchange: Tree.t ~ Wide exchange」、「Wilson Balding」の4つの項目の右側の数値をゼロに書き換えます（図6）。これにより、分岐年代推定中は樹形が固定されます。

解析の実行

　解析の実行を実行するためBEASTを起動します。「Choose File」をクリックし、上で作成したxmlファイルを選択します。その後、一番下にある「Run」をクリックすると解析が始まります。BEAGLEをインストールしている場合、「Use BEAGLE Library if available」にチェックを入れておくと解析が速くなります。BEAGLEが使えるかどうかは、「Use BEAGLE Library if available」と「Show list of ...」にチェックを入れて「Run」をクリックするとわかります。