Skip to content

卒業式

の季節ですね.

茨城大学でも本日,勉学を修めた若人が巣立っていきました.

茨大では半分以上が院に進みますが,岡西研では全員が卒業です.



式後に研究室に来た大江君とパシャリ.

一年間,クモヒトデを研究してくれて,ありがとうございました.

東京へ行っても元気でね!



がらんとした研究室を見ているとなんだか感傷的になってしまいますが,

また新入生が参入してくる四月は目の前です.



新たな日々に向けて,頑張りましょう!


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


MEGAによるモデルテスト

MEGAでは塩基(アミノ酸)置換の進化モデルテストを行うことができます.

進化モデルとは,解析するデータセットの中における,



①塩基(アミノ酸)の存在頻度

②塩基(アミノ酸)同士の置換確率



の組み合わせです.

最尤法では,節間の置換を計算するので,この進化モデルの選択は非常に重要です.



まずはアライメントした配列を用意します.今回はCOI領域の配列です.

Data→Phylogenetic Analysisを選択します.



タンパク質のコーディング領域かどうかを聞かれます.

今回はCOIなのでYesです.



このようなウィンドウが現れたらOKです.アライメントの中から,塩基置換が起きている部分だけを取り出したものです.

このPhylogenetic analysisモードは,MEGAで系統解析を行うための基本モードです.



Models→Find Best DNa/Protain Models (ML)を選びます.



この画面でモデルテストのセッティングができます.
 

Tree of use :テストの際に使う系統樹を指定できます.通常,デフォルトのNJ法でOKです.

Substitution type:Nucleotide(塩基)やAmino acid(アミノ酸)が選べます.

Gaps/Missing Data Treatment:ギャップやミッシングデータの取り扱いが選べます.Complete deletionで,ギャップが一つでもある列は解析殻除きます.Partial Deletionで,比較する配列ごとにギャップの扱いを決めます.Use all siteでギャップもミッシングデータも全て扱います.

Select Codon Positions:チェックを外したコドンを解析から除きます.コーディング領域の場合,特に3rdコドンが他とモデルが異なる場合があったりするので,きちんと全てのコドンについてのモデルテストを行っておくことをお薦めします.

Branch Swap Filter:系統樹の枝長に対する厳密性を決めるオプションのようです.複雑なモデルを使うと,色んな塩基置換パターンを考慮に入れられる反面,その分の解析上の煩雑さも増えてしまいます.この煩雑さを枝長と考え,各モデルごとに枝長と尤度を比較し,尤度の上昇と枝長の少なさのバランスが最も良いものをベストのモデルとして選ぶようです.この枝長は系統樹の探索によって決められていくのですが,その際に系統樹の一部の枝を入れ替えます.この時の入れ替えの「大胆さ」を決めるのがこのオプションのようです.従って,Strongにすると枝長の入れ替えがより消極的になり,解析時間は短くなりますが,考慮する系統樹は少なくなります.よりWeakにすると,枝長の入れ替えが大胆になり,解析時間は長くなりますが,考慮する系統樹が多くなるようです.ですので,時間に余裕があるときはよりWeakにするとよいでしょう.

※個人的な見解ですので,間違っていたら教えてください.
 

さて,長くなりましたが,セッティングを終えてComputeを選択するとモデルテストが始まります.



モデルを一つ一つ試しています.終わるまで気長に待ちましょう.



結果が出ました.注目すべきは,BIC, AICcです.これらの値が最も少ない(=尤度と煩雑さのバランスが最も良い)ものが,ベストなモデルとなります.

今回の解析では,GTR+G+Iがベストなモデルとして選択されました.他にも,f(A),f(T)…などはそのモデルの下での塩基の存在頻度で,r(AT)…などが置換確率となります.



このウィンドウの下の方に,それぞれのモデルに対する説明があります.

ここを参考にして,今後の最尤法解析などの際にパラメーターを設定します.



今日はここまで.


MEGAによるアライメント②

アミノ酸配列をコードしている領域の場合は,

DNA配列をアミノ酸配列に変換することで,

正しく配列が得られているか確認できます.
 

こちらは,アミノ酸をコードしているミトコンドリアのCOI領域

