參考：
利用orthofinder尋找單拷貝基因構(gòu)建系統(tǒng)發(fā)育樹 - 簡(jiǎn)書 (jianshu.com)
創(chuàng)建conda環(huán)境陈瘦，我個(gè)人命名為orthofinder目养，在該環(huán)境下用conda就可安裝orthofinder呐籽，muscle，mafft裕照，seqkit，phyml軟件壁熄。
其他軟件prottest,RAxML可以官網(wǎng)下載耻台，安裝就看README或者INSTALL。比如prottest紫皇，最后還要把phyml軟鏈到bin下慰安，比如RAxML，make命令選一個(gè)跟自己系統(tǒng)匹配的運(yùn)行一個(gè)就夠了聪铺，不用運(yùn)行全部四個(gè)化焕。
1 Othofinder運(yùn)行
首先激活環(huán)境，然后可以后臺(tái)運(yùn)行以下命令
orthofinder -f your/pep/fa/file/
我沒加線程參數(shù)铃剔，所以十個(gè)物種的大概運(yùn)行2-4個(gè)小時(shí)撒桨。生成文件中會(huì)有 Single_Copy_Orthologue_Sequences文件夾。
键兜。
2.對(duì)生成的單拷貝同源序列進(jìn)行批量處理凤类，比對(duì)。
在Single_Copy_Orthologue_Sequences文件夾下運(yùn)行
$ vi 01.muscle.sh
#!bin/bash
for i in *.fa
do muscle -in $i -out $i.1 #新版本的muscle命令變了普气。
#或者do muscle -align $i -output $i.1
done
運(yùn)行后多出一些后綴為.1的文件谜疤。*.1文件是比對(duì)好的序列文件。
3.提取保守序列
$ vi 02.gblock.sh
conda activate orthofinder
cd PATH/TO/Single_Copy_Orthologue_Sequences
for i in *.1
do Gblocks $i -b4=5 -b5=h -t=p -e=.2
done
會(huì)生成一堆后綴是*.1.2的文件。
4.把序列排序夷磕，合并
$ vi 03.seqkit.sh
conda activate orthofinder
cd OrthoFinder/Results_May20/Single_Copy_Orthologue_Sequences
for i in *.2
do seqkit sort $i >$i.3
seqkit seq $i.3 -w 0 > $i.3.4
done
這里參考教程新建了文件夾new 把.4后綴的文件移到文件夾里然后
$ mkdir new
$ mv *.4 new/
$ cd new
$ paste -d " " *.4 > all.fa
$ sed -i "s\ \\g" all.fa
得到all.fa文件履肃。
5.處理數(shù)據(jù)，將fa文件轉(zhuǎn)為phy格式
vi 04.fa2phy.py
#!usr/bin/python
import re
with open('all.fa', 'r') as fin:
sequences = [(m.group(1), ''.join(m.group(2).split()))
for m in re.finditer(r'(?m)^>([^ \n]+)[^\n]*([^>]*)', fin.read())]
with open('all.phy', 'w') as fout:
fout.write('%d %d\n' % (len(sequences), len(sequences[0][1])))
for item in sequences:
fout.write('%-20s %s\n' % item)
$ java -jar /your/software/to/prottest3-master/dist/prottest-3.4.2.jar -i all.phy -all-distributions -F -AIC -BIC -tc 0.5 -threads 24 -o prottest.out
得到蛋白模型坐桩。
6.RAxML構(gòu)樹
$ vi RAxML.sh
/your/path/to/01.Software/standard-RAxML-master/raxmlHPC-PTHREADS-SSE3 -T 16 -f a -x 123 -p 123 -N 1000 -m PROTGAMMAIJTTF -k -O -o Ensete_ventricosum,Musa_a,Musa_b \
-n all.rename.tre -s all_rename.fa
qsub提交任務(wù)就可尺棋。十個(gè)物種大概運(yùn)行1小時(shí)。
這里-o選項(xiàng)是定根撕攒，我的外群是芭蕉陡鹃。如果不定外群，結(jié)果會(huì)很奇怪抖坪。
還有就是萍鲸。-s參數(shù)后接的fa文件，其實(shí)是把每個(gè)物種的單拷貝同源序列合并成超長(zhǎng)一條擦俐。最好把ID改成相應(yīng)的物種脊阴。不然ID太長(zhǎng)也會(huì)報(bào)錯(cuò)。
然后我用的besttree復(fù)制到figtree軟件或者iTOL蚯瞧，就可以出圖了嘿期。
注：我只選了同一“目”下的一些物種，我要研究的五種物種的基因組埋合，加上近緣物種選了五種备徐，一共是十種。因?yàn)楹芙跛蹋詆blocks可以運(yùn)行蜜猾。因?yàn)槲蚁朐诟嗟奈锓N比如單子葉和雙子葉選了13種的時(shí)候gblocks好像選不出來保守基因，可能是單拷貝同源基因太少了振诬。（個(gè)人推測(cè)）
個(gè)人覺得這種單純的樹沒啥用蹭睡，后續(xù)還要加分化時(shí)間（r8s等）還有CAFE分析。
