0.软件的安装略
1.物种树的构建
参见 https://www.jianshu.com/p/336b65ca1b67
假设我们已经通过上述网页拿到了A-c.phy和相应的物种树文件

2.估计碱基替换速率(对应软件版本4.8, 4.9的参数和本文不同)

软件版本4.8

/data/01/user157/software/paml4.8/bin/baseml baseml.ctl

baseml.input.tree的内容:

6 1
((Bamboo,XD),((sppCa,sppJu)'@.001',(sppGa,sppGo)));
#需要标定时间，@.001表示大约在10万年
#标定的时间可以用timetree获得，输入两个物种，可以得到这两个物种的最晚分化时间

baseml.ctl内容如下, 前两行是上面提到的输入文件:

      seqfile = A-c.phy
     treefile = baseml.input.tree

      outfile = mlb       * main result file
        noisy = 3   * 0,1,2,3: how much rubbish on the screen
      verbose = 0   * 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
                    * 5:T92, 6:TN93, 7:REV, 8:UNREST, 9:REVu; 10:UNRESTu

        Mgene = 0   * 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff

        clock = 1   * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
    fix_kappa = 0   * 0: estimate kappa; 1: fix kappa at value below; 2: kappa for branches
        kappa = 2  * initial or fixed kappa

    fix_alpha = 0    * 0: estimate alpha; 1: fix alpha at value below
        alpha = 0.5   * initial or fixed 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
        nparK = 0   * rate-class models. 1:rK, 2:rK&fK, 3:rK&MK(1/K), 4:rK&MK 

        nhomo = 0   * 0 & 1: homogeneous, 2: kappa for branches, 3: N1, 4: N2
        getSE = 1   * 0: don't want them, 1: want S.E.s of estimates
 RateAncestor = 0   * (0,1,2): rates (alpha>0) or ancestral states

   Small_Diff = 7e-6
    cleandata = 0  * remove sites with ambiguity data (1:yes, 0:no)?
*        icode = 0  * (with RateAncestor=1. try "GC" in data,model=4,Mgene=4)
*  fix_blength = 1  * 0: ignore, -1: random, 1: initial, 2: fixed, 3: proportional
       method = 0  * Optimization method 0: simultaneous; 1: one branch a time

输出文件mlb中会有这么几行，记住这个数字, 这个就是替换率

Substitution rate is per time unit
   0.395351 +- 0.000685

3.第一次运行mcmctree
/data/01/user157/software/paml4.8/bin/mcmctree mcmctree.ctl
输出文件out.BV, 用作下面步骤的输入文件

mcmctreel.input.tree的内容:

6 1
((Bamboo,XD),((sppCa,sppJu)'>.001',(sppGa,sppGo)));
#需要标定时间，>.001表示大于10万年

          seed = -1
       seqfile = A-c.phy
      treefile = mcmc.input.tree
       outfile = out

         ndata = 1
       seqtype = 0  * 0: nucleotides; 1:codons; 2:AAs
       usedata = 3    * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
         clock = 2    * 1: global clock; 2: independent rates; 3: correlated rates
       RootAge = <1.0  * safe constraint on root age, used if no fossil for root.

         model = 0    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
         alpha = 0    * alpha for gamma rates at sites
         ncatG = 5    * No. categories in discrete gamma

     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?

       BDparas = 1 1 0    * birth, death, sampling
   kappa_gamma = 6 2      * gamma prior for kappa
   alpha_gamma = 1 1      * gamma prior for alpha

   rgene_gamma = 1 2.529397927411 1  * gamma prior for overall rates for genes ### 1/替换率
  sigma2_gamma = 1 10 1   * gamma prior for sigma^2     (for clock=2 or 3)

      finetune = 1: .05  0.1  0.12  0.1 .3  * auto (0 or 1) : times, rates, mixing, paras, RateParas, FossilErr

         print = 1
        burnin = 500000
      sampfreq = 5000
       nsample = 20000

cat out.BV > in.BV #用作后面的输入文件

4.第二次运行mcmctree,把上面的in.BV拷贝在第二次运行的目录下
/data/01/user157/software/paml4.8/bin/mcmctree mcmctree.ctl

