一、數(shù)據(jù)準(zhǔn)備

首先拐迁，選取不同物種的Protein數(shù)據(jù)集：Arabidopsis_thaliana.fa；Citrus_grandis.fa；Dimocarpus_longan.fa措伐；Durio_zibethinus.fa；Prunus_persica.fa军俊； Vitis_vinifera.fa侥加；Citrus_clementina.fa；Citrus_sinensis.fa粪躬；Diospyros_oleifera.fa担败；Malus_domestica.fa昔穴；Oryza_sativa.fa；Pyrus_communis.fa
然后進(jìn)行數(shù)據(jù)處理提前，去冗余吗货，只保留最長(zhǎng)轉(zhuǎn)錄本，去除可變剪切：
python3 removeRedundantProteins.py -i input.fa -o output.fa
removeRedundantProteins.py

import sys,getopt
def usage():
    print('usage:python3 removeRedundantProteins.py -i <in_fasta> -o <out_fasta> <-h>')
    return
def removeRedundant(in_file,out_file):
    gene_dic = {}
    flag = ''
    with open (in_file) as in_fasta:
        for line in in_fasta:
            if '>' in line:
                line1 = line.strip('>\n')
                line2 = line1.split('.')
                li = line2[0]
                flag = li
                try:
                    gene_dic[li]
                except KeyError:
                    gene_dic[li] = [line]
                else:
                    gene_dic[li].append(line)
            else:
                gene_dic[flag][-1] += line
    with open (out_file,'w') as out_fasta:
        for k,v in gene_dic.items():
            if len(v) == 1:
                out_fasta.write(gene_dic[k][0])
            else:
                trans_max = ''
                for trans in gene_dic[k]:
                    a = len(list(trans))
                    b = len(list(trans_max))
                    if a > b:
                        trans_max = trans
                out_fasta.write(trans_max)
def main(argv):
    try:
        opts, args = getopt.getopt(argv,'hi:o:')
    except getopt.GetoptError:
        usage()
        sys.exit()
    for opt, arg in opts:
        if opt == '-h':
            usage()
            sys.exit()
        elif opt == '-i':
            in_fasta_name = arg
        elif opt == '-o':
            outfile_name = arg
    try:
        removeRedundant(in_fasta_name,outfile_name)
    except UnboundLocalError:
        usage()
    return
if __name__ == '__main__':
    main(sys.argv[1:])

將處理好的數(shù)據(jù)置于一個(gè)文件夾中“Dataset”

二岖研、使用OrthoFinder查找直系同源基因

OrthoFinder這個(gè)軟件卿操，之前有一篇文章已經(jīng)介紹過(guò)了，這里就不在贅述孙援，這個(gè)軟件安裝十分友好害淤，直接conda安裝即可；
nohup orthofinder -f Dataset -M msa -S diamond -T iqtree -t 24 -a 24 2> orthofinder.log &
orthofinder參數(shù)詳情：
-t 并行序列搜索線程數(shù)（默認(rèn)= 16）
-a 并行分析線程數(shù)（默認(rèn)值= 1）
-M 基因樹推斷方法拓售】悖可選：dendroblast和msa（默認(rèn)= dendroblast）
-S 序列搜索程序（默認(rèn)= blast）選項(xiàng)：blast，mmseqs,础淤，blast_gz崭放，diamond（推薦使用diamond，比對(duì)速度很給力）
-A 多序列聯(lián)配方式鸽凶，需要添加參數(shù)-M msa時(shí)才有效币砂；（默認(rèn)= mafft）可選擇：muscle，mafft
-T 建樹方法玻侥，需要添加參數(shù)-M msa時(shí)才有效决摧，（默認(rèn) = fasttree）可選：iqtree，raxml-ng凑兰，fasttree掌桩，raxml
-s <文件> 可指定特定的根物種樹
-I 設(shè)定MCL的通脹參數(shù)（默認(rèn) = 1.5）
-x Info用于以othoXML格式輸出結(jié)果
-p <dir>將臨時(shí)pickle文件寫入到<dir>
-l 只執(zhí)行單向序列搜索
-n 名稱以附加到結(jié)果目錄
-h 打印幫助文本
如果只需要查找直系同源基因，只需接“-f” 參數(shù)即可姑食；此步也可建樹波岛，采用默認(rèn)的建樹方法fasttree，為無(wú)根樹音半。
nohup orthofinder -f Dataset &
如果添加-M msa -T iqtree設(shè)定制定參數(shù)则拷，可按照設(shè)定的參數(shù)使用最大似然法構(gòu)建有根的物種進(jìn)化樹，構(gòu)建的樹為STAG樹曹鸠。
nohup orthofinder -f Dataset -M msa -S diamond -T iqtree -t 24 -a 24 2> orthofinder.log &

