Skip to content

系統解析

MEGAによる最尤法系統樹推定

MEGAは視覚的にわかりやすいGUI形式ですので,

特に解説は必要ないと思いますが,備忘録的に.



アライメントした配列を用意します.COI配列です.

Data→Phylogenetic Analysisを選択.



アミノ酸のコーディング領域かどうかを聞かれます.

今回はCOIなのでYesです.



うまくいけばこのようなウィンドウが表示されます.

MEGAの基本ウィンドウ(ここではMEGA 6.06(6140226))でPhylogeny→Construct/Test Maximum Likelihood Tree…を選択.



ここで解析のパラメーターを設定します.



Phylogeny Test :どのように系統を検定するかです.各枝の信頼性が得たい場合はここでBootstrap Methodを選択しましょう.

No. of Bootstrap Replications:Bootstrap検定を何回行うかです.100回でも構いませんが,1000回くらいはしたほうが無難です.2-3000回やってる論文もたまに見かけます.

Substitution Type:塩基(Nucleotide)かアミノ酸(Amino Acid)が選べます.選べないと困ります.

Model/Method:塩基(アミノ酸)置換モデルが選べます.

Rates among Sites:座位ごとの置換の頻度の違いを,ガンマ分布に基づいて分類するか否かを設定できます.

No of Discrete Gamma Categories:ガンマカテゴリ数を設定できます.

Gaps/Missing Data Treatment:ギャップの扱いの設定です.Complete deletionとすると,一つでもギャップがあるサイト(列)を解析からのぞけます.

ML Heuristic Method:最尤法系統樹の探索方法を選べます.系統樹探索方法についてはまたどこかで書くかもしれませんが,局所解がたくさんあるようなデータセットだと,ここでの方法選びは結構重要だと思います.いくつか試してみても良いかもしれません.

Initial Tree for ML:系統樹探索を行う際の,最初の出発点の系統樹の作成法を選びます.上記した通り,これも結構重要だと思います.デフォルトのNJでも問題はないと思いますが,もし既にそれらしい系統樹があるのであれば,それを設定する事もできます.

Blanch Swap Filter:系統樹探索の際の枝の入れ替えの大胆さを決めます.Strongにすると枝長の入れ替えがより消極的になり,解析時間は短くなりますが,考慮する系統樹は少なくなります.よりWeakにすると,枝長の入れ替えが大胆になり,解析時間は長くなりますが,考慮する系統樹が多くなるようです.



パラメーター設定が終わったら,Computeをクリックして解析開始です.

Progressが100%になるまで,気長に待ちましょう.



解析が終わると,このようにTree Explorerに系統樹が表示されます.

各枝の上にブートストラップ確率(2桁の整数)と枝の下に枝長(有理数)が示されています.



このExplorerで色々系統樹をいじれます.例えば特定の枝を選択して,左上のコマンドの中からPlace root on Branchを選ぶと



外群を指定できます.



他にもCompress/Expand Subtreeで,枝を一つにまとめたり,



Fit Tree to Screenで,ウィンドウ内に系統を収めたり,



Flip Subtreeで枝の上下を入れ替えたりできます.

良い時代になりました.



この他にもたくさんのコマンドがありますが,感覚的に理解できると思いますので,いろいろ試してみてください.



今日はここまで.


iqtreeによる最尤法系統樹推定

※もっと良いやり方をご存知の方はそっと教えてください.

iq treeは,パーティション分けが可能な最尤法系統樹推定ソフトです.

コマンドプロンプト対応で非常に使い勝手が良く,しかも解析速度がべらぼうに速いのが特徴です.

RAxMLがコマンドプロンプトに対応した今,やり方自体はRAxMLとほぼ同じです.
 

塩基配列データ作成(MEGA, SeaVeiwなど)
  • 事前準備として解析用フォルダをつくりましょう(iqとします).

アライメントソフトでの作業

  • アライメントを行ったセッションを.phy形式で出力する(iq.phyとしましょう).
  • パーティション分けとそれぞれのモデル指定のファイルを作る(iq.nexとしましょう).以下は例と解説(赤字)です.

