官網(wǎng)：http://www.hyphy.org/
http://www.bork.embl.de/pal2nal/
參考：[hyphy]http://www.itdecent.cn/p/2e8f7f7d545a
[pal2nal]https://blog.csdn.net/qq_50637636/article/details/120226785?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2
1 所需要的
- 基因組序列（案例是野生稻和栽培稻）
- GFF3基因組注釋文件
- 蛋白質(zhì)序列
- CDS序列（有ATG起始密碼子，沒有TAG, TAA, TGA終止密碼子）
- 樹文件（.nwick格式文件）
- 在python2.7的conda環(huán)境下使用
2 步驟
- 得到野生稻和栽培稻的CDS文件
- 得到野生稻和栽培稻的蛋白質(zhì)文件
- 比對蛋白質(zhì)文件
- 用pal2nal軟件將已經(jīng)比對好的蛋白及其對應(yīng)的DNA多序列轉(zhuǎn)化為密碼子比對的程序。
舉例：
本來用CDS序列直接去比對，順序可能是：-AGTTG-GGAATAAT-TT--TAT-
但通過pal2nal后
蛋白質(zhì)序列：-SWE--LF-Y-
將CDS反向比對回去就變成：-AGTTGGGAA--TAATTT-TAT- （三個三個密碼子成對出現(xiàn)）
- 構(gòu)建系統(tǒng)發(fā)育樹
- 使用hyphy計算選擇壓力
2 安裝
pal2nal的安裝
conda create -n pal2nal
conda activate pal2nal
conda install pal2nal
hyphy的安裝
conda create -n python2.7 python=2.7 # 創(chuàng)建環(huán)境
conda activate python2.7 # 進入環(huán)境
conda install hyphy # 安裝hyphy
3 序列的準備
因為該軟件要求蛋白質(zhì)序列和CDS序列要保持一致性，所以我用gffread軟件從基因組文件中提取出CDS文件，再翻譯成蛋白質(zhì)文件。
gffread test.gff3 -g test.genomic.fasta -x test.cds -y test.pep
我之前試過直接用網(wǎng)上下的數(shù)據(jù)，會報錯，數(shù)據(jù)比較混亂。如下圖，第一個密碼子AGC應(yīng)該翻譯成S（絲氨酸），但卻翻譯成了E。
2fee0f2b9cb9278ade538f9ac30c514.png
4 序列的比對
4.1 蛋白質(zhì)序列的比對
mafft --auto --thread 10 test.pep > test.pep.aln # 線程數(shù)可以自己看著改
4.2 pal2nal的使用
perl /Path/To/pal2nal.pl pep.fas nuc.fas -output fasta
# 輸入文件，兩個，分別是比對好的PEP和原始的CDS序列
# -output 輸出文件格式，默認為clustal
4.3 過濾一些gap（可選）
為什么要過濾？
- 前后會產(chǎn)生很多gap，對選擇壓力的計算沒有用，徒增計算壓力，也會浪費時間；
- 有些序列可能會全都是gap，在后面構(gòu)樹的時候也會被濾掉的。這樣的話，在最后hyphy計算選擇壓力的時候，就會出現(xiàn)樹的支（即序列數(shù)）比文件的序列數(shù)少1的錯誤。
怎么過濾？
trimal -in test.pal2nal.fasta -out test.pal2nal.gappyout.fasta -gappyout
5 構(gòu)建系統(tǒng)發(fā)育樹
iqtree -s test.pal2nal.gappyout.fasta -m MFP -bb 1000 -nt AUTO --prefix test
# -s 輸入比對序列
# -m MFP 選擇最佳模型后構(gòu)樹
# -bb 1000 快速自展1000次
# 因為快速自展1000次支持值結(jié)果可能會偏大，所以時間比較寬裕的話可以加上-bnni
# -nt AUTO 線程數(shù)，給iqtree自行安排和選擇
# --prefix 輸出文件前綴
6 hyphy的使用
- 選擇方法
使用hyphy -i，一步一步選擇合適你目的的方法 - 開始計算(我這里用的方法是absrel)
hyphy absrel --alignment test.pal2nal.gappyout.fasta -tree test.treefile CPU=10 > hyphy.out
# 如果是新版本的話，hyphy會自動根據(jù)需求采用多線程
# 你可以加上foreground或者branchs信息，詳細的使用可以看官網(wǎng)。
- 輸出結(jié)果
Hyphy運行的時候，默認打印到屏幕上的結(jié)果是以markdown格式輸出的，而保存到本地文件的結(jié)果是以json格式輸出的(json格式可以很方便的用python的json模塊提取各種信息，例如pvalue和正選擇位點，適合多個任務(wù)批量操作)。默認是輸出到和多序列比對文件相同的文件夾，可以用--output來改變輸出位置。 - 結(jié)果可視化
可以去官網(wǎng)http://vision.hyphy.org/來可視化輸出結(jié)果。
