Note on HyPhy

2007 年 9 月 17 日　改訂
井上 潤
個人的なメモを書いています．HyPhy に関するコメントはこちらをご覧下さい．
アミノ酸 rate matrix の最尤推定
hp_condor.txt
petal では動作を確認．phoenix は動いたか?
アミノ酸 rate matrix の最尤推定
セルゲイ version
[inoue@petal ModelFitter]$ pwd
/home/inoue/HYPHY_Source/ModelFitter

mtREV_FitREV.bf: mtREV で樹長を最尤推定し，R matrix とアミノ酸組成 (Pi) を計算．
ExCyb57Act.nex: シーケンスファイル (tree 入り)．
Act.list: シーケンスファイルの所在を記載．
Act57: ExCyb57Act.nex が入っているディレクトリ．

たまに動かなくなる．Act57 をマック側と交換するだけで動いた．
/Users/inoue/HYPHY_Source/FromSergei_ModelFitter/QmatrixEstimation

[inoue1:FromSergei_ModelFitter/QmatrixEstimation/TrCybResults] inoue% pwd
/Users/inoue/HYPHY_Source/FromSergei_ModelFitter/QmatrixEstimation/TrCybResults
1RPi_Q.pl: R*(Pi) を行い，Q matrix を計算．
自家製
[inoue@petal ~/Hyphy_G8_petal_fol]$ pwd
/home/inoue/Hyphy_G8_petal_fol

EstAAMatrix_G8.txt
 mtmam の rate matrix を推定するために (下のコラム)，paml4/dat/mtmam.dat に従い，REVaa+dGamma(K=8) with the estimated gamma parameter to be 0.37 を設定してある．

mtCDNAmamAA.txt
 paml/example/mtCDNA/mtCDNAmam.nuc をアミノ酸に変換し nexus 形式で保存．tree の樹長は mtCDNAmam.trees の topology を K=8, alpha = 0.37 として codeml で計算した．
codeml と HyPhy で推定された rate matrix の比較
[inoue@petal ~/Hyphy_G8_petal_fol]$ pwd
/home/inoue/Hyphy_G8_petal_fol

上記 mtCDNAma のデータを codeml と HyPhy (バッチファイル: EstAAMatrix_G8.txt) で解析し，結果を mtmam.dat と比較した．

HyPhy では tree に樹長付きと樹形のみのものを用いた．rate matrix を推定した際に得られた系統樹は，樹長も含めて，codeml, HyPhy (樹長付き prior)，HyPhy (樹長なし prior) の間で差異は見られなかった．

さらにR も Pi も大きな違いが見られなかったため，HyPhy と自作の EstAAMatrix_G8.txt は，正常に作動していると思う．さらに，starting tree に樹長を含めてもそうでなくても，codeml とほぼ同じ rate matrix が得られることがわかった．
アミノ酸サイト別の尤度を推定
mtREV_Test.bf　バッチファイルで #include "mtREV_24_F.mdl" を行っている．
IncWAGF_Test.br　バッチファイル内部に WAGF matrix を書き込んでいる．

Model options では Fixed rates を選んでいる．variation や variation+HM モデルは目盛りを増やさないと動かないというエラーメッセージが出る．57 OTU, 350 AA, topology 固定でも 90MB メモリーを増やすように言われた．
[inoue@petal AAsiteLnL_fol]$ pwd
/home/inoue/HYPHY_Source/AAsiteLnL_fol
[inoue@petal AAsiteLnL_fol]$ cat mtREV_Test.bf

DataSet ds = ReadDataFile ("TestActCYB.nex");
DataSetFilter filteredData = CreateFilter (ds,1);

#include "mtREV_24_F.mdl";

Tree givenTree =((ChiChiMon,(ScyScyCan,RajRajRad)),(((PolPolSpa,PolPseGla),
(AciHusHus,(AciAciSte,AciAciTra))),(AlbAlbGlo,(AlbPteGis,NotNotChe))));

LikelihoodFunction lf = (filteredData,givenTree);
Optimize (res,lf);

fprintf (stdout, "\nOptimization results:\n", lf, "\n");

/* get a vector of site by site likelihoods*/

ConstructCategoryMatrix (siteBySite, lf, COMPLETE);

/* if there is more than 1 row, then there was rate variation in the model - so we ignore it for now.
I'll expand this example to show you how this can be extended */

if (Rows (siteBySite) == 1)
{
fprintf (stdout, "\nSite-by-site likelihoods:\n\n");
for (site = 0; site < Columns (siteBySite); site = site+1)
{
/* Format commands are simply to make the output align in columns nicely */

fprintf (stdout, "Site ", Format (site+1,5,0), ": ", Format (Log(siteBySite[site]),15,10), "\n");
}
}

END;
MacClade で作成した Nexus 形式から OTU を削り，simple な形にする

perl -p Exclude.pl infile > outfile

　Exclude.pl のスクリプト
　------------------------------------------------------------------------------------------------

s/Parajulis_poecilepterus_NC_009459[ ][ ].*?\r//g; #
s/Lophius_litulon_NC_008125[ ][ ].*?\r//g; # Lophiiformes
s/Caulophryne_pelagica_NC_004383[ ][ ].*?\r//g;

s/#NEXUS \r\[.*?]\r\r/#NEXUS \r[ $ARGV ]\r/g;
s/NTAX=337/NTAX=213/g;
s/[ ]SYMBOLS = " 1 2 3 4"[ ][ ]//g; s/ CHARSTATELABELS.*?;\r//g;
s/MATRIX\r\[.*?10[ ]{8}20.*?\]\r\[.*?\.[ ]{9}\..*?\]\r\r/MATRIX\r/g;
s/\r\rBEGIN CODONS;.*//g;

　------------------------------------------------------------------------------------------------