0.写在前面
SIFT本身内置了许多常见物种的评分数据库，在尝试构建自己的库之前，先在官网找找有没有自己需要的。但是有些库上次更新时间有点吓人。
SIFT的下载界面太难找了，官网：SIFT - 预测非同义/错义变体的影响 (a-star.edu.sg)。似乎服务器不太稳定，多刷新几次试试
Chatgpt机翻
SIFT4G_PATH：指向可执行文件sift4g的路径，即需求软件
PROTEIN_DB：指向蛋白质数据库（.fa或.fasta）的路径。建议使用UniProt 90
PARENT_DIR：生成所有输出的路径。用户必须有写入权限。
ORG：物种名称，例如homo_sapiens, rat等。应该是一个单词，没有特殊字符。
ORG_VERSION：最终的预测数据库将在<PARENT_DIR>/<ORG_VERSION>中
GENE_DOWNLOAD_SITE (仅限Ensembl)：下载基因注释（gtf）文件的网址
PEP_FILE (仅限Ensembl)：下载蛋白质序列文件的网址。这个文件是可选的，用于质量控制。
CHR_DOWNLOAD_SITE (仅限Ensembl)：下载基因组序列（.fa）的网址
DBSNP_ORGANISM_DOWNLOAD_SITE (仅限Ensembl)：下载dbSNP注释的网址文件夹（可选）
DBSNP_VCF_FILE (可选)：一个*.vcf.gz文件，用于用dbSNP rs id注释变异
GENETIC_CODE_TABLE：用于将DNA序列翻译成蛋白质的遗传密码（基于NCBI的整数值）*
GENETIC_CODE_TABLENAME (不使用)：字符串，用于提醒用户使用了哪种GENETIC_CODE_TABLE
MITO_GENETIC_CODE_TABLE：命名为Mt, chrM, 或Mito的fasta序列将使用这种遗传密码。（这是用于线粒体序列）。要禁用或编辑此功能，请编辑dna_protein_subs.pl中的函数chr_is_mito
MITO_GENETIC_CODE_TABLENAME (不使用)：字符串，用于提醒用户使用了哪种MITO_GENETIC_CODE_TABLE
PLASTID_GENETIC_CODE_TABLE：命名为Pt, PT, chrPt, chrPT, 或chloroplast的fasta序列将使用这种遗传密码。（这是用于叶绿体）。要禁用或编辑此功能，请编辑dna_protein_subs.pl中的函数chr_is_plastid
1.SNP
需要多个软件：
https://github.com/pauline-ng/SIFT4G_Create_Genomic_DB （构建数据库）
https://github.com/rvaser/sift4g （需求软件）
https://github.com/pauline-ng/SIFT4G_Annotator （注释）
配置要求：g++ (4.9+)，GNU Make，nvcc (2.*+) (optional)）
1.构建数据库
1.下载
git clone --recursive https://github.com/rvaser/sift4g.git sift4g
cd sift4g/
make #请确认具有7.3以上的gcc和g++，否则后续建库会报错（what(): regex_error）。这一步报错了就检查一下配置需求
git clone https://github.com/pauline-ng/SIFT4G_Create_Genomic_DB.git scripts_to_build_SIFT_db
蛋白质库作者建议是用UniProt 90，下载链接：https://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref90/uniref90.fasta.gz
文件很大（39G），国内建议使用镜像站点：ftp://ftp.ebi.ac.uk/pub/databases/uniprot/uniref/uniref90/
下载完还要解压缩，解压后75G，整了一个小时
1.小测试
强烈建议用这个小基因组先测试一下，运行过程中会出现各种报错，需求库，见招拆招即可。小数据跑起来快一点，哺乳动物基因组一般是24-48小时
cd test_files/candidatus_carsonella_ruddii_pv_config.txt
编辑配置文件以设置<PARENT_DIR>、<SIFT4G_PATH><PROTEIN_DB> ，需要完整路径
这里先用ensembl的物种蛋白质数据库
运行两次，一次使用该物种的蛋白质序列（配置环境，几分钟跑完），一次使用uniref90.fasta（确定能否正确运行，约70分钟），具体运行结果差异见文末
在软件主目录下
perl make-SIFT-db-all.pl -config test_files/candidatus_carsonella_ruddii_pv_config.txt --ensembl_download #
结果存储在：ASM1036v1.34文件夹下
2.ensembl
最简单便捷的方法
使用test_files/saccharomyces_cerevisiae-template.txt作为模版
设置网址链接到Ensembl的基因组和基因注释文件：GENE_DOWNLOAD_SITE, PEP_FILE, CHR_DOWNLOAD_SITE 可选：DBSNP_ORGANISM_DOWNLOAD_SITE
设置输出路径：PARENT_DIR, ORG, ORG_VERSION
设置遗传密码表：GENETIC_CODE_TABLE, MITO_GENETIC_CODE_TABLE
如果物种是脊椎动物，那么
GENETIC_CODE_TABLE=1
GENETIC_CODE_TABLENAME=Standard
MITO_GENETIC_CODE_TABLE=2
MITO_GENETIC_CODE_TABLENAME=Vertebrate Mitochondrial
对于非哺乳动物，可以试试参照官方对应的数据库，里面会有一份配置表。
4.设置SIFT4G路径：SIFT4G_PATH, PROTEIN_DB
perl make-SIFT-db-all.pl -config <config file> --ensembl_download
跳过下述根据gtf和gff配置数据库，开始检测数据库能否正常运行。下同
3.使用本地的fa和gtf文件构建数据库
0.对文件的预处理
如果你和我一样是从NCBI下载的数据，请注意，有的基因组fa文件和注释文件里面的染色体号不一致
pattern=$(sed 's/\(.*\)\t\(.*\)/s\/\1\/\2\/g;/g' CM-NC.table) #其中table文件是制表符分隔的两列文件，第一列是旧内容，第二列是新内容
sed "$pattern" test.fa > new.fa #运行起来也略慢
创建一个新文件夹，将 test_files/homo_sapiens-test.txt 复制过来作为模板
设置：<PARENT_DIR>、<ORG>、<ORG_VERSION>、<SIFT4G_PATH>，<PROTEIN_DB>
检查GENETIC_CODE_TABLE和MITO_GENETIC_CODE_TABLE设置正确。同上
可选设置：<DBSNP_VCF_FILE>
2.将各种文件移动到合适的位置
mkdir <PARENT_DIR> #即输出目录
mkdir <PARENT_DIR>/gene-annotation-src
mkdir <PARENT_DIR>/chr-src
mkdir <PARENT_DIR>/dbSNP
mv *.fa.gz <PARENT_DIR>/chr-src #基因组fa文件
mv *.gtf.gz <PARENT_DIR>/gene-annotation-src #gtf注释文件
mv *.vcf.gz <PARENT_DIR>/dbSNP #可选，dbsnp数据库，还是要注意染色体号的问题
mv *.pep.all.fa.gz <PARENT_DIR>/gene-annotation-src #可选：将压缩的蛋白质文件（.pep.all.fa.gz）放入<PARENT_DIR>/gene-annotation-src。此文件用于检查。这个命名方式应该是ensembl的物种的蛋白质序列。
上述均为压缩文件
3.运行
在运行之前，务必确认：
已成功运行测试物种
蛋白质数据库文件已经解压缩
基因组fa文件的染色体号与gtf文件一致
gtf文件，fa文件已压缩
dbsnp数据库（如果有）文件染色体号与上述一致
24小时后服务器有足够大的内存和CPU空间，CPU至少800，内存30G（uniref90）
perl make-SIFT-db-all.pl -config <config_file> #主目录下
在运行过程中，除非程序主动停止或显示All done！，尽量不要因为看到几个error就将其停止。
当显示All done!时，数据库构建完成
2.70GB基因组
使用NCBI下载的物种蛋白质库运行时间：31号00:10:39→(processing database 环节cpu占用升到了800) →（收尾环节猛报ValueError）→All done！1号11:21:56
使用uniref90.fasta运行时间：31号20:53：10→1号下午18点左右，内存占用飙升（即上述cpu占用上升环节），只能等有空再跑了。
这里服务器还有0.9G的空余，不太清楚是程序预留的空间，还是占用只有这么多，保险起见，还是终止了
重新运行：4号14:56:30→(确定了，这程序是有多少吃多少，最高达到过50G，还没摸到上限。常态20g，24个小时）→（还是有ValueError）→7号5:04:35
uniref90.fasta运行时间过长，占用资源过多，因此我有个大胆的想法，为什么不把每条染色体的注释文件拆分出来单独运行呢
结果见文末
其实这里我有一个疑惑，众所周知啊，NCBI的转录组id，蛋白质id，外显子id都是有一套独特命名规则，使用NCBI的各种id在uniref90.fasta中也检索不到对应的。但是这里对gtf文件并不需要进行特殊的id转换，也能正常构建数据库，不太清楚它的内部逻辑。（Candidatus_carsonella_ruddii的NCBIgtf注释文件并不能用于构建数据库，可能是NCBI上该物种的gtf文件不完善？）
此外，我还建议，从gtf文件中提取一条染色体乃至一小部分基因进行测试，以避免运行24个小时之后，报错的窘境。
如：
grep 'NC_040268.1' genomic.gtf >new.gtf
gzip new.gtf
再重命名一下genomic.gtf文件，破坏其后缀，运行即可
这里我用uniref90.fasta成功得到了chr.gz文件，但打分列有较多的NA，region文件是空的，且最后一步会出“key error”，换成NCBI下载的fa文件则没有问题。尽管高置信度的比例是一样的，均为56%。
4.使用本地的fa和gtf文件构建数据库
没实操，应该没问题
从https://github.com/gpertea/gffread下载gffread，
2.将gff文件转换为gtf文件
mv <gff3.gz> <PARENT_DIR>/gene-annotation-src
gunzip <*.gff3.gz>
gffread <*.gff3> -T -o [FILENAME].gene.gtf
gzip [FILENAME].gene.gtf
3.按上述步骤构建数据库即可
5.检测数据库是否正确配置
数据库存储在<PARENT_DIR>/<ORG_VERSION>中，这是在配置文件中设置的
每条染色体对应一个gz的文件
zcat <chr>.gz | more # does it look all right?
zcat <chr>.gz | grep CDS | more
SIFT评分位于该文件的10-12列，没有太多的NA即可
CHECK_GENES.LOG是统计文件，最后一行的百分比均较高，那么数据库就构建完成了
6.整理冗余文件
从注释流程来看，只需要数据库文件夹，其他的数据都可以删除，但是我认为适当的保留日志文件是有必要的，因此，建议删除：singleRecords（约100g大小）文件夹，SIFT_predictions文件夹（10g，不要打开，内部文件太多了）
希望删了之后，不会影响最后的注释。
2.注释
1.下载
wget https://github.com/pauline-ng/SIFT4G_Annotator/raw/master/SIFT4G_Annotator.jar
2.运行
java -jar <Path to SIFT4G_Annotator> -c -i <Path to input vcf file> -d <Path to SIFT4G 数据库目录> -r <Path to your result folder> -t
`-c` 表示使用坐标模式（Coordinate mode），即根据VCF文件中的染色体和位置信息来注释变异（必须的参数）
`-i` 表示输入vcf文件的路径
`-d` 表示SIFT4G数据库目录的路径
`-r` 表示结果文件夹的路径
`-t` 表示提取多个转录本的注释（Optional），可选参数
另外，如果是自主构建的数据库，确保存在非空.gz和区域文件，且数据库文件不能以“chr开头”
mv chr19.gz -> 19.gz; mv chr19.regions 19.regions; mv chr19_SIFTDB_stats.txt 19_SIFTDB_stats.txt
2.indel
找sift下载界面不小心找到的
划掉，怎么只有人类基因组
3.做个总结
可以看到最终结果，高置信度的比例只有55%，不管怎么说，都有点略低。用dbsnp库，蛋白质序列排列组合试试吧，再下一步我打算用ensembl的数据库做试试，看是不是NCBI的gtf文件的问题。
整体上来说，uniref90.fasta的数据结果还是要好一点的，尽管60也不是很高，那就大概率是gtf文件的问题，也没啥办法
官方的是14年构建的，不管是参考基因组还是注释都是远远落后于时代，但是有很高的高置信度比例，而且有线粒体的数据，这可能是gtf文件或蛋白质序列的问题
虽然前前后后花了我一个星期的时间，但其实整体软件运行下来，难度不是很大，最头疼的是构建一次数据库需要几十个小时。
