工作原因剛接觸病毒基因組，發(fā)現(xiàn)跟細(xì)菌不太一樣。讀到一個(gè)帖子，詳細(xì)說(shuō)明一些位點(diǎn)可信度的問(wèn)題，以及是否應(yīng)該過(guò)濾一些位點(diǎn)從而達(dá)到更真實(shí)的系統(tǒng)發(fā)育樹(shù)。說(shuō)的很詳細(xì)，鏈接如下
https://virological.org/t/issues-with-sars-cov-2-sequencing-data/473
但是真實(shí)操作的時(shí)候，有些是不能?chē)?yán)格按照這個(gè)方法做的。秉著構(gòu)建更合理的樹(shù)的原則（本人并不專(zhuān)業(yè)，是自學(xué)的而且剛?cè)腴T(mén)的），自己搞了一套方法。如果后面還要繼續(xù)做病毒，可能還會(huì)完善；之后就不再做這個(gè)方向了，也就不再完善了。以上，因此此方法僅供參考，同時(shí)歡迎提出更好的方法，糾正我的錯(cuò)誤。

1. 下載數(shù)據(jù)
GISAID下載即可，鼠標(biāo)點(diǎn)點(diǎn)點(diǎn)，很簡(jiǎn)單，不再贅述
GISAID網(wǎng)站：[https://www.gisaid.org/](https://www.gisaid.org/)
先需要注冊(cè)

2. 根據(jù)長(zhǎng)度過(guò)濾樣本
個(gè)人標(biāo)準(zhǔn)是完整基因組，長(zhǎng)度大于等于29000.同時(shí)也看了下長(zhǎng)度大于等于29400的，相差2個(gè)樣本而已，就按照自己喜歡的29000來(lái)了。

3. mafft多序列比對(duì)
/path/to/software/mafft-7.453-with-extensions/core/mafft --thread -1 --auto --keeplength --addfragments world_20200420.fasta EPI_ISL_402125.fasta > world_20200420.mafft.fasta

或者
數(shù)據(jù)多的時(shí)候這一步特別慢，嘗試用blast來(lái)做
#EPI_ISL_402125.fasta = wuhan-hu-1
makeblastdb -in world_20200420.fasta -dbtype nucl -out world_20200420.fasta
blastn -db world_20200420.fasta -query EPI_ISL_402125.fasta -out blast -num_alignments 30000 -num_descriptions 30000
/path/to/software/mafft-7.453-with-extensions/core/mview -in blast -out fasta blast -pcid reference -top all -moltype dna > blast.mview.fa

#去掉名字里的特殊字符（自己寫(xiě)的腳本）
perl rename.pl blast.mview.fa > blast.mview.rename.fa
#去掉不想要的序列（自己寫(xiě)的腳本，現(xiàn)在還不需要，后面需要去掉序列時(shí)候再去）
perl remove.pl blast.mview.rename.fa 2019 > blast.mview.rename.rmout.fa
#序列變成一行（自己寫(xiě)的腳本，個(gè)人習(xí)慣）
perl fa.pl blast.mview.rename.rmout.fa > blast.mview.rename.rmout.straight.fa

4. 腳本提取consensus序列上的variant
統(tǒng)計(jì)
1）哪些位點(diǎn)覆蓋度不足90%（例如，1000條序列，只有800個(gè)覆蓋了11011位點(diǎn)，該位點(diǎn)覆蓋度不足）
2）哪些位點(diǎn)突變類(lèi)型多，例如C->T占30%  C->A占30%  C->G還有20%  雜合度非常高
3）兩兩比較樣本，某些和其他樣本相比variant非常多的
4）建樹(shù)后，枝長(zhǎng)特別長(zhǎng)，還跟別的都不聚在一起的
這些位點(diǎn)/樣本都去掉。注：第四條不是此時(shí)去掉的，建樹(shù)后才去
修整完，文件是test.mafft.fasta

5. 建樹(shù)
IQTree，目前被偏愛(ài)的建樹(shù)寵兒
輸入文件可以是fasta，phylip等

#檢查best-fit model
iqtree -s test.mafft.fasta -m MF
#如果是snp串聯(lián)的文件
iqtree -s test.mafft.fasta -m MFP+ASC
#產(chǎn)生的log文件里有最優(yōu)模型，我測(cè)試了兩個(gè)都是GTR+F+R2

#建樹(shù)
iqtree -s test.mafft.fasta -m GTR+F+R2 -bb 1000 -bnni -nt AUTO -redo