SIFT4G:预测氨基酸取代是否会影响蛋白质功能

使用方法：
1.下载对应物种的SIFT数据库https://sift.bii.a-star.edu.sg/sift4g/

也可以自己制作SIFT数据库
2.下载SIFT4G(win,linux,mac均可，是一个Jre程序，依赖Java)
3.使用SIFT4G官方使用说明
 综合性突变危害预测软件
dbNSFP是一个综合性的有害突变预测软件库，包含多种算法的多个软件。详情

自己制作SIFT数据库

sift4g独立构建基因组，安装过程

root权限下的操作方法

1.检查是否安装perl(linux自带的有)

perl -Version

2.下载安装DBI

wget -C http://www.cpan.org/modules/by-module/DBD/DBI-1.643.tar.gz
tar -zxvf DBI-1.643.tar.gz
cd DBI-1.643
perl Makefile.PL
make
make test
make install

3.安装LWP

perl -MCPAN -e "install Bundle::LWP"
perl -MCPAN -e "install HTML::Formatter"
perl -MCPAN -e "install IO::Socket::SSL and OpenSSL"

4.安装bioperl

wget https://cpan.metacpan.org/authors/id/C/CJ/CJFIELDS/BioPerl-1.7.7.tar.gz
tar -zxvf BioPerl-1.7.7.tar.gz
cd BioPerl-1.7.7
perl Makefile.PL
make
make test
make install

5.安装Switch.pm

sudo apt-get install libswitch-perl

6.安装SIFT4G

git clone --recursive https://github.com/rvaser/sift4g.git sift4g
cd sift4g
make

6.1. sift4g依赖项安装
g++(4.9+)
GNU Make
nvcc(2.*+)(optional)
安装g++和make

apt-get install g++
apt-get install make

非root用户的操作方法

#conda安装最新版的gcc
conda install gcc_impl_linux-64 
#conda安装g++

1.下载源文件
[DBI]http://www.cpan.org/modules/by-module/DBD/DBI-1.643.tar.gz
[LWP]https://cpan.metacpan.org/authors/id/O/OA/OALDERS/libwww-perl-6.47.tar.gz
[BIoPerl]https://cpan.metacpan.org/authors/id/C/CJ/CJFIELDS/BioPerl-1.7.7.tar.gz
2.perl的包在make前要设置prefix

#当前工作目录为~/perl5
tar -zxvf DBI-1.643.tar.gz
cd DBI-1.643
perl Makefile.PL PREFIX=~/perl5
make
make test
make install
echo 'export PERL5LIB="$HOME/perl5/:$PERL5LIB" '>>~/.bashrc

其他包也采用上述方法安装，或者采用其他方法安装。

上述软件在安装过程中，bioperl的检查make test可能无法通过.但是不影响后续安装和使用。

开始构建陆地棉的数据库参考官方github

1.先准备相关的文件

genome基因组、gtf注释或者gff3注释、（protein和dbSNP可选文件）
NCBI或者ensemble都可以下载相关信息

2.创建文件夹，移动上述文件到对应的目录

cd test_files
cp homo_sapiens-test.txt Gossypium_hirsutum_config.txt
mkdir Gossypium_hirsutum
mkdir Gossypium_hirsutum/gene-annotation-src
mkdir Gossypium_hirsutum/chr-src
mkdir Gossypium_hirsutum/dbSNP
mv Gossypium_hirsutum.fa.gz Gossypium_hirsutum/chr-src
mv Gossypium_hirsutum.gtf.gz Gossypium_hirsutum/gene-annotation-src
mv Gossypium_hirsutum.protein.fa Gossypium_hirsutum/gene-annotation-src

特别提示，蛋白质文件一定要解压缩为fa或者fasta格式，如果是压缩文件直接运行，后续会报错。

`gunzip Gossypium_hirsutum.protein.fa.gz`

3.修改Gossypium_hirsutum_config.txt的配置文件里面的路径。

设置

PARENT_DIR=/mnt/e/sift4g/sift4g_Create/test_files/Gossypium_hirsutum
ORG=Gossypium_hirsutum_sapiens
ORG_VERSION=Gh.V1
SIFT4G_PATH=/mnt/e/sift4g/sift4g/bin/sift4g
PROTEIN_DB=/mnt/e/sift4g/sift4g_Create/test_files/Gossypium_hirsutum/gene-annotation-src/Gossypium_hirsutum.protein.fa

