rgene_gamma计算2023-10-02

rgene_gamma计算，准备如下两个文件，以及ctl配置文件

树文件：请注意，根节点需要一个简单的化石校准，比如“@4.5”（450Ma ago）。如果没有这个标记计算出来的结果就没有Substitution rate is per time unit,也就没办法计算后续。
序列文件：fas 或者 phy

---

核苷酸序列 DNA dataset

针对核苷酸序列 DNA dataset，需要使用baseml baseml.ctl

      seqfile =input_pmal.phy
      outfile = mlb        * main result file
     treefile = input_marktime.tree

        noisy = 3   * 0,1,2,3: how much rubbish on the screen
      verbose = 1   * 1: detailed output 详细输出, 0: concise output 简明输出
      runmode = 0   * 0: user tree;  1: semi-automatic;  2: automatic
                    * 3: StepwiseAddition; (4,5):PerturbationNNI 

        model = 7   * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
        Mgene = 0   * 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff

    fix_kappa = 0
        kappa = 2   * initial or given kappa

    fix_alpha = 0 
        alpha = 0.5  * initial or given alpha, 0:infinity (constant rate)
       Malpha = 0   * 1: different alpha's for genes, 0: one alpha
        ncatG = 5   * # of categories in the dG, AdG, or nparK models of rates

      fix_rho = 1  
          rho = 0.  * initial or given rho,   0:no correlation
        nparK = 0   * rate-class models. 1:rK, 2:rK&fK, 3:rK&MK(1/K), 4:rK&MK 
           
        clock = 1   * 0: no clock, unrooted tree, 1: clock, rooted tree
        nhomo = 1   * 0 & 1: homogeneous, 2: kappa's, 3: N1, 4: N2
        getSE = 1   * 0: don't want them, 1: want S.E.s of estimates
 RateAncestor = 0   * (1/0): rates (alpha>0) or ancestral states (alpha=0)
    cleandata = 0  * remove sites with ambiguity data (1:yes, 0:no)?

noisy控制在屏幕上显示的输出量,对于计算较大的数据集可以减少屏幕上输出内容，数字越大输出越少
runmode尽可能多地使用runmode=0，提供一个树结构并不是困难的事
model 模型没得说，看下图选择

image.png

7是常用GTR [REV (also known as GTR)]

Mgene与序列数据文件中的选项G组合使用，用于对来自多个基因或多个位点分区（如三个密码子位置）的数据进行组合分析。如果数据文件中未使用option G，请选择0。

option G

image.png

这里顺便说明一点，对于分歧时间估计，其本身就是一个推测时间。自己测试过很多参数的变化，对分歧时间上改变并不大，两个mcmctree的run之间差异估计跟小参数改变差不多，所以很多比较复杂的参数就选择了最简单的模式。PS：实际测试时rgene_gamma值的变化影响都不大，但是这是为数几个方便计算的参数。

ndata: 指定文件中独立数据集的数量，实际就是分区数量。
clock: clock=0意味着没有时钟，并且速率在不同分支之间完全可以自由变化；clock=1表示全局时钟，并且所有分支具有相同的速率；clock=2指定本地时钟模型，系统发育中的大多数分支符合时钟假设，并具有默认速率（r0=1）。但一些预定义的分支可能具有不同的速率，分支的速率是使用树文件中的分支标签（#和$）指定的比如(((1,2) #1, 3)；选项clock=3用于多基因或多分区数据的组合分析，允许分支速率在数据分区之间以不同的方式变化（为了说明数据分区之间的进化过程中的差异，您必须使用序列文件中的选项G以及控制文件中的控制变量Mgene）。
clock = 1 or 2，需要有根的树，如果化石信息被指定，就会计算absolute rate

fix_kappaand kappa指定K80、F84或HKY85中的κ是固定值还是根据数据迭代估算，如果fix_kappa=1，则另一个变量kappa的值为给定值，否则（fix_kapa=0），则kappa值用作迭代的初始估计。
fix_alpha and alpha 原理同上，是指跨位点可变替代率的伽马分布的形状参数α。
fix_rho and rho以类似的方式工作，涉及相邻站点速率的独立性或相关性，其中ρ（rho）是自离散伽马模型的相关参数
nparK指定了站点间可变和马尔可夫相关率的非参数模型：nparK=1或2表示站点的几个（ncatG）类独立率，而nparK=3或4表示相邻站点的率是马尔可夫相关的；nparK=1和3具有每个速率类别具有相等概率的限制，而nparK=2和4不具有
这一限制（杨1995）。变量nparK的优先级高于alpha或rho。(神仙打架，我看=0，应该就是不用这个参数)
nhomo nhomo = 1适用于均质模型（F81, F84, HKY85, T92, TN93, or REV），通过最大似然迭代估计频率参数πT、πC和πA（因为πG不是自由参数，因为频率和为1）。在默认选项nhomo=0的情况下，频率参数是通过观测频率的平均值来估计的。对于K80、F84和HKY85模型，nhomo=2对树中的每个分支使用一个转换/转换率比率（κ）

getSE是否需要估计参数的标准误差的估计，当然使用=1
RateAncestor通常使用0。=1强制程序执行两个额外的分析，如果不需要结果，可以忽略这两个分析。首先，在γ等站点的可变利率模型下，RateAncestor=1迫使程序使用经验贝叶斯程序计算序列中各个站点的利率（文件利率中的输出）。
cleandata=1意味着从所有序列中去除涉及模糊特征（未确定的核苷酸，如N、？、W、R、Y等，四个核苷酸以外的任何核苷酸）或比对间隙的位点。这样可以加快计算速度。cleaddata=0（默认值）使用这些站点。

