【问题标题】:How to extract chains from a PDB file?如何从 PDB 文件中提取链?
【发布时间】:2012-07-27 10:09:49
【问题描述】:

我想从 pdb 文件中提取链。我有一个名为 pdb.txt 的文件,其中包含 pdb ID,如下所示。前四个字符代表 PDB ID,最后一个字符是链 ID。

1B68A 
1BZ4B
4FUTA

我想 1) 逐行读取文件 2)从对应的PDB文件中下载每条链的原子坐标。
3) 将输出保存到文件夹中。

我使用以下脚本来提取链。但是这段代码只打印来自 pdb 文件的 A 链。

for i in 1B68 1BZ4 4FUT
do 
wget -c "http://www.pdb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId="$i -O $i.pdb
grep  ATOM $i.pdb | grep 'A' > $i\_A.pdb
done

【问题讨论】:

    标签: python bash bioinformatics biopython


    【解决方案1】:

    以下 BioPython 代码应该能很好地满足您的需求。

    它使用PDB.Select 仅选择所需的链(在您的情况下为一条链),并使用PDBIO() 创建仅包含该链的结构。

    import os
    from Bio import PDB
    
    
    class ChainSplitter:
        def __init__(self, out_dir=None):
            """ Create parsing and writing objects, specify output directory. """
            self.parser = PDB.PDBParser()
            self.writer = PDB.PDBIO()
            if out_dir is None:
                out_dir = os.path.join(os.getcwd(), "chain_PDBs")
            self.out_dir = out_dir
    
        def make_pdb(self, pdb_path, chain_letters, overwrite=False, struct=None):
            """ Create a new PDB file containing only the specified chains.
    
            Returns the path to the created file.
    
            :param pdb_path: full path to the crystal structure
            :param chain_letters: iterable of chain characters (case insensitive)
            :param overwrite: write over the output file if it exists
            """
            chain_letters = [chain.upper() for chain in chain_letters]
    
            # Input/output files
            (pdb_dir, pdb_fn) = os.path.split(pdb_path)
            pdb_id = pdb_fn[3:7]
            out_name = "pdb%s_%s.ent" % (pdb_id, "".join(chain_letters))
            out_path = os.path.join(self.out_dir, out_name)
            print "OUT PATH:",out_path
            plural = "s" if (len(chain_letters) > 1) else ""  # for printing
    
            # Skip PDB generation if the file already exists
            if (not overwrite) and (os.path.isfile(out_path)):
                print("Chain%s %s of '%s' already extracted to '%s'." %
                        (plural, ", ".join(chain_letters), pdb_id, out_name))
                return out_path
    
            print("Extracting chain%s %s from %s..." % (plural,
                    ", ".join(chain_letters), pdb_fn))
    
            # Get structure, write new file with only given chains
            if struct is None:
                struct = self.parser.get_structure(pdb_id, pdb_path)
            self.writer.set_structure(struct)
            self.writer.save(out_path, select=SelectChains(chain_letters))
    
            return out_path
    
    
    class SelectChains(PDB.Select):
        """ Only accept the specified chains when saving. """
        def __init__(self, chain_letters):
            self.chain_letters = chain_letters
    
        def accept_chain(self, chain):
            return (chain.get_id() in self.chain_letters)
    
    
    if __name__ == "__main__":
        """ Parses PDB id's desired chains, and creates new PDB structures. """
        import sys
        if not len(sys.argv) == 2:
            print "Usage: $ python %s 'pdb.txt'" % __file__
            sys.exit()
    
        pdb_textfn = sys.argv[1]
    
        pdbList = PDB.PDBList()
        splitter = ChainSplitter("/home/steve/chain_pdbs")  # Change me.
    
        with open(pdb_textfn) as pdb_textfile:
            for line in pdb_textfile:
                pdb_id = line[:4].lower()
                chain = line[4]
                pdb_fn = pdbList.retrieve_pdb_file(pdb_id)
                splitter.make_pdb(pdb_fn, chain)
    

    最后一点:不要为 PDB 文件编写自己的解析器。格式规范很丑(真的很丑),而且有缺陷的 PDB 文件数量惊人。使用像 BioPython 这样的工具来为您处理解析!

    此外,您应该使用与 PDB 数据库交互的工具,而不是使用 wget。他们将 FTP 连接限制、PDB 数据库不断变化的性质等因素考虑在内。我应该知道 - 我updated Bio.PDBList 负责数据库中的更改。 =)

    【讨论】:

    • 非常感谢您的代码和解释。但我不知道,如何运行你的代码?你能帮帮我吗?
    • 当然可以。我添加了一种特别适合您目标的方法。将正文中的行更改为您希望文件转到的目录。 Install Python(你可能已经有了,试试$ which python),然后安装BioPython。以.py 扩展名(例如extract.py)保存上述文件,然后运行$python extract.py pdb.txt。就是这样!
    • 如果你正在做更多与生物信息学相关的事情,我强烈建议你学习 Python。它在该领域非常流行(BioPythonPyMOL,是两个很好的例子),它是一种很好的通用语言。 Python docsThink Python 都是不错的起点。
    【解决方案2】:

    回答这个问题可能有点晚了,但我会给出我的意见。 Biopython 有一些非常方便的功能,可以帮助您轻松实现这样的想法。您可以使用自定义选择类之类的东西,然后使用原始 pdb 文件在 for 循环中为要选择的每个链调用它。

        from Bio.PDB import Select, PDBIO
        from Bio.PDB.PDBParser import PDBParser
    
        class ChainSelect(Select):
            def __init__(self, chain):
                self.chain = chain
    
            def accept_chain(self, chain):
                if chain.get_id() == self.chain:
                    return 1
                else:          
                    return 0
    
        chains = ['A','B','C']
        p = PDBParser(PERMISSIVE=1)       
        structure = p.get_structure(pdb_file, pdb_file)
    
        for chain in chains:
            pdb_chain_file = 'pdb_file_chain_{}.pdb'.format(chain)                                 
            io_w_no_h = PDBIO()               
            io_w_no_h.set_structure(structure)
            io_w_no_h.save('{}'.format(pdb_chain_file), ChainSelect(chain))
    

    【讨论】:

    • 如何在链文件中保留二级结构?这只会从每个链中获取原子,并且您会丢失 HELIX 和 SHEET 注释。
    • 我能想到的唯一方法是使用 get_header() 方法从 PDBParser 对象中提取该信息,然后将其添加到该文件稍后。这可能有助于biopython.org/docs/1.75/api/…
    【解决方案3】:

    假设您有以下文件 pdb_structures

    1B68A 
    1BZ4B
    4FUTA
    

    然后将您的代码放在 load_pdb.sh 中

    while read name
    do
        chain=${name:4:1}
        name=${name:0:4}
        wget -c "http://www.pdb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId="$name -O $name.pdb
        awk -v chain=$chain '$0~/^ATOM/ && substr($0,20,1)==chain {print}' $name.pdb > $name\_$chain.pdb
        #   rm $name.pdb
    done
    

    如果您不需要原始 pdb,请取消注释最后一行。
    执行

    cat pdb_structures | ./load_pdb.sh
    

    【讨论】:

    • 这个简单有效,但它忽略了HETATM记录,它们是链的一部分。既然已经有很多解析器,为什么还要为复杂的文件格式编写自己的解析器?
    • 另外,使用 PDB 检索程序比使用 wget 手动检索 URL 更好。除其他功能外,压缩是自动处理的。请参阅我回答的最后一段。
    • @zelleke,非常感谢您的回答。我尝试了您的代码。我不需要 pdb 文件中的所有链。我已经在文本文件中提到了链 ID。例如,我只需要来自“1B68”的“A”链。
    猜你喜欢
    • 1970-01-01
    • 2016-08-28
    • 1970-01-01
    • 2015-06-15
    • 1970-01-01
    • 1970-01-01
    • 2019-01-05
    • 2022-12-14
    • 1970-01-01
    相关资源
    最近更新 更多