实验目的：比较基因组学，从基因组获取单拷贝基因，并构建有根的进化树，并估计分化时间，最终画出有分化时间的进化树
实验流程：
orthofinder2获取单拷贝基因——muscle多序列比对——Gblocks寻找保守位点——iqtree2构建进化树——MCMCTREE计算分歧时间——
（1）单拷贝基因的查找 Orthofinder；
（2）多序列比对 Muscle / MAFFT / ClustalW / T-coffee， 本实验使用Muscle,准确度最高。 mafft的准确性真的低
（3）调取保守区域，并收尾连接，形成supergene ，使用软件Gblocks　　
（4）进化树构建 RaxML，PhyML或Mrbayes,本实验使用iqtree2
构建树的教程：https://www.yuque.com/wusheng/gw7a9p/mcc73y
（5）分化时间分析 divergence time，本实验使用mcmctree 是PAML 软件包的一个功能。分化时间进化树可视化使用MCMCTreeR
（6）基因扩张收缩分析 CAFE
（7）基因正选择codeML PAML中一个程序

1. 提取所有基因组每个基因的最长的转录本
   使用getLongProteins.py脚本获取，最后一个参数是指定相同的基因后的分隔符
   getLongProteins.py input.pep.fa out.pep.fa .
checkFa.py *.pep.fa检测fa文件序列中是否有其他字符，如果有的话会输出其他的字符，需要提前删除这些其他字符。
   2.使用Orthofinder进行分析
   把上面处理后的蛋白文件放入OrthoFinder/primary文件夹中，或者是其他名称也行，注意修改下面的-f后文件夹名称。一般尽量把蛋白序列命名为abbr.pep.fa这种格式。

第一次运行使用

orthofinder -t 16 -a 16 -f OrthoFinder/primary
-t 16 -a 16是指定使用的cpu进程数量为16
-f是指定输入的蛋白文件的位置

如果你分析完成了，想再添加新的物种，只需要新建一个文件夹run3,把新加的物种蛋白序列放入run3。再运行下面的命令即可

orthofinder -t 16 -a 16 -b OrthoFinder/Results_Jul08 -f run3
-b参数是指定之前分析的输出的文件夹
-f是指定新增的蛋白序列的文件夹

如果是想删除一个物种，同时新增一些物种，命令还是下面这样

orthofinder -t 16 -a 16 -b OrthoFinder/Results_Jul08 -f run3
但是注意，需要进入 OrthoFinder/Results_Jul08/WorkingDirectory这个文件夹中，找到SpeciesIDs.txt这个文件，在需要删除的物种前加上#符号即可。新增的物种蛋白序列还是放入run3

orthofinder输出的结果目录文件结构说明

Species_Tree 目录 里面是物种的进化树
SpeciesTree_rooted.txt 基于所有蛋白序列构建的进化树（可能不是真实的进化关系，因为涉及到染色体倒位等结构变异）
SpeciesTree_rooted_node_labels.txt 和上面一样的树，有label
Single_Copy_Orthologue_Sequences 里面是单拷贝直系同源基因的fa文件，每个fa文件里是每一个物种的一条序列。
Orthogroups目录
Orthogroups_SingleCopyOrthologues.txt 单拷贝基因的OG开头的编号的列表文件。总行数就是每个物种里单拷贝基因的数量。
Comparative_Genomics_Statistics目录
OrthologuesStats_one-to-one.tsv 这个文件里可以看到每个基因组和其他基因组比较，1对1的直系同源基因的数量。

如果singlecopy的基因数量比较少，这时候需要使用参数过滤，例如：在70%的物种里是单拷贝的基因的列表，这样能够多一些。使用getSingleCopy.py Results_XXXXX/Orthogroups/Orthogroups.GeneCount.tsv all 0.7即可获取对应的基因家族的名称。
从Results_XXXXX/Orthogroup_Sequence里找到对应家族的序列。

3.使用singlecopy2tree.bash脚本完成比对和保守位点提取
bash singlecopy2tree.bash OrthoFinder/Results_Jul08/Single_Copy_Orthologue_Sequences orthofinder_inputpath

4. 使用iqtree2构建单拷贝基因的进化树

##使用iqtree2构建进化树（此处需要16个cpu,其实少一点也可以，速度很快）
iqtree2 -s all.fa -T 16 -B 1000 -m MFP --bnni -o ${outgroup}

5.进化树进化时间估计
从http://www.timetree.org/查询两个物种之间的分化时间。输入拉丁学名，即可查询到。

image.png

#R8S安装方法
wget -q https://sourceforge.net/projects/r8s/files/r8s1.81.tar.gz \
  && tar -zxvf r8s1.81.tar.gz \
  && cd r8s1.81/src \
  && sed -i 's|/usr/include/sys/errno.h||' Makefile.linux \
  && sed -i 's/continuousML.o //' Makefile.linux \
  && sed -i 's/continuousML.o:/#continuousML.o:/' Makefile.linux \
  && make -f Makefile.linux

OrthoFinder/Results_Jul08/Orthogroups/Orthogroups.GeneCount.tsv：每一行对应一个基因家族在该基因组对应的基因数目；

获取所有的直系同源基因的ID列表

cat name.list|tr ">" "\t"  >SingleGeneID.tsv
awk '
{
    for (i=1; i<=NF; i++) {
        if (NR == 1) {
            transposed[i] = $i
        } else {
            transposed[i] = transposed[i] "\t" $i
        }
    }
}
END {
    for (i=1; i<=NF; i++) {
        print transposed[i]
    }
}' SingleGeneID.tsv > Specs.singlecopy.id
tail -1 Specs.singlecopy.id >header
sed -i '$ d' Specs.singlecopy.id
cat header Specs.singlecopy.id > Species.singlecopy.tsv
rm -rf SingleGeneID.tsv Specs.singlecopy.id header 

最终Species.singlecopy.tsv是所有的单拷贝基因的列表，第一行是每个物种的名称，每一列是对应的基因的ID。每一行之间是单拷贝直系同源基因。