#nexus

begin sets;

このコマンドで読み込みを開始します.

 

charset part1 = 1-430 431-1940\3 432-1940\3;

charset part2 = 433-1940\3;

charset part3 = 1941-2889;

ここで,モデルごとのパーティション分けを指定します.

ハイフンでつないだ配列が一つの領域で,それぞれをスペースで分けます.

アミノ酸指定領域の場合は,”日本円マーク”か”\”で各コドンをします.

ここでは,431番目から始まるアミノ酸指定の第一コドンと第二コドンは同じパーティションですが,

第三コドンは違うパーティションとしています.

 

charpartition mine = GTR+I+G:part1, TN93+I+G:part2, K2P+G: part3;

ここで,各パーティションごとのモデルを指定します.

よくモデルテストで見る形式をそのまま入力すればよいので楽です.

主な塩基置換モデルのリストはコチラ

 

end;

おしまい.


iqtreeを走らせる(コマンドプロンプト,iqtree)
  • 事前準備として,iqフォルダにiqtreeからDLしたiqtree.exeとiqtree-click.exeを入れておきます.

コマンドプロンプトでの作業

  • コマンドプロンプトを立ち上げ,iqフォルダまでのパスを通します.
  • iqtreeを走らせるためのコマンドを入力します.以下,例とコマンドの解説(赤字)です.

iqtree -s iq.phy -spp iq.nex -m TEST -bb 1000

“-s”で配列ファイルの読み込みを行います.

“-spp”でパーティション分けファイルの読み込み,

“-m” で検定方法を指定,

“-bb”でブートストラップ検定とその回数の指定を行います.


 

  • 解析が終わったら,iqフォルダ内に解析のログファイルや,系統樹ファイルiq.phy.treefileが生成されるはずです.
  • 参考までに,iqtreeでは他にも様々なコマンドがあります.コチラをご参照ください.

 

今日はここまで.


Expanded Bayesian Skyline Protによる集団の遷移年代推定

※もっと良いやり方をご存知の方はそっと教えてください.

細かいコマンドは勉強中です.こちらもいろいろ教えてください.

 

塩基配列データ作成(MEGA, BEAUti)

MEGAでの作業

  • 解析用フォルダを作る(Skyとしましょう).
  • 解析したい内群について,MEGAでアライメントを行ったセッションをfasファイルで出力する(Sky.fasとしましょう).
各種パラメータ設定(BEAUti)

BEAUtiでの作業

  • File→Import AlinmentでSky.fasを読み込む.Partitionタブでアライメントセッションファイルの名前が確認できればOK.

  • Site ModelタブのSubst Modelでモデルを設定する.

  • Clock Timeタブで,特にこだわりが無ければ”Strict Clock”を選び,Clock rateに塩基置換レートを入力する.例えば,1.25%/100万年であれば,0.0125と入力.

  • PriorsタブでTree.t: EBSP(アライメントのセッション名)から”Coalescent Extended Bayesian Skyline”を選ぶ.

  • Tree.t: EBSPの左の▶をクリックし,Population modelの横の鉛筆をクリックし,Factorを入力する.母系遺伝のミトコンドリアは.0.5と入力するのが良いようです.

  • MCMCタブのChain lengthで,MCMC連鎖の数を設定する.デフォルトは10,000,000ですが,多分30,000,000位はやったほうがよさそうです.この後のTracerで結果を見ながら,この値を調整するとよいでしょう.

  • Fileの並びのメニューのVeiwタブから,Show Operators panelを選ぶ.

  • もし複数の遺伝子データを使う場合は,Operatorsタブで,Bit Flip…, Sample off…, Scale…の値をそれぞれの数×30, 15, 15にするとよいようです.今回はデフォルトですが,もし二つの遺伝子を使う場合は60, 30, 30にするとよいでしょう.
  • これで準備OKです.File→Saveでxmlファイルを保存しましょう(Sky.xmlとしましょう).
