参考4: Generalized Fragment Picking in Rosetta: Design, Protocols and Applications
教程位于: $ROSETTA/demos/public/fragment_picking/BestFragmentsProtocol
高能提示:本文需要下载大量的数据库数据，因此请准备足够的硬盘空间(>=100GB)
前言
3~9个氨基酸的短肽片段(Fragments)是重头算预测蛋白和Loop结构预测的重要元件，本文将以简要概述如使用快速地从pdb数据库中获取有用的Fragments。
本文的内容分别如下介绍:
- Fragment_picker的计算原理;
- Fragment_picker的使用方法;
- make_fragments.pl自动生成输入参数;
一、Fragment_Picker计算原理:
Fragment_Picker是Rosetta中最新的短肽片段选择应用，相对于老版本NNMAKE，从精度上来说，基本持平，但是添加了新的模块与功能:
-
兼容了实验数据，如NMR chemical shift data;
-
完全可自定义的函数系统;
-
完全可根据已知信息，有选择性地选取Fragment。
总的来说FragmentPicker的计算分为3步:
1 创建打分函数系统，读取整个Vall数据库；(处理)
2 从数据库中生成短肽片段，打分，并收集打分较优者；(打分-选取)
3 最后在预收集的库中，生成最终短肽片段库并输出。(输出)
以下将简要的介绍数据库处理、打分和选取的细节:
1 Vall数据库的读取与片段选取
Vall数据库由许许多多的PDB结构组成，而每个PDB会被进一步被"剪切"成片段库(vall chunk)，其中vall chunk记录着每个Fragment的序列和结构信息(vall residue)。
vall数据库.jpgVall数据库的搜索与处理:
FragmentPicker在读取完整个数据库后，将逐步对每个VallChunk进行匹配，一但VallChunk满足需求(通过了VallChunkFilter过滤器的标准)，VallChunk将会被分解成连续的短肽片段库。其中的每条短肽将会被打分，评估它在我们提交序列中的匹配度，得分较高的这些短肽片段将被收集器(fragment collector)收集起来。
目前ROSETTA安装包自带的几个vall数据库如下:
- $ROSETTA/tools/fragment_tools/vall.jul19.2011.gz (推荐)
- $ROSETTA/tools/fragment_tools/vall.apr24.2008.extended.gz
- $ROSETTA/main/database/sampling/filtered.vall.dat.2006-05-05.gz
- $ROSETTA/main/database/sampling/small.vall.gz(用于测试)
注: 生成自己vall库的方法: 使用Rosetta自带的pdb2val工具，pdb2val工具位于: $ROSETTA/tools/fragment_tools/
2 短肽片段的打分
短肽片段的打分本质上是短肽片段中每个氨基酸得分(Ns)的线性加和，其中Si和Wi由wghts权重文件设置。
打分值代表选取的片段与输入序列的相似度，打分的区间为0~1, 0代表完美匹配，1代表完全不相关。其中打分项包括了RamaScore、SecondarySimilarity、SequenceIdentity、SecondaryIdentity、Constraint、ProfileScore。
下图为打分项与相关需要输入文件:
打分项.jpg- 其中.ss文件包括了3项主要的打分:SequenceSimilarity/SecondarySimilarity/RamaScore.
- 输入的.fasta文件主要用于SequenceSimilarity
- 输入的限制文件: 二面角/距离限制用于AtomPairConstraintScore和DihedralConstraintScore的计算
- PDB文件，仅仅用于对照作用。
注意: 在计算打分过程中，如果打分项的某一项不满足我们所设定的阈值，那么它就立马被淘汰，而不会去计算剩余的打分。
3 配额机制(Quota)
打分系统总是偏向于分值高的片段，因此片段库的二级结构可能会过度倾向于某种局部构象，因此打分机制中引入了"配额"机制，用于提升所选片段的多样性，并且有利于避免局部二级结构预测错误的问题。
配额机制的原理:
- 使用三个不同的二级结构预测器，对片段的二级结构属性进行预测。并生成3个不同的短肽片段库；
- 根据二级结构的预测的可能性，分别选取相应比例数量的短肽(打分最佳者)。
- 从3个软件的H(螺旋)、L(Loop)、E(β-折叠)库中，合成一个短肽库。
Quota目前推荐利用PsiPred, Porter和SAM作为二级结构预测器。
Quota实例.jpg此例中: PsiPred, Jufo(目前已经改用Porter作为预测器)和SAM的权重分别是0.6/0.2/0.2，代表最后生成的top200片段库60%由PsiPred提供，而剩下的40%由Jufo和SAM提供。
并且每个片段库中的短肽数量，也会根据预测二级结构的预测类型的占比进行分配，比如PsiPred预测39号氨基酸D有68.6%的概率为Coil结构，那么在最终的短肽库中，他的占比就是: 200*0.6*0.686~=31条，那么这31条序列都是在Coil二级结构分类下，打分最高的短肽。这就是Quota机制的运算机理。
二、Rosetta Fragment_picker的使用
Fragment Picker有三种基本的运行模式:
- Quota: 使用配额模式，生成de novo蛋白预测的短肽片段库；
- Best fragments: 根据限制或则已知信息来寻找最佳匹配的短肽片段库；
-
denovo design: 当蛋白序列发生变化时，不断地重新进行实时Fragment获取，被
flexibleLoopDesign.linuxgccrelease应用调用。(此处不详细深究)
1 Best fragments模式的使用
此模式下，会有偏向性选择与打分最佳的短肽片段，可能会导致二级结构采样有偏向性，只有在具有实验或则其他已知信息是才使用。
1.1 获取二级结构预测文件.ss
准备FASTA文件:2JSV.fasta, 用于二级结构的预测:
> 2JSV X
MQYKLILNGKTLKGETTTEAVDAATAEKVFKQYANDNGVDGEWTYDDATKTFTVTE
使用本地psipred对序列进行二级结构预测:(安装见本文目录，比较麻烦)
runpsipred 2JSV.fasta
下载解压结果文件后，其中以*.ss2的文件是我们所需要的，并将其重命名为:2jsv.psipred.ss2
1.2 创建wghts权重文件
.wghts文件是用于定义打分系统，来评估每一个我们搜索的Fragment的分值，程序会将Top200打分的Fragment用作输出。
一个权重文件至少应该有4列数据，分别是score name、priority、wght和max_allowed。
max_allowed代表最大值，如果一条短肽的打分超出该分值后，直接被忽略。
新建simple.wghts文件并填写以下内容:
# score name priority wght max_allowed extras
RamaScore 400 2.0 - predA
SecondarySimilarity 350 1.0 - predA
FragmentCrmsd 0 0.0 -
predA代表: 链接的ss文件，如下述-frags::ss_pred 2jsvX.psipred.ss2 predA,可以是任意的字段。
1.3 准备运行参数文件
新建参数文件flags, 并填写以下内容:
# Input databases
-in::file::vall $ROSETTA/tools/fragment_tools/vall.jul19.2011.gz # 最新的vall
# 初始输入文件
#-in::file::checkpoint input_files/2jsvX.checkpoint
-in::file::fasta 2jsvX.fasta
-in::file::s 2jsvX.pdb # 仅用于验证和测试Fragment选取的质量。
-frags::ss_pred 2jsvX.psipred.ss2 predA
# 设置权重文件
-frags::scoring::config simple.wghts
# What should we do?
-frags::bounded_protocol
# Fragment选择设置
-frags::frag_sizes 3 9
-frags::n_candidates 200
-frags::n_frags 200
# 输出设置
-out::file::frag_prefix frags
-frags::describe_fragments frags.fsc
1.4 运行fragment_picker
fragment_picker.mpi.macosclangrelease @flags
此时将输出文件:frags.200.3mers和frags.200.9mers用于输入；
此外fsc files记录着fragment的信息。
2 Quota模式的使用
增加短肽二级结构多样性的方法，多用于de novo预测。
2.1 获取二级结构预测文件.ss
需要分别从3个预测器中进行预测，获取.ss格式的文件，作为输入文件并使用后文的make_fragments.pl脚本来进行Fragment选取。
2.2 定义配额与weights文件
新建Quota.def file文件，默认使用psipred配额为0.6，其余两个分别为0.2。
#pool_id pool_name fraction
1 psipred 0.6
2 porter 0.2
3 sam 0.2
新建文件quota-protocol.wghts，使用以下默认设置即可。
# score name priority wght min_allowed extras
SecondarySimilarity 350 0.5 - psipred
SecondarySimilarity 300 0.5 - sam
SecondarySimilarity 250 0.5 - jufo
RamaScore 150 1.0 - psipred
RamaScore 150 1.0 - sam
RamaScore 150 1.0 - jufo
ProfileScoreL1 200 1.0 -
PhiPsiSquareWell 100 0.0 -
FragmentCrmsd 30 0.0 -
2.3 准备运行参数文件
新建参数文件flags, 并填写以下内容:
# Input databases
-in::file::vall $ROSETTA/tools/fragment_tools/vall.jul19.2011.gz # 最新的vall
# Query-related input files
#-in::file::checkpoint fragment_picking_with_quota/input_files/2jsvX.checkpoint
-in::file::s fragment_picking_with_quota/input_files/2jsvX.pdb
-frags::ss_pred 2jsvX.psipred.ss2 psipred 2jsvX.sam.ss2 sam 2jsvX.jufo.ss2 jufo
# the name root for the output fragment files
-out::file::frag_prefix frags
# Show score components for each selected fragment
-frags::describe_fragments frags.fsc
# Weights file
-frags::scoring::config quota-protocol.wghts
# 短肽长度设置
-frags::frag_sizes 9 3
# 从1000个候选物中挑选出200条序列
-frags::n_candidates 1000
-frags::n_frags 200
# Quota.def file defines the shares between difefrent quota pools. The total should be 1.0
-frags::picking::quota_config_file quota.def
# Get rid of homologues fragments
#-frags::denied_pdb 2jsvX.homolog_vall
2.4 运行
fragment_picker.mpi.macosclangrelease @flags
3 控制Fragment Picker的行为
用户可以通过使用限制参数来改变和调整Fragment生成。
3.1 使用限制文件
只有两类限制是支持的: AtomPairConstraint以及DihedralConstraints。一般这些数据来源于NMR、交联质谱等，常见类型为原子间距离的限制。
使用限制文件:
-
需要在flag参数文件中添加
-constraints::cst_file $file选项即可。$file代表输入的限制文件全名。 -
还需要在weight文件中添加
AtomPairConstraintsScore 500 2.0 -
注意: 限制定义的两个原子不能超出滑动窗口的范围。
3.2 使用Torsion bin限制
序列中每个氨基酸的二面角分布数据可以来源于:机器学习、同源蛋白的晶体结构(比较靠谱)。
使用Torsion class score:
-
需要在flag文件中额外添加
-in::file::torsion_bin_probs $ABEGO选项，$ABEGO代表输入包含每个氨基酸在ABEGO分布信息的文件。(本质上就是预先定义了每个氨基酸二面角的区间，只允许Fragment在这些区域出现)。 -
还需要在weight文件中，添加
TorsionBinSimilarity 500 1.0 0.1
$ABEGO格式:
# HEADER
.........
11 T 1.0000 1.0000 1.0000 1.0000 1.0000
12 L 1.0000 1.0000 1.0000 1.0000 1.0000
13 K 1.0000 1.0000 1.0000 1.0000 1.0000
14 G 1.0000 1.0000 1.0000 1.0000 1.0000
15 E 1.0000 1.0000 1.0000 1.0000 1.0000
16 T 1.0000 0.0000 0.0000 0.0000 0.0000
17 T 1.0000 1.0000 1.0000 1.0000 1.0000
18 T 1.0000 1.0000 1.0000 1.0000 1.0000
19 E 1.0000 1.0000 1.0000 1.0000 1.0000
20 A 1.0000 1.0000 1.0000 1.0000 1.0000
.........
注: 此处5个浮点数，分别代表ABEGO的概率。如ID:16 此处修改了它BEGO的4个概率。因此对于此氨基酸，强制骨架二面角为'A'区域的短肽被保留。
此处的概率代表惩罚权重。如，概率设置为1.0时，惩罚打分为0；概率设置为0.0时，惩罚打分项最高，由TorsionBinSimilarity 500 1.0 0.1决定。
三、使用make_fragments.pl + PSIPRED3 快速获取Fragments
仅限Linux系统使用，除非做针对性修改。自动化程度高。本质上也是还是调用fragment_picker应用。
注意: 需要Python2环境。
1 运行前先修复错误:
打开:$ROSETTA/tools/fragment_tools/install_dependencies.pl文件,修复即可
# 将269行文字进行替换:
原: system("wget -N
修: system("wget -N ftp://ftp.ncbi.nlm.nih.gov/blast/temp/update_blastdb.pl");
sudo yum install -y perl-Switch
2 安装第三方依赖包:
如果图省事直接把一下命令的skip_nr删掉即可，也可以替换为所需要数据库的字段，如uniref90、uniref50，默认为下载NCBI的nr数据库。
注意1: 如果选择了uniref数据库，那么会覆盖当前已经下载的NCBI nr库。并且数据库的名称也会改成nr库。
$ROSETTA/tools/fragment_tools/install_dependencies.pl standard uniref90
# 安装完毕之后的软件路径如下，无需额外设置环境:
blast: /usr/local/rosetta_src_2019.12.60667_bundle/tools/fragment_tools/blast
psipred: /usr/local/rosetta_src_2019.12.60667_bundle/tools/fragment_tools/psipred
csblast: /usr/local/rosetta_src_2019.12.60667_bundle/tools/fragment_tools/csblast
sparks-x: /usr/local/rosetta_src_2019.12.60667_bundle/tools/fragment_tools/sparks-x
注: 如果想更新到最新版本，参考附录的更新方法。
3 修改make_fragments.pl文件:
# 33 行，根据编译的环境修改实际的Fragment_picker app的名称
my $FRAGMENT_PICKER = $ENV{"FRAGMENT_PICKER"} || "$Bin/../../main/source/bin/fragment_picker.mpi.macosclangrelease";
4 运行脚本计算:
make_fragments.pl使用PSI-BLAST + PSIPRED3来生成短肽片段。
# cmd1 best模式:
$ROSETTA/tools/fragment_tools/make_fragments.pl -verbose -use_vall_files $vall_path -frag_sizes 3 9 -n_frags 200 -n_candidates 1000 -nohoms <fasta file>
# cmd2 quota模式:
$ROSETTA/tools/fragment_tools/make_fragments.pl -verbose -porterfile *.porter5.ss2 -nohoms <fasta file>
- 默认使用:vall.jul19.2011数据库，也可以通过
-use_vall_files控制读入的数据库； - -nofrags 启用后，生成所有需要的文件，但是不运行Fragment picker app，以便以我们对Fragment Picker的行为进行控制，如添加限制，引入更多的实验信息等；
- -sam -porter 启用后等于使用Quota模式(不推荐!!，因为已经很难配置porter和sam了)，
- -psipredfile -samfile -porterfile 来读取对应的.ss文件(推荐)；
- -frag_sizes <int> <int> <int>… , 代表生成这些长度的短肽片段(默认3，9)；
- -n_frags 200，控制最后输入短肽片段库的大小（默认为200）。
- -n_candidates 1000，Quota模式使用，避免数据库过小而产生的错误(默认1000)；
- <fasta file>, 我们所提交的序列文件；
- -nohoms 排除同源序列，仅用于测试时使用。
注: 此处脚本的weight打分文件的参数与短肽片段的长度相关。还会添加SPARKS-X的打分信息，比自己来设置可能会更加靠谱。
四、升级最新二级结构预测器
如果想使用Quota模式进行Fragment的选取, 就需要本地配置安装多个蛋白二级结构预测器，目前最新测试结果表明: PSIPRED4.02(Q3~=0.84)、PORTER5(Q3=0.83)的预测准确度较高，因此本文还提供了安装和使用的方法。此外: Jpred4、SPIDER3等也是开源软件，有兴趣的朋友可以尝试安装。本文仅介绍psipred.4.02与Porter5的安装和配置方法。
跑分.png1 升级至psipred.4.02
cd $ROSETTA/tools/fragment_tools/
# clean
rm -rf psipred
# 下载
wget http://bioinfadmin.cs.ucl.ac.uk/downloads/psipred/psipred.4.02.tar.gz
tar -zxvf psipred.4.02.tar.gz && cd ./psipred/src
# 安装编译
make && make install
2 安装PORTER5:(Python3环境)
依赖: python3、numpy的配置不详细阐述。
- python3
- numpy
- HHblits
- UniRef20
- PSI-BLAST+(可选，精度提升)
- UniRef90(可选，精度提升)
# 配置UniRef20, 自行调整数据库地址:
mkdir /mnt/sdd/blastdb/UniRef20 && cd /mnt/sdd/blastdb/UniRef20
nohup wget -c -t 0 > download.log &!
# 解压、重命名得到.index和.ffdata文件
nohup tar zxvf uniprot20_2016_02.tgz &!
# HHblits源码编译安装:
git clone
mkdir -p hh-suite/build && cd hh-suite/build
cmake -DCMAKE_INSTALL_PREFIX=. ..
make -j 4 && make install
# PORTER5下载:
git clone
# 设置初始化Porter5环境: 创建配置文件
vi /mnt/sdd/software/Porter5/scripts/config.ini
[DEFAULT]
psiblast = /mnt/sdd/software/ncbi-blast-2.9.0+/bin/psiblast
uniref90 = $ROSETTA/tools/fragment_tools/databases/nr
hhblits = hhblits
uniprot20 = /mnt/sdd/blastdb/UniRef20/uniprot20_2016_02/uniprot20_2016_02
# 设置环境变量:
echo "export PATH=$PATH:/mnt/sdd/software/hh-suite/build/bin" >> ~/.zshrc
echo "export PATH=$PATH:/mnt/sdd/software/hh-suite/build/scripts" >> ~/.zshrc
echo "export PATH=$PATH:/mnt/sdd/software/Porter5" >> ~/.zshrc
echo "export Porter5=/mnt/sdd/software/Porter5" >> ~/.zshrc
运行方法:
# HHblits only
python3 $Porter5/Porter5.py -i *.fasta --cpu 8 --fast
# HHblits only + PSI-BLAST, 准确度更高。
python3 $Porter5/Porter5.py -i *.fasta --cpu 8
最后输出.fasta.ss3文件，这里已经是.ss格式，但是不可以直接被make_fragment.pl调用，我们需要进一步转化它的格式:
此处使用笔者的脚本, 将脚本保存到:$ROSETTA/tools/fragment_tools文件夹下
#!/usr/bin/env python
#-*- coding: utf-8 -*
'''
# 转化为3位小数;
# 顺序为 L H E。
porter5 输出的顺序是 H E L
'''
import os,sys
with open(sys.argv[1],'r') as f, open(str(sys.argv[1])[:-4]+'.porter5.ss2','w') as f1:
f1.write('# PSIPRED VFORMAT (PSIPRED V2.6 by David Jones)\n\n')
lines = f.readlines()[1:]
for line in lines:
line = (line.strip('\n')).split('\t')
index = line[0]
aa = line[1]
ss_type = line[2]
H = float(line[3])
E = float(line[4])
L = float(line[5])
newline = index.rjust(4) + ' ' + aa + ' ' + ss_type + ' ' + "%.3f" % L + ' ' + "%.3f" % H + ' ' + "%.3f" % E + '\n'
f1.write(newline)
输入命令:
python $ROSETTA/tools/fragment_tools/convert_ss.py *.fasta.ss3
即可得到对应的二级结构格式.ss2
如果需要调整quota模式的权重, 修改make_fragment.pl, 此处可见psipred:sam:porter = 3:1:1。由于proter5的精度也不低了。我认为可以调整为1:1或则2:1
#1842行左右, sub get_ss_quota_vote函数
sub get_ss_quota_vote {
my $ss_pred_name = shift;
if ( $ss_pred_name eq 'psipred' ) {
return 1;
}
elsif ( $ss_pred_name eq 'sam' ) {
return 1;
}
elsif ( $ss_pred_name eq 'porter' ) {
return 1;
}
else {
die "Error: don't recognize ss_pred_name $ss_pred_name!\n";
}
}
五、升级与迁移备份
如果需要更新Rosetta之前，可以先将数据备份，避免大量的下载和格式化耗时;
# 备份数据库
mv $ROSETTA/tools/fragment_tools/databases /somewhere
# 恢复数据库
mv /somewhere $ROSETTA/tools/fragment_tools/database
六、使用ROBETTA在线服务
如果觉得自行配置Fragment Picker比较难，又没有具体的距离限制等实验数据，可以考虑直接使用服务器:(排队时间根据具体情况而定)
默认为Quota模式(6:2:2 / Top200 selected): PSIPRED3 + Jufo + SAM, 数据库: vall.jul19.2011.gz。