这一次的mcmctree.ctl的内容：(区别在于usedata = 2

         seed = -1
       seqfile = A-c.phy
      treefile = mcmc.input.tree
       outfile = out

         ndata = 1
       seqtype = 0  * 0: nucleotides; 1:codons; 2:AAs
       usedata = 2    * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
         clock = 2    * 1: global clock; 2: independent rates; 3: correlated rates
       RootAge = <1.0  * safe constraint on root age, used if no fossil for root.

         model = 0    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
         alpha = 0    * alpha for gamma rates at sites
         ncatG = 5    * No. categories in discrete gamma

     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?

       BDparas = 1 1 0    * birth, death, sampling
   kappa_gamma = 6 2      * gamma prior for kappa
   alpha_gamma = 1 1      * gamma prior for alpha

   rgene_gamma = 1 2.292111240743 1  * gamma prior for overall rates for genes  ### 1/替换率
  sigma2_gamma = 1 10 1   * gamma prior for sigma^2     (for clock=2 or 3)

      finetune = 1: .05  0.1  0.12  0.1 .3  * auto (0 or 1) : times, rates, mixing, paras, RateParas, FossilErr

         print = 1
        burnin = 500000
      sampfreq = 5000
       nsample = 20000

*** Note: Make your window wider (100 columns) before running the program.

然后就会输出FigTree.tre，也就是最后的结果了

5.可以重复第四步，结果稳定即可
mcmc.ctl的最后几个参数可以尝试更换

6.我们就拿到了Cafe的输入文件（raw）
内容如下：

#NEXUS
BEGIN TREES;

        UTREE 1 = ((Bamboo: 0.507727, XD: 0.507727) [&95%={0.266, 0.907}]: 0.055557, ((sppCa: 0.010037, sppJu: 0.010037) [&95%={0.004, 0.020}]: 0.010716, (sppGa: 0.012062, sppGo: 0.012062) [&95%={0.005, 0.024}]: 0.008690) [&95%={0.009, 0.041}]: 0.542531) [&95%={0.286, 0.978}];

END;

需要删掉空格和置信区间，空格一定要删干净，这个就是cafe要求的输入文件格式
我命名为FigTree.nwk

#grep  "UTREE 1 =" FigTree.tre | sed -E -e "s/\[[^]]*\]//g" -e "s/[ \t]//g" -e "/^$/d" -e "s/UTREE1=//" > FigTree.nwk

((Bamboo:0.507727,XD:0.507727):0.055557,((sppCa:0.010037,sppJu: 0.010037):0.010716,(sppGa:0.012062,sppGo:0.012062):0.008690):0.542531);

7.cafe的安装可以用conda，这里略，然后准备第二个输入文件，也是OrthoFinder的结果

awk -v  OFS="\t" '{if($1=="Orthogroup"){print"Descript",$1,$2,$3,$4,$5,$6,$7}else{print"(null)",$1,$2,$3,$4,$5,$6,$7}}' Orthogroups.GeneCount.tsv > GeneCounts.tsv

python /data/01/user158/kuangzhuoran/software/CAFE5/docs/tutorial/clade_and_size_filter.py -i GeneCounts.tsv -o gene_family_filter.txt -s

cafe5 -i gene_family_filter.txt -t FigTree.nwk -o out -c 20

8.输出的文件都仔细看一下, 提取某个物种显著扩张收缩的基因

cat Base_family_results.txt | grep "y"|cut -f1 >p0.05.significant
#提取显著的OG

head -n 1 Base_change.tab > tmp1
grep -f p0.05.significant Base_change.tab > tmp2
cat tmp1 tmp2 > Base_p0.05change.tab
#提取对应OG的expand和expansion数

cat Base_p0.05change.tab | cut -f1,4 | grep "+[1-9]" | cut -f1  > sppJu.significant.expand
#cut -f1,4 的这个4，就是对应的节点/物种，我想要研究的物种
#grep "+[1-9]" 就是提取显著扩张的

cat Base_p0.05change.tab | cut -f1,4 | awk '{if($2!="+0") print}' | grep "-" | cut -f1 > sppJu.significant.expansion
#cut -f1,4 的这个4，就是对应的节点/物种，我想要研究的物种
#grep "-" 就是提取显著收缩的

#这个Orthogroups.tsv也是OrthoFinder的结果
grep -f judaei.significant.expansion Orthogroups.tsv | cut -f7 | sed "s/ /\n/g" | sed "s/\t/\n/g" | sed "s/,//g" | sort | uniq > sppJu.significant.expansion.gene
#cut -f7 就是就是对应的节点/物种
#这样就拿到了显著收缩的基因

9.如果是提取某个节点显著收缩与扩张的基因

awk 'NR!=1 && $13>0 {print $0}' Base_count.tab | cut -f1 > node11.orthogroups
#比如我要提取Node11节点

grep -f node11.orthogroups Orthogroups.txt |sed "s/ /\n/g" | grep -E "carmeli|galili|golani|judaei" | sort | uniq > node11.genes
#这个就是背景库

cat Base_p0.05change.tab | cut -f1,13 | grep "+[1-9]" | cut -f1  > node11.significant.expand
grep -f node11.significant.expand Orthogroups.tsv | cut -f7,8,9,10 | sed "s/ /\n/g" | sed "s/\t/\n/g" | sed "s/,//g" | sort | uniq > node11.significant.expand.gene

cat Base_p0.05change.tab | cut -f1,13 | awk '{if($2!="+0") print}' | grep "-" | cut -f1 > node11.significant.expansion
grep -f node11.significant.expansion Orthogroups.tsv | cut -f7,8,9,10 | sed "s/ /\n/g" | sed "s/\t/\n/g" | sed "s/,//g" | sort | uniq > node11.significant.expansion.gene
#Node11显著扩张的基因
#cut -f1,13 的这个13，就是Node11节点对应的列

#node11.genes、node11.significant.expand.gene、node11.significant.expansion.gene就是后续做富集要用到的