从单个PDB文件中提取多个蛋白质链

Question

我有一个PDB文件，其中包含多个链，但没有chainid。我想用R来指定chainid，这样我就可以分析单个蛋白质链并找到每个蛋白质链中的特定位点。

我正在使用Rpdb从下面提取文件和示例数据（来自单个pdb文件的每个链的顶部几行）。

REMARK 99 Chain ID : 1 
REMARK 99 Residues : 593 
REMARK 99 Atoms : 4782 
REMARK 99 File  : final.sc.pdb 
ATOM  1 N MET  1  17.471 -55.657 42.605 1.00 0.00    
ATOM  2 CA MET  1  17.516 -55.479 41.136 1.00 0.00    
ATOM  3 CB MET  1  16.328 -56.188 40.460 1.00 0.00    
ATOM  4 C MET  1  17.525 -54.045 40.745 1.00 0.00    
ATOM  5 O MET  1  17.991 -53.186 41.492 1.00 0.00    
ATOM  6 CG MET  1  14.961 -55.764 41.001 1.00 0.00   C 
ATOM  7 SD MET  1  14.550 -56.460 42.632 1.00 0.00   S 
ATOM  8 CE MET  1  12.951 -55.613 42.782 1.00 0.00   C 
ATOM  9 N THR  2  17.012 -53.760 39.535 1.00 0.00    
ATOM  10 CA THR  2  16.993 -52.420 39.040 1.00 0.00    
ATOM  11 CB THR  2  16.552 -52.347 37.612 1.00 0.00       
TER 
REMARK 99 Chain ID : 1 
REMARK 99 Residues : 531 
REMARK 99 Atoms : 4211 
REMARK 99 File  : final.sc.pdb 
ATOM  1 N MET  1  55.179 17.162 2.445 1.00 0.00    
ATOM  2 CA MET  1  55.489 16.069 3.613 1.00 0.00    
ATOM  3 CB MET  1  55.199 16.623 5.019 1.00 0.00    
ATOM  4 C MET  1  53.890 15.434 3.310 1.00 0.00    
ATOM  5 O MET  1  52.902 15.782 3.971 1.00 0.00    
ATOM  6 CG MET  1  56.062 17.833 5.341 1.00 0.00   C 
ATOM  7 SD MET  1  55.937 18.517 7.006 1.00 0.00   S 
ATOM  8 CE MET  1  56.886 17.217 7.874 1.00 0.00   C 
ATOM  9 N ALA  2  53.854 14.445 2.424 1.00 0.00    
ATOM  10 CA ALA  2  52.895 13.660 2.231 1.00 0.00    
ATOM  11 CB ALA  2  53.134 12.918 0.924 1.00 0.00    
ATOM  12 C ALA  2  52.253 12.986 3.391 1.00 0.00    
ATOM  13 O ALA  2  51.034 12.834 3.347 1.00 0.00 
TER 

列名进行RPDB添加（注：chainid，插入和segid没有值）：

recname eleid elename alt resname chainid resid insert  x1  x2  x3 occ temp segid 

有谁知道一种方式，说chainid的补充？ 谢谢！

Answer 1

通过使用“TER”明确了蛋白质链的开头和结尾，我能够做什么工作了，但如果任何人有一个更好/更顺畅/更快的方式请让我知道：

#works for pdb file with two chains 
pdb.input.table=read.delim(file.choose(),sep="",header=F) 

#pdb chain splitting 
chainAstart=1 
chainAend=which(pdb.input.table=="TER")[1] 
chainBstart=which(pdb.input.table=="TER")[1]+1 
chainBend=which(pdb.input.table=="TER")[2] 

new.chain.id=c(rep("A",chainAend),rep("B",chainBend-chainAend)) 

pdb.dock.input=cbind(pdb.input.table,new.chain.id)