工作原因刚接触病毒基因组，发现跟细菌不太一样。读到一个帖子，详细说明一些位点可信度的问题，以及是否应该过滤一些位点从而达到更真实的系统发育树。说的很详细，链接如下
https://virological.org/t/issues-with-sars-cov-2-sequencing-data/473
但是真实操作的时候，有些是不能严格按照这个方法做的。秉着构建更合理的树的原则（本人并不专业，是自学的而且刚入门的），自己搞了一套方法。如果后面还要继续做病毒，可能还会完善；之后就不再做这个方向了，也就不再完善了。以上，因此此方法仅供参考，同时欢迎提出更好的方法，纠正我的错误。
1. 下载数据
GISAID下载即可，鼠标点点点，很简单，不再赘述
GISAID网站：[https://www.gisaid.org/](https://www.gisaid.org/)
先需要注册
2. 根据长度过滤样本
个人标准是完整基因组，长度大于等于29000.同时也看了下长度大于等于29400的，相差2个样本而已，就按照自己喜欢的29000来了。
3. mafft多序列比对
/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
或者
数据多的时候这一步特别慢，尝试用blast来做
#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
#去掉名字里的特殊字符（自己写的脚本）
perl rename.pl blast.mview.fa > blast.mview.rename.fa
#去掉不想要的序列（自己写的脚本，现在还不需要，后面需要去掉序列时候再去）
perl remove.pl blast.mview.rename.fa 2019 > blast.mview.rename.rmout.fa
#序列变成一行（自己写的脚本，个人习惯）
perl fa.pl blast.mview.rename.rmout.fa > blast.mview.rename.rmout.straight.fa
4. 脚本提取consensus序列上的variant
统计
1）哪些位点覆盖度不足90%（例如，1000条序列，只有800个覆盖了11011位点，该位点覆盖度不足）
2）哪些位点突变类型多，例如C->T占30% C->A占30% C->G还有20% 杂合度非常高
3）两两比较样本，某些和其他样本相比variant非常多的
4）建树后，枝长特别长，还跟别的都不聚在一起的
这些位点/样本都去掉。注：第四条不是此时去掉的，建树后才去
修整完，文件是test.mafft.fasta
5. 建树
IQTree，目前被偏爱的建树宠儿
输入文件可以是fasta，phylip等
#检查best-fit model
iqtree -s test.mafft.fasta -m MF
#如果是snp串联的文件
iqtree -s test.mafft.fasta -m MFP+ASC
#产生的log文件里有最优模型，我测试了两个都是GTR+F+R2
#建树
iqtree -s test.mafft.fasta -m GTR+F+R2 -bb 1000 -bnni -nt AUTO -redo
网友评论