のDNA配列をアライメントしたものです.



ちなみに,このOTUの名前の付け方は長すぎる上に,

()とか系統解析状優しくない記号が使われている悪い例です.



各配列の一番上の行にアスタリスク(*)が示してある列は,

全てのOTUで塩基置換が同じ事を示しています.

従って,*が無い列はどこかのOTUで塩基置換が起きています.

この*を見ると,二列おきに無くなっている傾向が見られると思います.

これは,この配列がアミノ酸をコードするコドンは塩基3つで一組となっており,

その三番目が置換しやすい事を表しています.



左上の方のDNA Sequencesの横のTranslated Protein Sequencesのタブをクリックすると,



このようなウィンドウがポップアップするので,

Yesを選択するとDNA配列がアミノ酸配列に変換されます.

Noを選択すると,
 

このように,コドンがコードするアミノ酸の遺伝コードを選択する事ができます.

分類群に合わせて変更しましょう.



遺伝子コード表は,Data→Select Genetic Code Tableからも選択できます.
 

これが置換後のアミノ酸配列です.色付きの文字は全てアミノ酸です.

灰色の*は,ストップコドンでこれが見られる場合は,

アミノ酸変換が上手くいっていない事がほとんどです

(勿論,複数の遺伝子が発現するためのストップコドンの場合もあります).



というわけで,一番左の列を削除してもう一度アミノ酸変換してみました.

うーん,まだストップコドンがみられます.
 

更にもう一列削除してアミノ酸変換しました.

ストップコドンもみられず,配列も結構「揃って」います.

恐らく,これで正しくアミノ酸配列に変換できたものと思われます.
 

今日はここまで.


MEGAによるアライメント

MEGA (Molecular Evolutionary Genetics Analysis software)は,

アライメント,モデルテスト,系統解析までを一手にこなせる優れた系統解析ソフトウェアです.

昔はアライメントや簡単な系統解析が出来るくらいでしたが,現在はver. 7までアップデートを重ね,

モデルテスト,最尤法,祖先解析,分子分岐年代測定などが行えるようになっています.

また,開発者が日本人という事で,日本語でも扱えるため,大変有用なツールとなっています.

DLはこちら
 

備忘録的に,このMEGAの使い方を記録していきます.

MEGAをインストールしたら,MEGAのアイコンをダブルクリックすると,以下のウィンドウが起動します.
 

これがMEGAの基本画面です.
 

まずアライメントです.”Align”→”Edit/Built Alignment”を選択.
 

Creat a new alignmentをチェックし,OKをクリックすると,以下のようなExplorerが起動します.
 

“Edit”→”Insert Sequence From File”を選択します.
 

DNA(もしくはアミノ酸)配列のファイルをPCから読み込みます.読み込める形式は色々ありますが,

このように,メモ帳に直接塩基配列をペーストしただけのものでもOKです.
 

配列が読み込めました.色が付いている部分が塩基配列です.A, T, G, Cでそれぞれ色分けされています.
 

“Alignment”から,”Clustal W”と”Muscle”の二つのアライメント方式が選べます.

今回はClustal Wを選択してみます.
 

ここでは,アライメントのパラメータの設定ができます.

Gap Opening PenaltyやGap Extension Penaltyの値を高く設定すれば,なるべくギャップが少なく,連続しないアライメントになります.

DNA Weight Matrixでは,塩基とアミノ酸の置換に対するスコアを選択できます.

Transition Weightで,転位と転換の差を設定できます.

他にも,Keep Predefined Gapsでギャップを完全に無視したアライメントや,

Specify Guide Treeで,アライメントの指標となる系統樹を指定できるようです.
 

もしデフォルトのセッティングでアライメントが上手くいかない場合は,この値などをいじってみましょう.
 

とりあえずデフォルトの設定でOKを押すと,アライメントが始まります.
 

アライメントができました.うまく揃っている部分と,ごちゃごちゃになっている部分があります.

これはリボソームRNAの配列なので,立体構造をとったときの活性部位と不活性部位です.

ごちゃごちゃの部分はギャップが多いエリアを除いてしまうか,

立体構造を考慮して系統解析を行います.
 

今日はここまで.