氨基酸序列 Amino Acid Dataset

氨基酸序列 Amino Acid Dataset，需要使用codeml codeml.ctl

      seqfile = input.fas
     treefile = input_newick.tree
      outfile = mlb           * main result file name

        noisy = 9  * 0,1,2,3,9: how much rubbish on the screen
      verbose = 1  * 0: concise; 1: detailed, 2: too much
      runmode = 0  * 0: user tree;  1: semi-automatic;  2: automatic
                   * 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise

      seqtype = 2  * 1:codons; 2:AAs; 3:codons-->AAs
    CodonFreq = 2  * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table

*        ndata = 1
        clock = 1  * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
       aaDist = 0  * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a
   aaRatefile = paml_dir/dat/jones.dat  * only used for aa seqs with model=empirical(_F)
                   * dayhoff.dat, jones.dat, wag.dat, mtmam.dat, or your own

        model = 3
                   * models for codons:
                       * 0:one, 1:b, 2:2 or more dN/dS ratios for branches
                   * models for AAs or codon-translated AAs:
                       * 0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F
                       * 6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)

      NSsites = 0  * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
                   * 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
                   * 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
                   * 13:3normal>0

        icode = 0  * 0:universal code; 1:mammalian mt; 2-10:see below
        Mgene = 0
                   * codon: 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff
                   * AA: 0:rates, 1:separate

    fix_kappa = 0  * 1: kappa fixed, 0: kappa to be estimated
        kappa = 2  * initial or fixed kappa
    fix_omega = 0  * 1: omega or omega_1 fixed, 0: estimate 
        omega = .4 * initial or fixed omega, for codons or codon-based AAs

    fix_alpha = 1  * 0: estimate gamma shape parameter; 1: fix it at alpha
        alpha = 0. * initial or fixed alpha, 0:infinity (constant rate)
       Malpha = 0  * different alphas for genes
        ncatG = 5  * # of categories in dG of NSsites models

        getSE = 0  * 0: don't want them, 1: want S.E.s of estimates
 RateAncestor = 1  * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)

   Small_Diff = .5e-6
    cleandata = 1  * remove sites with ambiguity data (1:yes, 0:no)?
*  fix_blength = -1  * 0: ignore, -1: random, 1: initial, 2: fixed
       method = 0  * Optimization method 0: simultaneous; 1: one branch a time

* Genetic codes: 0:universal, 1:mammalian mt., 2:yeast mt., 3:mold mt.,
* 4: invertebrate mt., 5: ciliate nuclear, 6: echinoderm mt., 
* 7: euplotid mt., 8: alternative yeast nu. 9: ascidian mt., 
* 10: blepharisma nu.
* These codes correspond to transl_table 1 to 11 of GENEBANK.

Amino acid sequences用seqtype = 2，model，alpha, aaRatefile
model = 3指定氨基酸取代的模型：=0是假设速率相等，泊松模型，对于任何氨基酸取代（Bishop和Friday 1987）；=1用比例模型，其中氨基酸的变化率与该氨基酸的频率成比例，=2/3指定了经验模型，而经验氨基酸替换率矩阵在aaRatefile指定的文件中给出。
the empirical models（aaRatefile）：Dayhoff (dayhoff.dat), JTT(jones.dat)， Revised JTT（jones-dcmut.dat）, WAG (wag.dat), mtMAM (mtmam.dat), mtREV24 (mtREV24.dat), cpREV(cpREV10.dat、cpREV64.dat)，LG（lg.dat），mtART（mtArt.dat）
fix_alpha, alpha:使用fix_alpha=0，根据数据迭代估算。

codon sequences 用seqtype = 1，model，aaDist，CodonFreq，NSsites，icode
CodonFreq = 2CodonFreq指定了密码子替换模型中的平衡密码子频率。
假设频率相等（对于标准遗传密码，每个频率为1/61，CodonFreq=0），
根据三个密码子位置的平均核苷酸频率（CodonFreq=2）
自由参数（CodonFreq=3）。其中涉及的参数数量密码子频率的模型为0、3、9和60（对于通用代码），对于CodonFreq=0，分别为1、2和3。
aaDist = 0指定假设氨基酸距离相等（=0）还是使用Grantham's matrix（=1），在名为OmegaAA.dat的文件中指定了取代类型的数量以及哪对氨基酸变化属于哪种类型。您可以使用该模型来拟合“保守”和“带电”氨基酸取代的不同ω比率。
NSsites = 0指定站点模型，这个参数指定很麻烦，如下可以参考。

您可以在一个批次中运行多个Nsites模型，如下所示，在这种情况下，将使用ncatG的默认值，如Yang等人（2000b）所述。
NSsites = 0 1 2 3 7 8
分支站点模型是通过使用变量模型和NSsites来指定的。
Model A: model = 2, NSsites = 2, fix_omega = 0

---

baseml/codeml生成 mlb，其中结尾的后面十几行以下内容。

Substitution rate is per time unit
    0.386508 +- 0.007505

根据A Step-by-Step Tutorial: Divergence Time Estimation with Approximate Likelihood Calculation Using MCMCTREE in PAML (9 June 2011)，获得m&s的估值。
由此可以获得 rgene_gamma =1 1/0.38