【比较基因组学】01、基因家族聚类分析

        为进行多物种的系统发育分析，进行基因家族的聚类工作是必不可少的。

物种数据准备        

        首先，需要根据拟研究的对象去准备物种数据。如果做植物的话建议尽量选择在 phytozome 上下载，注册账号简单快捷，最关键的是该数据库格式非常规范，对后续分析有很大帮助。

        根据选择物种信息下载数据，并基于数据信息提取最长 cds 以用于后续基因家族聚类分析的输入。

        这里需要下载物种的 .fa .cds .pep .gff 文件，数据包括：

        # 基因注释文件 xxx.gff3

        # cds序列文件 xxx.cds.all.fa

        # 蛋白序列文件 xxx.pep.all.fa

        # 基因组序列文件 xxx.fa

        基于 gff3 注释文件提取最长 cds 序列 ID ，并过滤 gff3 文件（以 Acorus_americanus 为例，简写为属名第一个字母加种名前三位字母）：

perl ~/software/gff_ensembl_longest.pl   Acorus_americanus_v1.1.gene_exons.gff3 \

Aame_longest.gene2mrna_id  Aame_longest.gff3

        得到最长基因及 mRNA ID 后，通过 cds 文件和 pep 文件提取该物种最长 cds 文件和 pep 文件：

        观察到在原 cds 文件和 pep 文件中每条序列的命名方式与提取的 mRNA ID 不同，需要分别做出修改。当然，如果实际情况命名格式一样的话那最好。

        首先获取用于蛋白序列提取的 ID 文件：

less -S Aame_longest.gene2mrna_id | awk '{print $2}' | sed 's/.v1.1$/.p/g' > Aame_longest.mrna_id\(4pep\)

        再获取用于 cds 序列提取的 ID 文件：

less -S Aame_longest.mrna_id\(4pep\) | sed 's/.p$//' > Aame_longest.mrna_id\(4cds\)

        最后基于以上两个 ID 文件用 seqtk 软件去提取最长 cds 序列和 pep 序列：

seqtk subseq Acorus_americanus_v1.1.cds.fa  Aame_longest.mrna_id\(4cds\) > Aame_longest.cds.fasta

seqtk subseq Acorus_americanus_v1.1.protein.fa  Aame_longest.mrna_id\(4pep\) > Aame_longest.pep.fasta

        其他物种信息处理方式类似。

        提取完成所有输入物种的最长 cds 序列和 pep 序列后，将蛋白序列作为输入，进行基因家族聚类。

        基因家族聚类常用软件为 OrthoFinder 和 orthoMCL，orthoMCL 使用较为繁琐；而 orthofinder 集成度更好，可自动完成建树等工作，这里使用 orthofinder 进行基因家族聚类。

        将所有输入物种的蛋白序列放入 /GeneFamilyCluster/ 文件夹下。蛋白序列文件名称统一为 “物种简写.fasta”。OrthoFinder 会根据 fasta 文件前缀汇总统计表格。

        最后运行 orthofinder ：

orthofinder \

-f /home/zijie/workspace/01.GeneFamilyCluster/  \  #输入数据目录

-S diamond \  #比对方法，blast, mmseqs, blast_gz, diamond，使用 diamond 提升序列比对速度

-M msa \  #基因树推断方法，dendroblast，msa 软件默认使用 dendroblast 进行比对，使用 msa 即多序列比对的方法进行物种树的构建提高准确性

-T fasttree \  #建树软件选择，iqtree, raxml-ng, fasttree, raxml

-t 30\  #线程数，默认为56

        运行完成后，主要聚类结果位于 Orthogroups 目录下。 MultipleSequenceAlignments 为多序列比对结果，其中包含单拷贝基因supergene矩阵。