關(guān)于構(gòu)建系統(tǒng)進(jìn)化樹煌茬，有很多種做法，常見的有利用物種全部的蛋白序列物延，構(gòu)建STAG物種樹；也有使用單拷貝直系同源基因構(gòu)建的物種進(jìn)化樹仅父，關(guān)于這一點(diǎn)叛薯，OrthoFinder查找同源基因浑吟，可以輸出直系單拷貝同源基因的序列結(jié)果，后續(xù)也可使用其他構(gòu)樹軟件及算法進(jìn)行進(jìn)化樹構(gòu)建耗溜。關(guān)于建樹方法组力，則有距離矩陣法、最大簡(jiǎn)約法抖拴、最大似然法以及貝葉斯燎字；當(dāng)然目前主流采用的基本為最大似然法和貝葉斯，其中貝葉斯算法計(jì)算量巨大阿宅，耗時(shí)最久候衍，其構(gòu)建的樹也認(rèn)為最為“逼真”，但文章中使用較多的還是最大似然法洒放，其耗時(shí)也需蠻久蛉鹿。

三、OrthoFinder輸出結(jié)果解讀

OrthoFinder輸出的結(jié)果會(huì)在OrthoFinder文件夾下面的以日期命名的文件夾中往湿，如：~/OrthoFinder/Results_May08

orthofinder輸出結(jié)果

其中欧漱，我們可以用Orthogroups.GeneCount.tsv來(lái)作為CAFE的輸入文件，分析基因家族的擴(kuò)張與收縮葬燎；使用SpeciesTree_rooted.txt作為推斷的物種樹误甚，并使用r8s，從中提取超度量樹（ultrametric tree）即時(shí)間樹谱净；

四窑邦、用推斷的物種樹提取超度量樹（ultrametric tree）

1、軟件準(zhǔn)備（r8s）

###Install:
wget http://ceiba.biosci.arizona.edu/r8s/r8s1.81.tar.gz 
tar -xzvf  r8s1.81.tar.gz 
cd r8s1.81/src
make -f Makefile.linux

2壕探、準(zhǔn)備r8s輸入文件r8s_ctl_file.txt

