1. 前言
系统发育分析是我硕士的课题，也是我生信入门的第一课。如今Phylosuite已经可以完成系统发育分析的大部分内容，所以这一块的内容说难不难，但是要想做好也并不简单。就单单鉴定OG（Orthogroup Gene）这一步来说，从初期的Orthofinder再到后来的Hasmtr和Y & S，鉴定的手段不断完善，得到的信息也越来越多，随后的建树也逐渐从单基因树过渡到多基因树，直至现在的物种树。虽然分析的手段不短完善，但是数据集的选择在很大程度上会影响最后的结果，随着基因组数据的不短扩充，系统发育分析也应该逐渐发展至下一个维度，也就是基于全基因组数据进行系统发育研究。SOI就是一款专门为基因组数据做系统发育分析开发的软件，目前来看算法非常前沿，结果相对可靠。
2. 软件介绍
SOI是由昆植所的张仁纲老师开发的，张老师开发的软件最近我用的非常多吗，质量上乘，并且Github的问题回答速度非常及时。SOI的核心算法就是将共线性区块和直系同源基因进行了良好的结合，根据一个共线性区块中的直系同源基因的占比来评估这个共线性区块的同源性，很好的将共线性区块（定位）和同源基因（同源性）的优点结合。
3.软件下载
git clone https://github.com/zhangrengang/orthoindex.git
cd orthoindex
mamba env create -f OrthoIndex.yaml
mamba activate OrthoIndex
python3 setup.py install
soi -h
4. 输入文件
输入文件就是主要是4个文件：（1）每个物种的GFF文件（2）PEP文件（3）CDS文件（4）染色体长度文件。其实也就是WGDI的输入文件（如果不知道WGDI的输入文件格式，请自行去学习WGDI官方文档），并且SOI是对WGDI的格式完全兼容的。需要注意的是：the GENE ID is needed to label with SPECIES ID (e.g. Angelica_sinensis|AS01G00001) for uniqueness and compatibility (legacy from OrthoMCL). The CHROMosome ID should also be unique (e.g. As1, As2) to avoid conflicts (legacy from MCscanX). The IDs can be easily labeled for your own data.也就是每个物种的基因名要加上物种信息，张老师也提供了下面的脚本对GFF、CDS和PEP文件进行重命名。
SP=Angelica_sinensis
# using OrthoMCL command for fasta files:
orthomclAdjustFasta $SP $SP.pep 1
# or using sed:
sed ‘s/>/>’$SP’|/’ $SP.pep > $SP.fasta # always using the separator ‘|’
SP2=As
# using awk for MCscanX/WGDI gff files or JCVI bed files:
cat $SP.gff0 | awk -v sp=$SP -v OFS=”\t” ‘{$2=sp”|”$2;print $0}’ | perl -pe ‘s/^\D+/’$SP2’/’ > $SP.gff
cat $SP.bed0 | awk -v sp=$SP -v OFS=”\t” ‘{$4=sp”|”$4;print $0}’ | perl -pe ‘s/^\D+/’$SP2’/’ > $SP.bed
# 示例数据
# Gff
Artemisia_annua1 Artemisia_annua|ac1schr1g00001 714305 721459 – 1 chr1g00000451
Artemisia_annua1 Artemisia_annua|ac1schr1g00002 763149 767688 + 2 chr1g00000511
Artemisia_annua1 Artemisia_annua|ac1schr1g00003 789279 791221 – 3 chr1g00000551
Artemisia_annua1 Artemisia_annua|ac1schr1g00004 1061154 1062734 + 4 chr1g00000811
Artemisia_annua1 Artemisia_annua|ac1schr1g00005 1116199 1121766 + 5 chr1g00000871
Artemisia_annua1 Artemisia_annua|ac1schr1g00006 1122409 1123990 + 6 chr1g00000881
Artemisia_annua1 Artemisia_annua|ac1schr1g00007 1124611 1125263 + 7 chr1g00000891
Artemisia_annua1 Artemisia_annua|ac1schr1g00008 1126419 1130152 – 8 chr1g00000911
Artemisia_annua1 Artemisia_annua|ac1schr1g00009 1143062 1144152 – 9 chr1g00000931
Artemisia_annua1 Artemisia_annua|ac1schr1g00010 1158015 1161127 – 10 chr1g00000961
# PEP
>Artemisia_annua|ac1schr1g00002 chr1g00000511
MMGRSSSSDENGLKKGPWTPEEDQKLIDYIERNGHGSWRALPSLAGLNRCGKSCRLRWTN
YLRPDIKRGNFTEDEEKLIIHLHSHLGNKWSAIATHLPGRTDNEIKNFWNTHLKKKLLQM
GIDPVTHQPRTDHLDILAHLPQLLAIAANNLGATSLDLNIRDLLNPSLFDMNITALMSHS
DATQIAKFQLLQNILQVLSPNAPTNTTSIPQNIDPFYAMQLSDYLKINPQHLQNLQELAM
NVTPSASSLPPFLTHDHFRSPTYQTLNLDNDQASKGFSEQDINRAHGVTVGKEEYLGFES
VPALVPASPEHLDQSYSGKGQSPGANVNIGLGNSCNLKQKQSEHLSSSDPTSTATTIDAN
WPEFMDDEASGCYWKEIIERSSPPWS
>Artemisia_annua|ac1schr1g00003 chr1g00000551
MLIWIFRTKSMECASALHQPDGVGLQGSDIGSLSEINGVSIALVARTGILS
# CDS
>Artemisia_annua|ac1schr1g00001 chr1g00000451
ATGTCAGCAATTGTTTACTATGATGATATGCGAGTGAGCCCAGAGGTCATTGATCCTCCT
CAAATCGAAGATGCAATGGATAAAGAACATCTACATGTCCATGACCATGACCCAGCTCAA
AATGCTACAAAACCAAATCAACTTGCCACTAGCAGTGTTCGGGACCTATTGGAGTGTCCT
GTTTGCTTGAATGCCATGTACCCTCCTATCCATCAGTGCTCCAATGGTCACACATTATGT
TCTGGTTGCAAACCTCGGGTGCACAATCGATGTCCAACATGCAGACATGAACTAGGAAAC
ATTAGATGTCTTGCACTCGAGAAGGTTGCTGCGTCCCTTGAACTTCCTTGTAAATATCAG
AGTTACGGGTGTGTGGGAATATTTCCTTATTACAACAAATTAAAACATGAATCTCAATGT
GTTTTCCGACCCTATAATTGTCCTTATGCGGGTTCAGAATGTACGGTTATTGGTGACATC
CCTCATCTTGTGGCCCATTTAAAAGATGGCCATAAAGTCGACATGCATAATGGGAGTACA
# 染色体信息
Artemisia_annua1 238524853 8098
Artemisia_annua2 204548056 6575
Artemisia_annua3 194009994 6900
Artemisia_annua4 183757890 5631
Artemisia_annua5 180395630 5602
Artemisia_annua6 173031239 5290
Artemisia_annua7 158731363 5035
Artemisia_annua8 133201674 2627
Artemisia_annua9 81773804 1822
5. 运行
因为我是为了做系统发育分析，所以只提供系统发育分析的流程。
# 运行Orthocinder
orthofinder -f OrthoFinder/ -M msa -t 60
# 运行diamond，生成WGDI配制文件，运行WGDI，生产ctl文件
cat Species_pairs.txt | while read LINE
do
arr=($LINE)
SP1=${arr[0]}
SP2=${arr[1]}
prefix=”${SP1}_${SP2}”
conf=”${prefix}.conf”
echo “[collinearity]
gff1 = ${SP1}.gff
gff2 = ${SP2}.gff
lens1 = ${SP1}_length.txt
lens2 = ${SP2}_length.txt
blast = ${prefix}.blast
blast_reverse = no
multiple = 1
process = 50
evalue = 1e-5
score = 100
grading = 50,40,25
mg = 25,25
pvalue = 1
repeat_number = 10
positon = order
savefile = ${prefix}.icl
[ks]
cds_file = Conbine.cds
pep_file = Conbine.pep
align_software = muscle
pairs_file = ${prefix}.icl
ks_file = ${prefix}.ks” > “$conf”
# blast
diamond blastp -q ${SP1}.pep -d ${SP2}.pep -o ${prefix}.blast –more-sensitive -p 10 –quiet -e 0.001
# call synteny and calculate Ks
wgdi -icl $conf
# make ctl
echo 1500 >> $prefix.ctl
echo 1500 >> $prefix.ctl
cut -f 1 ${SP1}.gff| sort -u | paste -sd “,” >> $prefix.ctl
cut -f 1 ${SP2}.gff| sort -u | paste -sd “,” >> $prefix.ctl
done
这是示例文本，单击 “编辑” 按钮更改此文本。