使用方法: 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文件,直接在文件中这一行前面加上#
即可。
# 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 ```