检查GENETIC_CODE_TABLE和MITO_GENETIC_CODE_TABLE是否正确设置。
可选设置：<DBSNP_VCF_FILE>

如果没有dbSNP文件，直接在文件中这一行前面加上`#`即可。

# DBSNP_VCF_FILE=Homo_sapiens.vcf.gz

上面修改配置文件时，里面的路径一定要写完整的绝对路径。不要使用相对路径，不然程序会找不到文件地址。

绝对路径格式/mnt/e/sift4g/genome.fa.gz

如果没有gtf文件，有gff文件，使用`gffread`工具转换即可。请一定要确认你最终的gtf文件的第9列包含gene_biotype信息，不然生成的数据库无法使用。

#检测gff3文件中第9列包含的信息
awk 'BEGIN{FS=OFS="\t"} $3=="gene"{split($9, a, ";"); for(i in a){split(a[i], b, "="); if(++c[b[1]]==1) print b[1]}}' Gossypium_hirsutum.gene.gff3

上述命令会输出你的gff3文件的第9列里包含的
gffread下载安装

gffread-0.12.2预编译的二进制文件
gffread-0.12.2源码
如果下载二进制，可以直接解压使用。
如果下载源码，本地编译，安装即可。
或者
github安装也可以

git clone https://github.com/gpertea/gffread
cd gffread
make release

三选一即可，安装完成gffread。

将gff转换为gtf

gunzip Gossypium_hirsutum.gff.gz
gffread Gossypium_hirsutum.gff -T -o Gossypium_hirsutum.gtf
gzip Gossypium_hirsutum.gtf

4. 生成本地sift4g数据库（推荐使用有root权限的Ubuntu）

perl make-SIFT-db-all.pl -config ./test_files/Gossypium_hirsutum_config.txt

构建过程可能需要1-24h,当显示All done!即代表构建完成。（个人电脑i9的跑了8天8夜）

经过分析发现，其中使用sift4g的这一步可以优化，sift4g支持多线程，默认是8个，可以手动增加进程数量到你的硬件可以承受的范围，以加快分析速度。

修改make-SIFT-db-all.pl的第109行内容中，“ -d ”前面添加上`-t 96`.96是你指定的进程数。109行修改后如下所示。本人亲测在24个cpu上运行96个进程，完全正常。

my $sift4g_command = $meta_hash{"SIFT4G_PATH"} .  " -t 96 -d " . $meta_hash{"PROTEIN_DB"} . " -q " . $meta_hash{"PARENT_DIR"} . "/all_prot.fasta --subst " .  $meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"SUBST_DIR"} . " --out " .  $met    a_hash{"PARENT_DIR"} . "/" . $meta_hash{"SIFT_SCORE_DIR"} . " --sub-results " ;

print $sift4g_command . "\n";

#`$siftsharp_command`;
`$sift4g_command`;

如果你研究的物种基因组非常大，可以手动运行第109行-114行命令这一步。手动运行完成后，把这段代码和109行之前的命令使用`#`注释掉，继续运行make-SIFT-db-all.pl就可以正常后续分析了。

手动运行109-114行的方法

#使用conda安装sift4g,相信我使用conda安装是最正确的选择，80%的错误是sift4g没有成功安装造成的，它需要特定的GCC,G++.5%的错误是第一步没有正确修改路径造成的。15%的错误是gtf文件格式不正确造成的。
#先使用conda创建新的虚拟环境，再进入虚拟环境
conda activate sift4g
conda install sift4g

外面运行sift4g的命令

conda activate sift4g
sift4g -t 96 -d /share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/gene-annotation-src/Gossypium_hirsutum.protein.fasta -q /share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/all_prot.fasta --subst /share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/subst --out /share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/sift4_predictions --sub-results

如果你的基因组非常大，可以分成两部分分别运算，再合成一块。
分别从all_prot.fasta提取前面的n条染色体和subst文件里对应的n条染色体对应的基因的.subst文件。
修改输出目录sift4_predictions为新的文件夹part1_predictions。
最后把两次运行的--out后的输出目录part1_predictions,part2_predictions里的alignments.txt，按照先后顺序合并成为一个alignments.txt。把part1_predictions,part2_predictions里所有的.fasta和.SIFTprediction，以及合并后的alignments.txt移到sift4_predictions里。注释make-SIFT-db-all.pl文件115行之前的命令全部注释，注意变量要不能注释。再次运行修改后的make-SIFT-db-all.pl即可。

