基因组上有很多变异，不同变异的效果不同。
除了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表示是三字母氨基酸缩写，需要转为单字母氨基酸缩写。

预测变异的效果