python cafetutorial_prep_r8s.py -i SpeciesTree_rooted.txt -o r8s_ctl_file.txt -s 6650255 -p 'Oryza_sativa,Arabidopsis_thaliana' -c '152'
參數(shù)：
-i path_tree_file: path to .txt file containing tree in NEWICK format
-s n_sites: number of sites in alignment that was used to infer species tree
-p list_of_spp_tuples: list of tuples (each tuple being two species IDs whose mrca's age we are constraining; e.g., [('ENSG00','ENSPTR'),('ENSFCA','ENSECA')]
-c list_of_spp_cal_points: list of flats, one for each tuple in list_of_spp_tuples (e.g., [6.4,80])
-s 即用于推斷物種樹的比對(duì)序列堿基數(shù)目冈钦；
-p 已知物種樹中的一對(duì)物種；
-c 已知一對(duì)物種的分化年限：
可在timetree網(wǎng)站查詢：為152 mya

Pairwise Divergence Time For Oryza sativa & Arabidopsis thaliana

__author__ = "Fabio H. K. Mendes"
import os
import argparse
def r8s_prep(path_tree_file, sites_n, list_of_spp_tuples, list_of_spp_cal_points, path_output_file):
    tree_str = str()
    with open(path_tree_file, "r") as tree_file:
        tree_str = tree_file.readline().rstrip()
    fixage_lines = list()
    with open(path_output_file, "w") as output_file:
        output_file.write(
            "#NEXUS\nbegin trees;\ntree nj_tree = [&R] " + tree_str + "\nEnd;\nbegin rates;\nblformat nsites=" + \
            sites_n + " lengths=persite ultrametric=no;\ncollapse;\n")
        for idx, pair in enumerate(list_of_spp_tuples):
            node_name = pair[0][-3:] + pair[1][-3:]
            output_file.write(
                "mrca " + node_name + " " + pair[0] + " " + pair[1] + ";\n")
            fixage_lines.append("fixage taxon=" + node_name + " age=" + list_of_cal_points[idx] + ";\n")
        for line in fixage_lines:
            output_file.write(line)
        output_file.write(
            "divtime method=pl algorithm=tn cvStart=0 cvInc=0.5 cvNum=8 crossv=yes;\ndescribe plot=chronogram;\ndescribe plot=tree_description;\nend;"
            )
if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=__doc__, prog="mcl2rawcafe.py")
    parser.add_argument("-i", "--input-file", action="store", dest="input_file", required=True, type=str, help="full path to .txt file containing tree in NEWICK format")
    parser.add_argument("-o", "--output-file", action="store", dest="output_file", required=True, type=str, help="full path to file to be written (r8s input file)")
    parser.add_argument("-s", "--sites-n", action="store", dest="sites_n", required=True, type=str, help="number of sites in alignment used to estimate species tree")
    parser.add_argument("-p", "--pairs-species", action="store", dest="spp_pairs", required=True, type=str, help="")
    parser.add_argument("-c", "--calibration-points", action="store", dest="cal_points", required=True, type=str, help="")
    args = parser.parse_args()
    list_of_spp_tuples = list()
    for pair in args.spp_pairs.split(" "):
        list_of_spp_tuples.append(tuple(pair.split(",")))
    list_of_cal_points = args.cal_points.split(",")
    print "\nRunning cafetutorial_clade_and_size_filter.py as a standalone...\n"
    r8s_prep(args.input_file, args.sites_n, list_of_spp_tuples, list_of_cal_points, args.output_file)

3李请、運(yùn)行r8s并提取超度量樹

/home/Tools/r8s1.81/src/r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt
#提取超度量樹
$  tail -n 1 r8s_tmp.txt | cut -c 16- > twelve_spp_r8s_ultrametric.txt

五瞧筛、基因家族擴(kuò)張與收縮分析

1厉熟、軟件安裝（CAFE）

conda install cafe

2、cafe輸入數(shù)據(jù)處理

awk -F"\t" '{print "Desc""\t"$1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11"\t"$12"\t"$13}' Orthogroups.GeneCount.tsv > cafe-input.data 
sed -e '1 s/(null)/Desc/' -e '1 s/Orthogroup/Family/' -i cafe-input.data
python cafetutorial_clade_and_size_filter.py -i cafe-input.data -o filtered_cafe_input.txt -s 2> filtered.log 

cafetutorial_clade_and_size_filter.py

__author__ = "Fabio H. K. Mendes"
import os
import argparse
def clade_filter(mcl_dump, clade_str):
    lines_to_keep_list = list()
    spp_idx_dict = dict()
    clades_list = list()
    if clade_str: # if clade filter was specified
        clades_list = clade_str.split(" ")
    with open(mcl_dump, "r") as input_file:
        for line_n, line in enumerate(input_file):
            line = line.rstrip()
            tokens = line.split("\t")
            spp_info = tokens[2:]
            if line.startswith("Desc"):
                spp_idx_dict = dict((sp, idx) for idx,sp in enumerate(spp_info))
                continue
            if clades_list:                
                clades_ok_list = list()
                for clade in clades_list:
                    spp_list = clade.split(",")
                    clade_count = sum(1 for sp in spp_list if int(spp_info[spp_idx_dict[sp]]) >= 1)
                    if clade_count >= 2:
                        clades_ok_list.append(1)
                if sum(clades_ok_list) == len(clades_list):
                    lines_to_keep_list.append(line_n)               
            # just keeping lines where >=2 species (among all of them) have gene copies
            clade_count = sum(1 for sp_count in spp_info if int(sp_count) >= 1)
            if clade_count >= 2:
                lines_to_keep_list.append(line_n)              
    return set(lines_to_keep_list)
def size_filter(mcl_dump, lines_to_keep_set):
    lines_to_remove_set = set()
    size_cutoff = 100
    fam_size = int()
    with open(mcl_dump, "r") as input_file:
        for line_n, line in enumerate(input_file):
            line = line.rstrip()
            if line.startswith("Desc"):
                continue
            elif line_n not in lines_to_keep_set and len(lines_to_keep_set) > 0:
                continue
            tokens = line.split("\t")
            spp_info = tokens[2:]
            for gene_count in spp_info:
                if int(gene_count) >= size_cutoff:
                    lines_to_separate_set.add(line_n)                 
    lines_to_keep_set -= lines_to_separate_set
    return lines_to_keep_set, lines_to_separate_set