MCMCによる解析(BEAST)

BEASTでの作業

  • BEASTを立ち上げ,Choose FileでSky.xmlを選ぶ.後はデフォルトでOKです.Runをクリックして解析を始めましょう.

  • 上手くいけば,カリカリ解析が始まるはずです.苦労したファイルが走るのを見るとなんか感動しますよね.
  • 解析が終了して,SkyフォルダにSky.logとEBSP.logが生成されればOKです.
集団推移の可視化(Tracer, R)

Tracerでの作業

  • Tracerを立ち上げ,Import Trace FileでSky.logを選ぶと,ファイルが読み込まれて,Traces:のところにlogの解析データが出力されます.

  • このうち,ESSの値が重要なようで,これが200以下の項目がある場合はMCMCの連鎖回数を増やすとよいでしょう.

  • 一番下のsum(indicators.alltrees)をクリックし,右側のEstimatesのタブをクリックすると,ヒストグラムが表示されます.

  • Rを立ち上げ,Skyフォルダ内に移動し,以下のように入力すれば,少し計算時間があった後,集団の遷移図がプロットされます.

> source (“plotEBSP.R”)

> plotEBSP(“EBSP.log”, useHPD=FALSE, log=”y”)

 

  • でました.横軸が年代(百万年)で,縦軸が集団サイズです.左が現在ですので,この集団はごく最近に少しだけ集団サイズの減少を経験しているようです.
  • x軸の表示範囲を調整するには,Rに以下のように(0~0.03の範囲で表示する場合)入力します.

> plotEBSP(“EBSP.log”, useHPD=FALSE, xlim=c(0, 0.03))

今日はここまで.


コマンドプロンプトの使い方

非常に基本的なソフトなので,使い方というほどでもないのですが,

最近はこれで動かせる系統解析ソフトも多いのでご紹介までに.



起動方法①

  • 画面左下のウィンドウズマークをクリック.
  • 画面左下の下矢印をクリック
  • アプリの一覧から「コマンドプロンプト」を探し,クリック

起動方法②

  • 検索画面(画面右下に矢印を持っていって,虫眼鏡マークをクリック)を開き,「cmd.exe」と入力.
  • Enterキーを押す.


  • このような画面が現れたらOK.ちなみに私はタスクバーにピン止めしています.


基本的なコマンド

  • コマンドプロンプトでは,まず解析したいフォルダに移動します.いわゆるパスを通す,というやつですが,そのコマンドが”cd”です.”>”の後に”cd”と入力し,その後に移動したいフォルダのパスを入力すれば,そのフォルダに移動できます.
  • 簡単なやり方は,フォルダの検索窓の横のところにカーソルを持っていき,クリックして下の画面のように青くします.


  • この状態でコピーし,コマンドプロンプトの画面で右クリックをすると,
  • このように比較的簡単にパスがペーストできます.後はこのパスの先頭部分に”cd “を入れればOKです.指定したフォルダに移動できます.


他にも以下のコマンドを知っておいて損はないでしょう.

“dir”:ディレクトリ内の全てのフォルダやファイルを表示.

“tree”:ディレクトリ構造をツリー状に表示.

“ren”:ファイル名を変更.フォルダ内の拡張子を一気に変更する際に便利です.



今日はここまで.


ArelquinによるAMOVA解析

※もっと良いやり方をご存知の方はそっと教えてください.

参考

http://hanzawalab.blog.fc2.com/blog-entry-43.html

http://nonomasu.blog.fc2.com/blog-entry-7.html

塩基配列データ作成(MEGA, DNAsP)

MEGAでの作業

  • MEGAでアライメントを行ったセッションをnexusファイルで出力する(Amo.nexusとしましょう).この段階で,AMOVA解析で扱う集団ごとにまとまるように配列の順番を入れ替えておくと後が楽.