5. 检查生成的数据库

cd Gossypium_hirsutum/Gh.V1
zcat chr1.gz|more  
zcat chr1.gz|grep CDS|more
more CHECK_GENES.LOG

上述只是查看了chr1,根据实际染色体数量，分别查看。sift数据应在第10-12列中，如果有很多是NA，则可能有问题。

注意：只有cds有预测的值，5’-UTR和3‘-UTR都是NA,这是正常的。

check_genes.log文件最后一行应该总结了整个基因组的预测。如果最后3个不同列的百分比很高，则数据库构建完成。
check_genes.log文件分为如下4列

Chr
Gene with SIFT scores
Pos with SIFT scores
Pos with Confident Scores

6.用数据库注释vcf文件

6.1下载SIFT4G_Annotator.jar下载地址

6.2 运行之前先对vcf文件进行排序（！VCF文件必须按染色体和位置排序才能正确注释。）

一般情况下，输出的vcf都是按照染色体的顺序排好的。

6.3 运行程序，进行注释

运行之前，先解压自己构建或者从sift4g下载的对应的基因组版本的数据库文件。如果sift4g数据库的物种参考基因组版本和你用的不一致。使用picard或GATK其他基因组版本转换工具把你的变异文件转换到数据库的基因组版本。转换需要一个chain文件。
或者按照上述文件，自行构建该版本的数据库。

java -Xmx128G -jar SIFT4G_Annotator.jar  -c -i < 输入vcf文件的路径 > -d < SIFT4G数据库目录的路径 > -r < 结果文件夹的路径 > -t

参数说明：
-t To extract annotations for multiple transcripts (Optional)可选参数

7.sift4g结果的可视化

基因组可视化工具很多
基于R语言的也有很多
CMplot或者ggplot或者ChromHeatMap，当然可以画圈图了。
可视化参考全基因组的reads覆盖
 CMplot包的使用

cat 94_chr.snp_SIFTannotations.xls||grep DELETERIOUS|grep -v Mt|awk '{print $6"\t"$1"\t"$2}' >94_chr.badSNP.tab

R语言里CMplot可视化

rm(list=ls())
setwd("~/cotton/")
install.packages("CMplot")
library(CMplot) #V3.6.2
mydata<-read.table("94_chr.badSNP.tab",header=FALSE)
names(mydata) <- c("snp","chr","pos")
head(mydata)
# snp         chr       pos
# snp1_1    1        2041
# snp1_2    1        2062
# snp1_3    1        2190
#可以自行修改window size默认的是1e6即1MB.
CMplot(mydata,plot.type="d",bin.size=1e6,col=c("darkgreen","yellow", "red"),file="jpg",memo="snp_density",dpi=300) 
##图片会直接输出到工作目录，注意保存重命名

报错信息请到官方github查看解决方案

注意查看`all_prot.fasta`，这个文件。

运行这行命令cat all_prot.fasta|grep ">" |sort|uniq -d，如果有行输出，说明你的all_prot.fasta里面有重复的行，就是显示的这些行。
如果需要重新运行，一定要删除all_prot.fasta这个文件。

报错信息处理

我第一次使用的是压缩的protein文件，结果报错。

/mnt/e/sift4g/sift4g/bin/sift4g -d /mnt/e/sift4g/sift4g_Create/test_files/Gossypium_hirsutum/gene-annotation-src/GossypiumHirsutum_ASM98774v1_protein.fa -q /mnt/e/sift4g/sift4g_Create/test_files/Gossypium_hirsutum/all_prot.fasta --subst /mnt/e/sift4g/sift4g_Create/test_files/Gossypium_hirsutum/subst --out /mnt/e/sift4g/sift4g_Create/test_files/Gossypium_hirsutum/SIFT_predictions --sub-results
** Checking query data and substitutions files **
* processing queries: 100.00/100.00% *

** Searching database for candidate sequences **
[ERROR:src/chain.c:101]: chain is empty after encoding, see scorerEncode function

上述报错之后，需要修改最开始运行的perl脚本之后，重新运行。从前面运行到此处，已经运行4*24h。
修改make-SIFT-db-all.pl这个文件后，继续运行脚本。

