一、phy文件的获得

phy文件由cds和pro序列比对得到，cds和pro需要保持id一致

（1）cds序列的获取

从数据库下载到总的cds序列后，需要对其进行相关处理

cds文件的基因id的转化

awk '/^>/ {if (match($0, /\\[locus_tag=([^]]+)\\]/, arr)) print ">" arr[1]; else print $0} !/^>/ {print}' input.fa > output.fa
##这里是保留locus_tag部分的ID，有些是protein_id则需把locus_tag换成protein_id

获取所有物种的cds序列汇集成all.fa

cat *cds > all.fa

pro序列，由orthofinder获得的单拷贝序列。

cds序列则是由pro序列id提取的，因此我们先要提取所有的单拷贝序列id然后提取响应的cds序列。

for file in *fa; do
	
	**base_name="${file%.fa}"**
				
	grep '>' $file | sed 's/>//g' > $base_name.id
				
	seqkit grep -f $base_name.id all.fa > $base_name.cds

done

对pro序列进行多序列比对并比对到cds序列上

for file in *.fa; do
    # 提取文件名前面的部分（去掉 .aln 后缀）
	base_name="${file%.fa}"
    
	mafft --maxiterate 1000 --localpair "$file" > "${base_name}.aln"

	pal2nal.pl "${base_name}.aln" "${base_name}.cds" -nogap -nomismatch -output fasta > "${base_name}.phy"

    **
    done

保证phy序列跟树文件的id一致

for i in *phy ; do sed -i -E 's/^(>...).*/\\1/' $i ; done

二、树文件的获得

可直接使用基因家族扩张与收缩中的树文件，想要检验哪一支，是不是受到了正选择，我们需要在树结构对应的位置上加上#1，表示foreground，其他的为background。

14      1
(Cre,(((Mco,Cva),(Psp,(Hsp,(Apr,Pwi)))),((Ebi,(Cpr,Csu #1)),(Agl,(Tsp4,(Tsp5,Tsp6)))))'B(.662,1)');

三、批量进行正向选择鉴定