DnaSPでの作業

  • 解析用のフォルダ(Amoとしましょう)を作り,その中にAmo.nexusを入れておく.

  • DnaSPを立ち上げ,File→Open Data FileでAmo.nexusを読み込む.
     

 

  • Data Informationが表示されれば,読み込み成功.Data InformationをCloseで閉じる.
     

 

  • Data→Formatでデータの形式を指定する.例えばミトコンドリアの16Sの場合は,DNA, Haploid, Mitochondorialにチェックを入れてOKをクリック.


 

  • Data→Define Sequence SetsでOTUのグループ分けを指定する.List of All Sequenesから,”>>”でIncluded  ListにOTUを移動し,Add new Sequence Setで名前を付ける.
  • 全てのグループに名前を付けたら,Update All Entries(赤文字になっているはず)をクリック.
  • Save/Export Data as…→Arlequin File Formatを選ぶ.
     

  • Not considered, Included,  Arlequin Haplotype Listにチェックが入っている事を確認し,OK.
  • .hapと.arpの二つを保存する.Amo.hap, Amo.arpとしましょう..hapは.arpのバッチファイル的なものらしく,必ずペアでAmoに保存しておく.
  • Amova解析の際には各集団を更に上位の集団にまとめる必要があるので,以下のテキストをAmo.arpの末尾に貼り付ける.赤字は任意に変えられる部分.

[[Structure]]

StructureName=”任意の名前
NbGroups=4

Group={
DnaSPで定義したグループ名1
}

Group={
DnaSPで定義したグループ名2
DnaSPで定義したグループ名3
}

