#分子模拟#使用shell脚本分割pdb文件中蛋白配体

上期(原创蛋白配体分割工具​)使用python创建来了一个分割由Pymol生成的pbd文件的工具.本期使用shell脚本实现相同的效果.

原始pdb文件是类似这样的:

Pymol生成的pdb文件

在这个pdb文件中,TER开头且不以h结尾的行是真正的配体之间的分隔符(图中蓝色箭头),TER开头以h结尾的行是假的分隔符(红色叉).所以图中的pdb文件分成了四部分,第一部分是第一行开始到第一个分隔符属于蛋白,剩下被分隔符分割的部分属于各个配体文件.

下面以一个pbd文件为例,这个pdb文件可以在(https://pan.baidu.com/s/1bOINmu)

下载.shell脚本如下:

#! /bin/bash
# 将这个pdb文件的总行数赋值给变量.
pro_rows=`wc -l protein.pdb|awk '{print $1}' -`
# 首先获得匹配到以TER开头且不以h结尾的行,记录行号并管道给下个awk.
awk 'BEGIN{OFS="\t"; }{if($0~ "^TER.*" && !($0~ "h\r"))print NR}' protein.pdb | awk -v pro_rows=$pro_rows '{
    if(NR==FNR){a[NR]=$1}    # 将分隔符的行号存储到数组
    if(NR>FNR){b[FNR]=$0}    # 将pdb文件每一行存入数组
}END{
    a[0]=1;    #pbd文件首行存入数组
    la=length(a);a[la]=pro_rows;    #pdb文件尾行存入数组
    # 注意,这里有个坑爹的地方,在awk的for循环里length函数貌似无效.所以先赋值给变量然后再在for循环里调用.
    la=length(a);
    #for循环嵌套输出结果.
    for(i=0;i<la;i++){
        for(j=a[i];j<a[i+1];j++){
            print b[j] > "result_"i".txt"
        }
    }
}' - protein.pdb 

结果会生成11个文件.result_0.txt是蛋白文件,其他是配体文件.

分割成的11个文件

脚本中的protein.pdb是文件名,需要处理其他文件直接替换文件名即可.

更多原创精彩内容敬请关注生信杂谈：