2016-11-22 98 views
0

我试图从PDB输入文件中提取NAD配体ID的chainwise。我想保存输出文件如下:如果输入文件是1AHI.pbd,包含四条链A,B,C和d,输出应该是独立的文件根据每行的内容将PDB文本文件拆分为多个文件

1AHI_A.txt 
1AHI_B.txt 
1AHI_C.txt 
1AHI_D.txt 

低于我的脚本不能给出预期的输出。可能是脚本中的一些逻辑问题。我也收到任何错误。

from glob import glob 

in_loc = 'C:/Users/Documents/NAD/NAD/result/test_result_file/' 
out_loc = 'C:/Users/Documents/NAD/NAD/result/test_result_file/output/' 

def test(): 
    fnames = glob(in_loc+'*.pdb') 

    for each in fnames: 
    # This is the new generated file out of input file (.txt). 
     formatted_file = each.replace('pdb', 'txt') 

     formatted_file = formatted_file.replace(in_loc, out_loc) 

    # This is the input file 
     in_f = open(each, 'r') 

    # A new file to be opened. 
     out_f = open(formatted_file, "w") 

    # Filtering results from input file 
     try: 
      out_chain_list = filter_file(in_f) 
      for each_line in out_chain_list: 
       out_f.write(each_line) 

     # Closing all the opened files. 
      out_f.close() 
      in_f.close() 

     except Exception as e: 
      print('Exception for file: ', each, '\n', e) 
      out_f.close() 
      in_f.close() 


def filter_file(in_f): 
    ligand_id = ['NAD'] 
    chain_ids = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z'] 
    previous_chain_id = None 
    chain_list = [] 
    out_chain_list = [] 

    for each in fnames: 
     for line in map(str.rstrip, fnames): 
      if line[:6] != "HETATM": 
       continue 
      chainID = line[21:22] 
      ligandID = line[11:14].strip() 

      if ligandID in ligandID and chain_id in chain_ids: 

       if chain_id != previous_chain_id: 
        c_ls = [] 

       if c_ls: 
        out_file_name = in_f.name.replace(in_loc, out_loc) 
        out_file_name = out_file_name.replace('.pdb', '_'+previous_chain_id+'.txt') 
        out_file = open(out_file_name, "w") 
        for l in c_ls: 
        out_file.write(l) 
        out_file.close() 

       chain_list.append(line) 
       previous_chain_id = chain_id 

       out_chain_list += c_ls 

    return out_chain_list 

test() 

实施例:

输入文件:像

HETATM15202 PA NAD A 501  44.008 102.331 5.491 1.00 11.48   P 
HETATM15203 O1A NAD A 501  43.295 103.140 6.507 1.00 11.48   O 
HETATM15204 O2A NAD A 501  42.939 101.407 4.919 1.00 11.48   O 
HETATM15205 O5B NAD A 501  45.052 101.397 6.166 1.00 11.48   O 
HETATM15247 PA NAD B 501  36.790 111.512 38.592 1.00 11.25   P 
HETATM15248 O1A NAD B 501  37.248 110.563 37.565 1.00 11.25   O 
HETATM15249 O2A NAD B 501  35.692 110.795 39.337 1.00 11.25   O 
HETATM15250 O5B NAD B 501  36.174 112.802 37.915 1.00 11.25   O 
HETATM15292 PA NAD C 501  100.016 130.669 21.776 1.00 12.28   P 
HETATM15293 O1A NAD C 501  99.311 131.864 22.293 1.00 12.28   O 
HETATM15294 O2A NAD C 501  101.501 131.009 21.932 1.00 12.28   O 
HETATM15295 O5B NAD C 501  99.727 130.510 20.238 1.00 12.28   O 
HETATM15337 PA NAD D 501  78.237 158.792 22.383 1.00 11.99   P 
HETATM15338 O1A NAD D 501  79.297 157.907 21.808 1.00 11.99   O 
HETATM15339 O2A NAD D 501  78.807 160.217 22.362 1.00 11.99   O 
HETATM15340 O5B NAD D 501  78.069 158.416 23.905 1.00 11.99   O 

预期输出:

1AHI_A.txt:A链输出端(OUT放文件)