def filter_print(mcl_dump, lines_to_keep_set, lines_to_separate_set, output_file_name):
    if len(lines_to_keep_set) == 0 and len(lines_to_separate_set) == 0:
        exit("No filtering was done! Exiting...\n")
    with open(output_file_name, "w") as output_file:
        with open("large_"+output_file_name, "w") as output_file2:
            with open(mcl_dump, "r") as input_file:
                for line_n, line in enumerate(input_file):
                    line = line.rstrip() + "\n"
                    if line_n == 0:
                        output_file.write(line)
                        output_file2.write(line) 
                    elif line_n in lines_to_keep_set and len(lines_to_keep_set) >= 1:
                        output_file.write(line)
                    elif line_n not in lines_to_separate_set and len(lines_to_keep_set) == 0:
                        output_file.write(line)
                    # has to be if, not elif
                    if line_n in lines_to_separate_set:
                        output_file2.write(line)
        # cleaning up in case size filtering was not done
        if len(lines_to_separate_set) == 0:
            os.unlink("large_"+output_file_name)
if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=__doc__, prog="cafetutorial_clade_and_size_filter.py")
    parser.add_argument("-i", "--input-file", action="store", dest="input_file", required=True, type=str, help="full path to mcl's output dump file")
    parser.add_argument("-o", "--output-file", action="store", dest="output_file", required=True, type=str, help="full path to file to be written")
    parser.add_argument("-cl", "--clade-filter", action="store", dest="clade_str", default=None, required=False, type=str, help="list of clades (separated by white spaces) comprised of species identifiers (separated by comma) that must have at least two species with gene copies for a given gene family")
    parser.add_argument("-s", "--size-filter", action="store_true", dest="size_filter", required=False, help="option to perform size filtering")
    args = parser.parse_args()
    lines_to_keep_set, lines_to_separate_set = set(), set()
    # applying size filter (if no groups are specified, just the lines where just 1 species has gene counts are removed
    lines_to_keep_set = clade_filter(args.input_file, args.clade_str)
    if args.size_filter:
        lines_to_keep_set, lines_to_separate = size_filter(args.input_file, lines_to_keep_set)
    filter_print(args.input_file, lines_to_keep_set, lines_to_separate_set, args.output_file) # .add(0) to get header back

3较幌、配置cafe文件

vim cafetutorial_run.sh

cafetutorial_run.sh

tree即為r8s提取的超度量樹揍瑟；

4、運(yùn)行cafe

mkdir -p reports
cafe cafetutorial_run.sh

5乍炉、結(jié)果統(tǒng)計(jì)

python cafetutorial_report_analysis.py -i reports/report_run.cafe -o reports/summary_run
summary_run_node.txt：統(tǒng)計(jì)每個(gè)節(jié)點(diǎn)中擴(kuò)張绢片，收縮的基因家族數(shù)目；
summary_run_fams.txt：具體發(fā)生變化的基因家族

六岛琼、CAFE_fig可視化

python3 /home/Tools/CAFE_fig/CAFE_fig.py resultfile.cafe -pb 0.05 -pf 0.05 --dump test/ -g svg --count_all_expansions
輸出svg格式的文件底循，可導(dǎo)入AI編輯美化；

CAFE_fig運(yùn)行報(bào)錯(cuò)：（module 'ete3' has no attribute 'TreeStyle'）
報(bào)錯(cuò)解決：
vim /home/Tools/CAFE_fig/CAFE_fig.py

CAFE_fig.py

程序還在運(yùn)行此叠，后續(xù)貼出結(jié)果圖。

參考

