记录一下比较基因组一部分内容。

下载相关物种的ref_fa及gff

可从NCBI等数据库下载

去除gff文件中gene,cds,mRNA多余字符
sed -i 's/gene-//g' xxx.gff
sed -i 's/rna-//g' xxx.gff
sed -i 's/cds-//g' xxx.gff

提取其对应的cds和pep文件

利用软件得到其cds pep序列文件 这里用到了iTools_code gffread也可以
/data/01/user152/software/iTools_Code/iTools Fatools getCdsPep-Ref xxx.fna -Gff xxx.gff -OutPut xxx

对两个文件进行处理及过滤

基于gff3提取最长cds序列ID，并过滤gff3文件
perl gff_filter_longest.pl xxx.gff xxx_gene_mrna_cds.ids xxx_final.gff

提取最长cds ID列表
具体查看列数为2还是3
awk '{print $2}' xxx_gene_mrna_cds.ids > xxx_mRNA.id

基于最长cds ID信息提取cds和蛋白质序列文件
seqtk subseq xxx.cds.fa  xxx_mRNA.id > xxx.L.cds.fa
seqtk subseq xxx.pep.fa  xxx_mRNA.id > xxx.L.pep.fa

cds和pep不匹配的过滤；pep长度小于50aa的过滤；终止密码子不考虑；cds非3倍数过滤
perl PEP_CDS_Flt.pl xxx.L.cds.fa xxx.L.pep.fa
mv xxx.L.cds.fa.flt xxx.F.cds.fa
mv xxx.L.pep.fa.flt xxx.F.pep.fa

将序列前添加物种信息
sed 's/^>/>xxx|/g' xxx.F.cds.fa > xxx.cds.clean.fa
sed 's/^>/>xxx|/g' xxx.F.pep.fa > xxx.pep.clean.fa

选做下面两步
由于mRNA和cds 序列ID不一致，因此需要将gff文件中mRNA id 替换成cds ID
awk '{print $2"\t"$3}' xxx_gene_mrna_cds.ids > xxx_mRNA.mapid
 更改gff文件中mRNA id成cds id
perl map_gff_ids xxx_mRNA.mapid xxx_final.gff

运行orthofinder

将所有处理好的*.pep.clean.fa放在一个文件夹中
orthofinder -f pep-data -S diamond -M msa -T fasttree -t 15