實(shí)驗(yàn)?zāi)康模罕容^基因組學(xué)，從基因組獲取單拷貝基因，并構(gòu)建有根的進(jìn)化樹，并估計(jì)分化時(shí)間，最終畫出有分化時(shí)間的進(jìn)化樹
實(shí)驗(yàn)流程：
orthofinder2獲取單拷貝基因——muscle多序列比對——Gblocks尋找保守位點(diǎn)——iqtree2構(gòu)建進(jìn)化樹——MCMCTREE計(jì)算分歧時(shí)間——
（1）單拷貝基因的查找 Orthofinder；
（2）多序列比對 Muscle / MAFFT / ClustalW / T-coffee， 本實(shí)驗(yàn)使用Muscle,準(zhǔn)確度最高。 mafft的準(zhǔn)確性真的低
（3）調(diào)取保守區(qū)域，并收尾連接，形成supergene ，使用軟件Gblocks　　
（4）進(jìn)化樹構(gòu)建 RaxML，PhyML或Mrbayes,本實(shí)驗(yàn)使用iqtree2
構(gòu)建樹的教程：https://www.yuque.com/wusheng/gw7a9p/mcc73y
（5）分化時(shí)間分析 divergence time，本實(shí)驗(yàn)使用mcmctree 是PAML 軟件包的一個功能。分化時(shí)間進(jìn)化樹可視化使用MCMCTreeR
（6）基因擴(kuò)張收縮分析 CAFE
（7）基因正選擇codeML PAML中一個程序

1. 提取所有基因組每個基因的最長的轉(zhuǎn)錄本
   使用getLongProteins.py腳本獲取，最后一個參數(shù)是指定相同的基因后的分隔符
   getLongProteins.py input.pep.fa out.pep.fa .
checkFa.py *.pep.fa檢測fa文件序列中是否有其他字符，如果有的話會輸出其他的字符，需要提前刪除這些其他字符。
   2.使用Orthofinder進(jìn)行分析
   把上面處理后的蛋白文件放入OrthoFinder/primary文件夾中，或者是其他名稱也行，注意修改下面的-f后文件夾名稱。一般盡量把蛋白序列命名為abbr.pep.fa這種格式。

第一次運(yùn)行使用

orthofinder -t 16 -a 16 -f OrthoFinder/primary
-t 16 -a 16是指定使用的cpu進(jìn)程數(shù)量為16
-f是指定輸入的蛋白文件的位置

如果你分析完成了，想再添加新的物種，只需要新建一個文件夾run3,把新加的物種蛋白序列放入run3。再運(yùn)行下面的命令即可

orthofinder -t 16 -a 16 -b OrthoFinder/Results_Jul08 -f run3
-b參數(shù)是指定之前分析的輸出的文件夾
-f是指定新增的蛋白序列的文件夾

如果是想刪除一個物種，同時(shí)新增一些物種，命令還是下面這樣

orthofinder -t 16 -a 16 -b OrthoFinder/Results_Jul08 -f run3
但是注意，需要進(jìn)入 OrthoFinder/Results_Jul08/WorkingDirectory這個文件夾中，找到SpeciesIDs.txt這個文件，在需要刪除的物種前加上#符號即可。新增的物種蛋白序列還是放入run3

orthofinder輸出的結(jié)果目錄文件結(jié)構(gòu)說明

Species_Tree 目錄 里面是物種的進(jìn)化樹
SpeciesTree_rooted.txt 基于所有蛋白序列構(gòu)建的進(jìn)化樹（可能不是真實(shí)的進(jìn)化關(guān)系，因?yàn)樯婕暗饺旧w倒位等結(jié)構(gòu)變異）
SpeciesTree_rooted_node_labels.txt 和上面一樣的樹，有l(wèi)abel
Single_Copy_Orthologue_Sequences 里面是單拷貝直系同源基因的fa文件，每個fa文件里是每一個物種的一條序列。
Orthogroups目錄
Orthogroups_SingleCopyOrthologues.txt 單拷貝基因的OG開頭的編號的列表文件。總行數(shù)就是每個物種里單拷貝基因的數(shù)量。
Comparative_Genomics_Statistics目錄
OrthologuesStats_one-to-one.tsv 這個文件里可以看到每個基因組和其他基因組比較，1對1的直系同源基因的數(shù)量。

如果singlecopy的基因數(shù)量比較少，這時(shí)候需要使用參數(shù)過濾，例如：在70%的物種里是單拷貝的基因的列表，這樣能夠多一些。使用getSingleCopy.py Results_XXXXX/Orthogroups/Orthogroups.GeneCount.tsv all 0.7即可獲取對應(yīng)的基因家族的名稱。
從Results_XXXXX/Orthogroup_Sequence里找到對應(yīng)家族的序列。

3.使用singlecopy2tree.bash腳本完成比對和保守位點(diǎn)提取
bash singlecopy2tree.bash OrthoFinder/Results_Jul08/Single_Copy_Orthologue_Sequences orthofinder_inputpath

4. 使用iqtree2構(gòu)建單拷貝基因的進(jìn)化樹

##使用iqtree2構(gòu)建進(jìn)化樹（此處需要16個cpu,其實(shí)少一點(diǎn)也可以，速度很快）
iqtree2 -s all.fa -T 16 -B 1000 -m MFP --bnni -o ${outgroup}

5.進(jìn)化樹進(jìn)化時(shí)間估計(jì)
從http://www.timetree.org/查詢兩個物種之間的分化時(shí)間。輸入拉丁學(xué)名，即可查詢到。

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：每一行對應(yīng)一個基因家族在該基因組對應(yīng)的基因數(shù)目；

獲取所有的直系同源基因的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是所有的單拷貝基因的列表，第一行是每個物種的名稱，每一列是對應(yīng)的基因的ID。每一行之間是單拷貝直系同源基因。