OrthoFinder
timetree
http://www.chenlianfu.com/?tag=genomecomparison
http://www.reibang.com/p/146093c91e2b
r8s
【OrthoFinder】
Emms, D.M., Kelly, S. OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy. Genome Biol 16, 157 (2015) (https://doi.org/10.1186/s13059-015-0721-2)
Emms, D.M., Kelly, S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol 20, 238 (2019) https://doi.org/10.1186/s13059-019-1832-y)
【CAFE v.4.2.1】
Han, M. V., Thomas, G. W. C., Lugo-Martinez, J., and Hahn, M. W. Estimating gene gain and loss rates in the presence of error in genome assembly and annotation using CAFE 3. Molecular Biology and Evolution 30, 8 (2013)
【iqtree v. 1.6.12】
Lam-Tung Nguyen, Heiko A. Schmidt, Arndt von Haeseler, and Bui Quang Minh (2015) IQ-TREE: A fast and effective stochastic algorithm for estimating maximum likelihood phylogenies. Mol Biol Evol, 32:268-274. https://doi.org/10.1093/molbev/msu300
【modelFinder】
Subha Kalyaanamoorthy, Bui Quang Minh, Thomas KF Wong, Arndt von Haeseler, and Lars S Jermiin (2017) ModelFinder: Fast model selection for accurate phylogenetic estimates. Nature Methods, 14:587–589. https://doi.org/10.1038/nmeth.4285
【R8s v. 1.81】
Sanderson M J. R8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock. Bioinformatics, 2003, 19(2): 301-302.
【STAG tree】
Emms D.M. & Kelly S. STAG: Species Tree Inference from All Genes (2018), bioRxiv https://doi.org/10.1101/267914

補(bǔ)充材料

直系同源低拷貝核基因（orthologous low-copy nuclear genes, LCN）：在進(jìn)化過(guò)程中随珠，新基因通常來(lái)自事先存在的基因灭袁，新基因的功能從先前基因的功能進(jìn)化而來(lái)。新基因的原材料來(lái)自基因組區(qū)域的重復(fù)窗看，這種重復(fù)可包括一個(gè)或多個(gè)基因茸歧。作為物種形成的伴隨事件而被重復(fù)，并繼續(xù)保持相同功能的基因显沈，稱為直系同源基因（orthologous）软瞎。新的基因功能可由在單個(gè)物種的基因組中發(fā)生的重復(fù)引起的。在一個(gè)基因組內(nèi)部的重復(fù)導(dǎo)致旁系同源基因（paralogous gene）拉讯。
最大似然法（maximum likelihood method）：使用概率模型涤浇，尋找能夠以較高概率產(chǎn)生觀察數(shù)據(jù)的系統(tǒng)發(fā)生樹。
外群的選擇：大多數(shù)的種系發(fā)生重建方法會(huì)產(chǎn)生無(wú)根樹魔慷，但是觀察樹的拓?fù)浣Y(jié)構(gòu)無(wú)法識(shí)別樹根應(yīng)在哪一分支上只锭。實(shí)際上，對(duì)于要證實(shí)哪一個(gè)分類單元的分支先于其他的分類單元院尔，樹根必須確定蜻展。在無(wú)根樹中設(shè)定一個(gè)根，最簡(jiǎn)單的方法是在數(shù)據(jù)集中增加一個(gè)外群（outgroup）邀摆。外群是一種分類操作單元纵顾，且有外部信息表明外群在所有分類群之前就已分化。研究演化歷史栋盹，一般選擇比目標(biāo)序列具有較早進(jìn)化歷史的序列作為外類群施逾。
Bootstrap support: bootstrap是統(tǒng)計(jì)學(xué)上一種非參數(shù)統(tǒng)計(jì)方法，通過(guò)有放回的隨機(jī)抽樣，構(gòu)建分類回歸樹汉额。Jackknife與bootstrap類似沪饺，只是每次抽樣時(shí)會(huì)去除幾個(gè)樣本，像小刀一樣切去一部分闷愤。所謂bootstrap法就是從整個(gè)序列的堿基（氨基酸）中任意選取一半，剩下的一半序列隨機(jī)補(bǔ)齊組成一個(gè)新的序列件余。這樣讥脐，一個(gè)序列就可以變成許多序列，一個(gè)序列組也就可以變成許多個(gè)序列組啼器。根據(jù)某種算法（距離矩陣法旬渠、最大簡(jiǎn)約法、最大似然法）端壳，每個(gè)多序列組都可以生成一個(gè)進(jìn)化樹告丢。將生成的許多進(jìn)化樹進(jìn)行比較，按照多數(shù)規(guī)則（majority-rule）就會(huì)得到一個(gè)最“逼真”的進(jìn)化樹损谦。