报错python文件无法找到，原因是本机没有配置pythonpath.
注意：这个脚本里，用到3个python脚本文件，建议修改本机的pythonpath为python2.
如果你默认的python是python3，也可以直接去修改调用python命令的地方，修改python为python2,至少要修改2个文件perl文件的python调用。
其中一处是：make-SIFT-db-all.pl的137行的python修改为python2,
make-sift-scores-db-batch.pl的65行的python修改为python2

报错原因是gcc版本低于4.9，要求是gcc 4.9+

/share/softwares/sift4g/sift4g/bin/sift4g: /lib64/libstdc++.so.6: version `CXXABI_1.3.8' not found (required by /share/softwares/sift4g/sift4g/bin/sift4g)
/share/softwares/sift4g/sift4g/bin/sift4g: /lib64/libstdc++.so.6: version `GLIBCXX_3.4.21' not found (required by /share/softwares/sift4g/sift4g/bin/sift4g)

解决办法，重新安装gcc8.3.0 从镜像下载 https://mirrors.aliyun.com/gnu/gcc/gcc-8.3.0/ ，并更新动态库。安装gcc8.3.0请严格按照此链接执行。安装时间可能需要几个小时，内存低于128g的慎重安装。
内存不足，还是安装低版本的吧。https://mirrors.aliyun.com/gnu/gcc/gcc-5.5.0/，安装上述链接，安装依赖，之后安装低版本gcc5.5.0

tar -zxvf gcc-8.3.0.tar.gz
#依赖太多,用conda,或者自己下载也行
conda install gcc
cd gcc-8.3.0
./configure --disable-multilib --enable-languages=c,c++ --prefix=$HOME/local
make -j5
make -j install

修改环境变量

export LD_LIBRARY_PATH=$HOME/local/lib64`
export PATH= /home/username/myinstall/gcc-4.8/bin: /home/username/myinstall/gcc-4.8/lib64$PATH  
export LD_LIBRARY_PATH= /home/username/myinstall/gcc-4.8/lib64: /home/username/myinstall/gcc-4.8/lib$LD_LIBRARY_PATH

报错信息汇总，在redhat上遇到的，程序一直执行不报错，但是并没有生成Gh.V1的数据库。

输出的信息如下：

start siftsharp, getting the alignments
/share/softwares/sift4g/sift4g/bin/sift4g -d /share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/gene-annotation-src/Gossypium_hirsutum.protein.fasta -q /share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/all_prot.fasta --subst /share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/subst --out /share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/SIFT_predictions --sub-results

这一步比较耗时，请耐心等待。在/share/softwares/sift4g/scripts_to_build_SIFT_db/test_files/Gossypium_hirsutum/SIFT_predictions目录里应该有预测的变异是可接受或有害的文件。如果没有，请检查你的gtf文件格式是否正确。

如果你生成

作者的脚本核心是：perl+python+shell三种脚本语言混合使用。

快速查询perl文件中包含`python`的行的方法

egrep python *.pl

关于GCC的安装，强力建议root安装，因为目前为止，各种方法都试过了。但是安装的并不成功。请参考此地址下载安装gcc 9.3.0

如果你的sift4g报错是

sift4g: /lib64/libstdc++.so.6: version `CXXABI_1.3.8' not found (required by sift4g)
sift4g: /lib64/libstdc++.so.6: version `GLIBCXX_3.4.21' not found (required by sift4g)

请使用find / -name "libstdc++.so.6"找出你系统中最新版本的路径，就是下面两行命令前面的路径，后面的路径是对应上面报错的路径。

gcc安装完成后，一定要做下面的操作

cp /share/soft/ohpc/pub/compiler/gcc/8.3.0/lib64/libstdc++.so.6 /lib64/
ln -s /share/soft/ohpc/pub/compiler/gcc/8.3.0/lib64/libstdc++.so.6.0.25 /lib64/libstdc++.so.6.0.25

最后编辑于：2020.11.19 13:55:45

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 204,189评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,577评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,857评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,703评论 1赞 276
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,705评论 5赞 366
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,620评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,995评论 3赞 396
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,656评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,898评论 1赞 298
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,639评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,720评论 1赞 330
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,395评论 4赞 319
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,982评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,953评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,195评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 44,907评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,472评论 2赞 342