Group={
DnaSPで定義したグループ名4

Group={
DnaSPで定義したグループ名5
DnaSPで定義したグループ名6
}


  • これでArlequin用の準備完了 .

※Arelquin上でもっと簡単にグループ設定する方法がありました↓

AMOVA解析(Arlequin)
  • Arelquinを起動し,”Open project”でAmo.arpを選択.左側の枠にSamplesやGroupsが表示されたらOK.

  • グループを設定する場合は,Structure Editorをクリックし,Population samplesの横のGroup(デフォルトでは0が表示されている)をダブルクリックし,任意のグループ名を入れると.右側の枠に階層構造として出力される.最後はUpdate ProjectをクリックすればOK.Projectタブでグループ分けが出来たかどうか確認できる.
  • “Settings”でAMOVAを選び,Standard AMOVA…をチェックし,プルダウンメニューの下の方で塩基置換モデルを選べる.普通はKimura 2Pが無難.
  • “Start”をクリックして解析.うまくいけば,Amo内にログなどが入ったAmo.resというフォルダが作成される.この中のAmo.xmlにAMOVA解析結果がHTML形式で保存されている.

更新:2017/3/11


RAxMLによる最尤法系統樹推定

※もっと良いやり方をご存知の方はそっと教えてください.

塩基配列データ作成(MEGA, SeaView)

MEGAでの作業

  • MEGAでアライメントを行ったセッションを出力する.配列の結合を行う場合は,その後SeaViewで読み込むので,”.nexus” or “.nex”でも”.fasta” or “.fas”とかでもOK.行わない場合は”.phylip or “.phy”で出力.
  • ↑のファイルのOTU名は,それぞれのファイルで同じにしておく(”Ophi281_COI”とか”Ophi281_16S”じゃダメ).また,配列内のギャップ(”-“で指定する事が多い)と混同する事が多いため,OTU名には”-“は使わないほうがいい.

SeaViewでの作業(配列を結合する場合)

  • FileタブのOpen ***(***はファイルのフォーマット)を選び,結合させたい配列その1を選ぶ.
  • 同じようにその2,その3を選ぶ.それぞれ別ウィンドウで表示されればOKだが,不安な場合はFileタブのNew windowを選び,新しいウィンドウで配列を読み込めばOK.
  • FileタブのConcatenateを選ぶと,結合先の配列選ぶウィンドウが表示されるので,結合したい配列を選ぶ.by nameを選び,”OK”をクリック.
  • 上手くいけば,各OTUの後に配列が付加される.
  • これで準備完了.最終的な結合配列をFileタブのSave asからphylip形式で保存する(RA.phyとしましょう).
最尤法系統樹作成(RAxML)
  • 解析用の適当なフォルダ(RAとしましょう)を作り,その中にRA.phyを入れておく.
  • RA内にパーティション分け用のファイルを作っておく.例えば1-300 bpまでが16S, 301-1200までがCOIの場合は,以下のようなテキストデータファイルを作ります(partition.txtとしましょう).

DNA, 16S=1-300
DNA, COI=301-1200


  • RA内に”raxmlHPC.exe”を入れておく.これがRAxMLの起動ファイルです.これで準備完了.
  • コマンドプロンプトを起動(ウィンドウズメニューから探すか,ファイル検索で”cmd.exe”を検索)し,”cd 「RAのアドレス」”を入力してRA内に移動します.「RAのアドレス」はフォルダのアドレスバーの部分を左クリックすると”C:\Users\****\Desktop\研究\分子系統解析\RA”みたいなのが選択できるので,コピーしてコマンドプロンプト内で右クリックするとペーストできます.
  • “>”の後ろにRAxMLの起動コマンドを入力します.このコマンドは色々あるのですが,例えば

raxmlHPC -f a -m GTRGAMMA -p 12345 -x 12345 -# 1000 -s RA.phy -q partition.txt -n out

  • と入力すれば解析が始まり,GTRGAMMAモデルでブートストラップ1000回行った最尤法系統樹が得られます.
  • このコマンドなどについては井上潤さんのページにすごく詳しいです.
  • 出力されたファイルのうち,”RAxML_bipartitions.out”が系統樹のファイルです.TreeViewやFigtreeなどで閲覧して見てみましょう.

更新:2017年3月5日


ハプロタイプネットワーク構築

※もっと良いやり方をご存知の方はそっと教えてください.

塩基配列データ作成(MEGA, DnaSP)

MEGAでの作業

  • MEGAでアライメントを行ったセッションを,”.nexus”で出力する.

DnaSPでの作業

  • 出力したファイルをDnaSPで読み込む(File→Open Data File).
  • Data Informationウィンドウが表示されたら読み込み完了.Closeでウィンドウを閉じる.この時,Nは読み込むが,その他の記号(YとかWとか)は読み込まないようなので,おとなしくNに変換するか,その配列もしくはOTUは削除する.
  • Generate→Haplotype Data Fileで,Haplotype/DNA… ウィンドウを表示させる.GenerateのRoehl Data File (Network software)をチェックし,.rdfファイルを出力する.
  • この時Output of:… というウィンドウが表示されるが,これはハプロタイプとOTU名の対応が書かれた重要なデータなので,メモ帳にコピペなどをしておく(Haplotype/DNA… ウィンドウでNEXUS Haplotype Data Fileにチェックを入れば,.nexとして出力できる).
ハプロタイプネットワーク作成(Network)

Network での作業

  • Data Entry→Import rdf fileで,.rdfファイルを読み込む.
  • File selection ウィンドウでDNA Nucleotide dataがチェックされている事を確認し,Continueをクリック→.rdfファイルを開く.
  • RDF-Editorウィンドウで各種設定を確認・変更可能.設定が終わったらSaveをクリックし,元の.rdfファイルを上書き,あるいは別名で保存.
  • RDF-EditorウィンドウをExitで閉じ,Calculate Network→Network Calculations→Median JoiningでMedian Joiningウィンドウが開く.
  • Median JoiningウィンドウでFile→Openで.rdfファイルを再び開く.
  • Median JoiningウィンドウでCalculate networkをクリックすると計算が開始され,.outファイルが出力される.
  • メインウィンドウのDraw networkでDraw Networkウィンドウを表示させ,File→Openで/outファイルを開く.
  • OK→Yes→Continue→Finaliseでハプロタイプネットワークが完成.各ハプロタイプを右クリックすると,円グラフの組成がいじれる.

 

更新:2017年3月4日