基因组上有很多变异，不同变异的效果不同。
除了wet-lab，现在也有很多软件预测。接下来，我打算用SIFT4G、Polyphen-2、SnpEff和provean预测。
因为老板想法改变，剩下的两个软件就···不弄了
0. vcf文件准备
msa2vcf可以把多序列比对转为vcf文件。较其他工具，它可以处理gap，但无法转换为reference genome上的位置，需要自己写程序转换。
1. SnpEff
1.1 检查database是否正确
软件里存了很多database，大部分物种都在，用之前可以先检查一下感兴趣物种在不在。
# 列出已构建的database
java -jar snpEff.jar databases
不过，即使有，也需要检查database创建是否有误，比如genetic code是否正确
java -Xmx4g -jar snpEff.jar -v Lactobacillus_plantarum_gca_001005805
编码方式居然是standard code，事实上应该是bacteria code ，也就是genetic code 11。
1.2 build自己的database
具体参考：https://pcingola.github.io/SnpEff/se_buildingdb/
如何build，要根据自己手边的数据。接下来，我采取的是genbank方式。
分三步：1）下载genbank数据，2）配置文件，3）run
1.2.1 下载genbank数据
拿着chromosome或scaffold的编号，去ncbi的nucleotide库中搜，下载数据。
下完之后，将所有chromosome的cat起来，命名为genes.gbk。
一定得叫genes.gbk
1.2.2 配置文件
在snpEff.config中添加配置信息：
# Lactobacillus plantarum ps128
Lactobacillus_plantarum_ps128.genome : Lactobacillus plantarum
Lactobacillus_plantarum_ps128.chromosomes : NZ_LBHS01000003.1, NZ_LBHS01000004.1, NZ_LBHS01000008.1, NZ_LBHS01000009.1, NZ_LBHS01000010.1, NZ_LBHS01000011.1, NZ_LBHS01000005.1, NZ_LBHS01000006.1, NZ_LBHS01000007.1, NZ_LBHS01000001.1, NZ_LBHS01000002.1
Lactobacillus_plantarum_ps128.NZ_LBHS01000003.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000004.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000008.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000009.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000010.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000011.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000005.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000006.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000007.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000001.1.codonTable : Bacterial_and_Plant_Plastid
Lactobacillus_plantarum_ps128.NZ_LBHS01000002.1.codonTable : Bacterial_and_Plant_Plastid
每个chromosome都要弄，chromosome的id不能搞错。比如NZ_LBHS01000003.1不能写成NZ_LBHS01000003。
在snpEff.config目录下做以下操作:
mkdir -p data/Lactobacillus_plantarum_ps128
mv genes.gbk data/Lactobacillus_plantarum_ps128/
1.2.3.run
java -jar snpEff.jar build -genbank -v Lactobacillus_plantarum_ps128
检查日志文件，是否有报错。
如果没有什么问题，在data/Lactobacillus_plantarum_ps128文件夹下会出现snpEffectPredictor.bin文件，创建成功。
1.3 预测
java -Xmx4g -jar snpEff.jar -v Lactobacillus_plantarum_ps128 -ud 0 0_align.formatted.vcf | ./scripts/vcfInfoOnePerLine.pl> 0_align.prediction.vcf
-Xmx4g：给程序分配4G内存，视自己情况而定
-ud: 设置upstream和downstream interval size。如果变异与gene A的距离小于interval size，预测时会包括对gene A的影响。
2. Provean
注意点1：
与SnpEff相比，Provean只能预测部分蛋白质编码基因的variation。建议Provean放在SnpEff之后，这样子可以借用SnpEff的结果作为Provean的输入，比较简单。
提取SnpEff结果中的非intergenic region的变异，作为Provean输入数据。
注意点2：
因为研究物种是细菌，所以只能选用PROVEAN Protein工具，缺点是不能批量，每次只能预测一个蛋白质序列。
From Provean Home
2.1 预测
为了方便整理，我为每个蛋白质都建了一个文件夹，里面放了两个文件：
- sequence.fa：存储一条蛋白质序列
- var：存放该蛋白质上的变异（必须是HGVS表示，SnpEff结果文件中有）
provean.sh -q sequence.fa -v var
SnpEff中的HGVS表示是三字母氨基酸缩写，需要转为单字母氨基酸缩写。
网友评论