HETATM15202 PA NAD A 501  44.008 102.331 5.491 1.00 11.48   P 
HETATM15203 O1A NAD A 501  43.295 103.140 6.507 1.00 11.48   O 
HETATM15204 O2A NAD A 501  42.939 101.407 4.919 1.00 11.48   O 
HETATM15205 O5B NAD A 501  45.052 101.397 6.166 1.00 11.48   O 

1AHI_B.txt:B链输出(输出文件)

HETATM15247 PA NAD B 501  36.790 111.512 38.592 1.00 11.25   P 
HETATM15248 O1A NAD B 501  37.248 110.563 37.565 1.00 11.25   O 
HETATM15249 O2A NAD B 501  35.692 110.795 39.337 1.00 11.25   O 
HETATM15250 O5B NAD B 501  36.174 112.802 37.915 1.00 11.25   O 

1AHI_C.txt:C链输出(输出文件)

HETATM15292 PA NAD C 501  100.016 130.669 21.776 1.00 12.28   P 
HETATM15293 O1A NAD C 501  99.311 131.864 22.293 1.00 12.28   O 
HETATM15294 O2A NAD C 501  101.501 131.009 21.932 1.00 12.28   O 
HETATM15295 O5B NAD C 501  99.727 130.510 20.238 1.00 12.28   O 

1AHI_D.txt: D链输出(输出文件)

HETATM15337 PA NAD D 501  78.237 158.792 22.383 1.00 11.99   P 
HETATM15338 O1A NAD D 501  79.297 157.907 21.808 1.00 11.99   O 
HETATM15339 O2A NAD D 501  78.807 160.217 22.362 1.00 11.99   O 
HETATM15340 O5B NAD D 501  78.069 158.416 23.905 1.00 11.99   O 

我希望th是会帮助你理解的。

我想提取NAD(在输入文件的第三列)并以链式方式保存输出文件单独文件。

+0

这里的大多数人不会熟悉PDB文件和配体 - 可以你展示了一个数据的例子来说明你想要做什么?预期的产出是多少?您取而代之的是什么?尽量减少这个最简单的例子,你可以。另外,我不清楚你是否看到错误信息,但如果是的话,它们是什么? – nekomatic

回答

0

这里的一些代码,将通过fileName读取指定的PDB文件,并写入找到的每个配位体ID的输出文件,其中包含从输入文件与HETATM开头,且包含NADA线和字母Z在指定位置:

import os 

fileName = 'C:\\path\\to\\file.pdb' 
baseName = os.path.splitext(fileName)[0] # get filename without extension 

# Read the lines from the input file that have 'HETATM' and 'NAD' in the 
# correct positions into a list: 
with open(fileName) as f: 
    linesNAD = [li for li in f if li[0:6] == 'HETATM' and li[17:20] == 'NAD'] 

# Build a set of the chains found in the lines we have read in: 
chains = {li[21] for li in linesNAD if li[21] in 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'} 

# For each chain, create an output file and write the lines for that chain: 
for chain in chains: 
    outfileName = baseName + '_' + chain + '.txt' 
    with open(outfileName, 'w') as o: 
     o.writelines([li for li in linesNAD if li[21] == chain]) 

如果您想处理多个输入文件,陷阱错误等,您将可以添加到此处。请注意列表和集合理解的使用,它们通常是转换或过滤一组数据的非常简洁的方式。如果你需要关于这些结构的更多细节,请查看Python帮助。

我已经测试了这个在Python 3.5上,但它应该在一些早期版本中工作 - 我认为集合理解在2.7中可用,但我不确定何时引入with结构。如果它在你的版本中不起作用,就像你在原始代码中那样打开文件(但是不要忘记以后他们会这么做)

+0

谢谢,我有2000个输入文件。我怎样才能使用目录输入和输出文件?我新的Python中,请建议我 – Jone

+0

我认为,你已经在你的代码中已经做到这一点的方式应该工作。如果它给你带来麻烦,可以在Python命令行或编写一个测试脚本(例如,只是打印每个文件的路径而不是尝试处理它。在Python中尝试一下很容易,如果你是新的,这是最好的学习方式。如果您再次遇到问题,请在此处搜索您的具体问题(因为它可能已被回答),如果找不到,请提出有关该问题的新问题。 – nekomatic

相关问题