美文网首页新的蛋白功能预测
SIFT4G:预测氨基酸取代是否会影响蛋白质功能

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

作者: wo_monic | 来源:发表于2020-03-13 22:26 被阅读0次

SIFT4G

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

自己制作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用户的操作方法

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工具转换即可。

gffread下载安装

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夜）

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，则可能有问题。
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 -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包的使用

报错信息请到官方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

报错信息汇总，在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目录里，如果有alignments.txt，查看下文件里包含所有的染色体数据。可以按照我下面的代码修改make-SIFT-db-all.pl，之后重新运行这个脚本即可正常生成数据库。

主要是修改第121行之前的命令全部加上#不执行，再重新运行修改后的脚本即可生成Gh.V1里的数据库。

#!/usr/bin/perl -w

# create the SIFT databases, must be run in scripts directory because it's lookin for metadocs

use strict;
#require 'common-utils.pl';
#use Getopt::Std;
use File::Basename;
use Cwd qw(abs_path);
my $directory_of_script = dirname(abs_path(__FILE__));
require $directory_of_script . '/common-utils.pl';

use Getopt::Long qw (GetOptions);

my %options  = ();
#getopts ("mhf:", \%options);

my $help = "";
my $ensembl_download = "";
my $metafile = "";
GetOptions ('help|h' => \$help,
            'ensembl_download' => \$ensembl_download,
            'config=s' => \$metafile );

if ($help || !$metafile || !$metafile) {
        print "make-SIFT-db-all.pl\n";
        print " -h  help\n";
        print " -config  [required config file]\n";
        print " -ensembl_download \n";
}

my $meta_href = readMeta($metafile);
my %meta_hash = %{$meta_href};

make_dir ($meta_hash{"PARENT_DIR"});
make_dir ($meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"ORG_VERSION"});

for my $key (keys %meta_hash) {
        my $value = $meta_hash{$key};
#       print $value . "\n";
        if ($key =~ /_DIR$/ && $key ne "PARENT_DIR") {
#               print "making $value\n";
                my $dir = $meta_hash{"PARENT_DIR"} . "/" . $value;
                if (! -d $dir) { mkdir($dir, 0775) or die $dir . $!; }
        }
}
my $siftscore_dir = $meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"SINGLE_REC_WITH_SIFTSCORE_DIR"};
if (! -d $siftscore_dir) { mkdir ($siftscore_dir, 0775); }

if ($ensembl_download) {
        # not manual, this is from Ensembl and needs to download
        print "downloading gene annotation\n";
        #`perl download-annotation-srcs.pl $metafile`;
        print "done downloading gene annotation\n";

        print "downloading fasta files\n";
        #`perl download-fasta.pl $metafile`;
        print "done downloading DNA fasta sequences";

        print "download dbSNP files\n";
        #`perl download-dbSNP-files.pl $metafile`;
        if (exists ($meta_hash{"DBSNP_VCF_FILE"}) && $meta_hash{"DBSNP_VCF_FILE"} ne "")
        {
                print "splitting dbSNP file\n";
                #`perl split-dbSNP-by-chr.pl $metafile`;
        }
} # end if downloading files from Ensembl

print "converting gene format to use-able input\n";
#`perl gff_gene_format_to_ucsc.pl $metafile`;
#`perl ensembl_gene_format_to_ucsc.pl $metafile`;
print "done converting gene format\n";

# decompress chromosome files so that Bio::DB can process them
my $fasta_dir =  $meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"CHR_DOWNLOAD_DEST"} ;
=pod
if (glob ("$fasta_dir/*.gz")) {
        $fasta_dir .= "/*.gz";
        `gunzip $fasta_dir`;
        # check that all DNA files finished downloading, otherwise do it again.
        # because all future programs won't work properly otherwise.
        if ($?) {
                die "DNA files do not exist or did not unzip properly\n";
        }
} elsif (glob ("$fasta_dir/*.fa") || glob ("$fasta_dir/*.fasta")) {
        # *.fa or *.fasta files already exist
} else {
        die "no DNA fasta files in $fasta_dir\n";
}
=cut
print "making single records file\n";
#`perl make-single-records-BIOPERL.pl $metafile`;
print "done making single records template\n";

print "making noncoding records file\n";
#`perl make-single-records-noncoding.pl $metafile`;
print "done making noncoding records\n";

print "make the fasta sequences\n";
#`perl generate-fasta-subst-files-BIOPERL.pl $metafile`;
print "done making the fasta sequences\n";

print "start siftsharp, getting the alignments\n";

my $combine_prot_fasta_command = "for file in " .  $meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"FASTA_DIR"} . "/*.fasta ; do cat \"\$file\" >> " . $meta_hash{"PARENT_DIR"} . "/all_prot.fasta; done";

#`$combine_prot_fasta_command`;


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

print $sift4g_command . "\n";

#`$sift4g_command`;

if ($? != 0) {
        exit (-1);
}

print "done getting all the scores\n";

print "开始后续分析了，2020.11.14\n";
print "populating databases\n";
# check quality
# Pauline to do , by chromosame in script

`perl make-sift-scores-db-batch.pl $metafile`;

print "checking the databases\n";
system ("perl check_genes.pl $metafile");

my @chromosomes =  getChr ($meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"GENE_DOWNLOAD_DEST"});

foreach my $chr (@chromosomes) {

        my $db_file = $meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"SINGLE_REC_WITH_SIFTSCORE_DIR"} . "/" . $chr .  "_scores.Srecords.with_dbSNPid.sorted";
        my $check_outfile = $meta_hash{"PARENT_DIR"}  .  "/" . $meta_hash{"ORG_VERSION"} . "/" . $chr .  "_SIFTDB_stats.txt";
        `python2 check_SIFTDB.py $db_file $check_outfile`;
}

my $pep_check_file = $meta_hash{"PARENT_DIR"} . "/ENS_peptide_check.txt";
`perl checkENSPep.pl $metafile ENS $pep_check_file`;

my $final_outfolder = $meta_hash{"PARENT_DIR"} ."/". $meta_hash{"ORG_VERSION"};
system ("cp  $metafile $final_outfolder");

# make some space
my $zip_dir = $meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"CHR_DOWNLOAD_DEST"} . "/*" ;
print "zipping up $zip_dir\n";
system ("gzip $zip_dir");
my $rm_dir = $meta_hash{"PARENT_DIR"} . "/" . $meta_hash{"SINGLE_REC_WITH_SIFTSCORE_DIR"} . "/*";
system ("rm $rm_dir");

print "All done!\n";
exit (0);

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

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

egrep python *.pl

相关文章

网友评论

    本文标题：SIFT4G:预测氨基酸取代是否会影响蛋白质功能

    本文链接：https://www.haomeiwen.com/subject/uljqrhtx.html

    延伸阅读

    深度阅读

    • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

    栏目导航

    热点阅读

    新的蛋白功能预测

    关于我们|服务条款|联系我们|SIFT4G:预测氨基酸取代是否会影响蛋白质功能|投稿指南|网站地图|RSS订阅|排版工具|手机版

    提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

    Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

    备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

    本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

    所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！