20180711-统计PDB中的蛋白质种类、膜蛋白文件个数及信息等

20180710完成这份工作。简单，但是完成了还是很开心。在我尝试如何使用pickle保存数据后，尝试保存PDB文件中“HEADER”中的信息。文件均保存于实验室服务器（97.73.198.168）/home/RaidDisk/BiodataTest/zyh_pdb_test/tests路径下。本文将记录流程并分享统计结果。

先插入一段代码看看

from PDBParseBase import PDBParserBase

import os

import json

import datetime

from DBParser import ParserBase

import pickle 

def length_counter(seqres_info):

    pdb_id = seqres_info['pdb_id']

    number = 0

    for item in seqres_info.keys():

        if item != 'pdb_id' and item != 'SEQRES_serNum':

            number += int(seqres_info[item]['SEQRES_numRes'])

        pass

    #print(pdb_id)

    #print(number)

    id_and_lenth = []

    id_and_lenth.append(pdb_id)

    id_and_lenth.append(number)

    return id_and_lenth

def find_all_headers(rootdir,saveFilePath,saveFilePath2):

    #with the help of PDBParserBase,just put HEADER inf into a pickle with the form of list

    pdbbase = PDBParserBase()

    start = datetime.datetime.now()

    count = 0

    f = open(saveFilePath,'wb')

    f2 = open(saveFilePath2,'wb')

    header_info_all = []

    for parent,dirnames,filenames in os.walk(rootdir):

        #print(dirnames)

        #print(filenames)

        for filename in filenames:

            count = count + 1

            try:

                PDBfile = filename

                #print(filename)

                header_info = pdbbase.get_header_info(os.path.join(parent,filename))

                #print(header_info)

                header_info_all.append(header_info)

                if count %1000== 0:

                    print(count)

                    end = datetime.datetime.now()

                    print (end-start)

                pass

            except:

                print(filename)

                end = datetime.datetime.now()

                print (end-start)

                pickle.dump(filename, f2)

                #pdb1ov7.ent

            else:

                if count %1111 == 0:

                    print(count)

    pickle.dump(header_info_all, f,protocol=2)

    end = datetime.datetime.now()

    print (end-start)

    print("Done")

    return header_info_all    

def hd_ctr_clsfcsh(listpath):

    #header_counter_classfication

    #return a dic ,the key is the classfication and the value is the num

    with open(listpath, 'rb') as f:

        header_list = pickle.load(f)

    #print(header_list)

    dic = {}

    for item in header_list:

        classification =  item['HEADER_classification']

        print (classification)

        if 'MEMBRAN' in classification:

            if classification in dic.keys():

                dic[classification] = dic[classification] +1

            else:

                dic[classification] = 1

        else:

            pass

    dict= sorted(dic.items(), key=lambda d:d[1], reverse = True)

    print(dict)

    return dict

def get_filenames(listpath,keyword,resultpath):

    with open(listpath, 'rb') as f:

        header_list = pickle.load(f)

    filenames = []

    for item in header_list:

        classification =  item['HEADER_classification']

        if keyword in classification:

            #print(item)

            #print (item['pdb_id'])

            #print (classification)

            id = item['pdb_id']

            filenames.append(id)

        else:

            #print ('dad')

            pass

    print(len(filenames))

    print(filenames)

    f2 = open(resultpath,'wb')

    pickle.dump(filenames, f2,protocol=2)

    print("done")

    return filenames

if __name__ == "__main__":

    rootdir = "/home/BiodataTest/updb"

    #1 is to save list 2 is to save some wrang message

    saveFilePath = "/home/BiodataTest/test_picale/header_counter.pic"

    saveFilePath2 = "/home/BiodataTest/test_picale/header_counter_wrang.pic"

    #all_header_list = find_all_headers(rootdir,saveFilePath,saveFilePath2)

    HEADER_LIST_FILE = "/home/BiodataTest/test_picale/header_counter_list.pic"

    hd_ctr_clsfcsh_dic = hd_ctr_clsfcsh(HEADER_LIST_FILE)

    resultpath = "/home/BiodataTest/test_picale/Membrane_Filename_list.pic"

    #wanted_filenames = get_filenames(HEADER_LIST_FILE,'MEMBRANE',resultpath)

昨天写的代码，第二天就忘记了。

函数一：find_all_headers(rootdir,saveFilePath,saveFilePath2)，输入文件路径，找到PDB文件中所有的“HEADER”信息，其中包括，文件日期、蛋白质种类、蛋白质名字。存入pickle文件中，保存成一个list。难点在于pickle的第三个参数的值是“2”。

函数二：hd_ctr_clsfcsh(listpath):#header_counter_classfication，将刚才输入的list输入，统计出每个分类的数量。输出的是一个字典，键是分类名字，值是该分类的数量。

函数三：get_filenames(listpath,keyword,resultpath)，根据关键词获取给定文件的文件名。在这里，我将所有分类这个信息中含有“MEMBRANE”的蛋白质名字找到，保存成列表输出到文件中。（蛋白质名字就是文件名字）

枯燥的代码介绍完了，来点好看的：大饼图。

这个饼状图画出了PDB数据库中蛋白质种类的分布，其实是不准确的，比如有的蛋白分别属于膜蛋白和转运蛋白。但是标注为“MEMBRANE PROTEIN, TRANSPORT PROTEIN”，那我们把它归为一类。

我们统计了140946个文件，总共有96(100>x）+400（100>x>9）+1606(10>x>1) +1863(x = 1) = 2965个种类。我们选取前1%看看它们分别是什么：

就是它们。总共105057，占总量140946的74.53%。也就是前百分之一的种类占了百分之七十五的数据量。（看起来好残酷的样子哦）

第二部分是有关自己的课题，膜蛋白究竟有多少呢？额，1836个，凄惨的很。所以我扩大了搜索范围，只要包含了“MEMBRANE”的就当作数据提取出来了。注意这个单词的末尾有个“E”，没有"E"的话，数量会增加5个。因为有五个文件叫做“MEMBRANCE PROTEIN”.

找到了它们之后，我们总共获得2335个膜蛋白文件名字，之后将它们解压好，放在预期的文件夹里去。

from PDBParseBase import PDBParserBase

import os

import json

import datetime

from DBParser import ParserBase

#import DataStorage as DS

import time, os,datetime,logging,gzip,configparser,pickle

def  UnzipFile(fileName,unzipFileName):

    """"""

    try:

        zipHandle = gzip.GzipFile(mode='rb',fileobj=open(fileName,'rb'))

        with open(unzipFileName, 'wb') as fileHandle:

            fileHandle.write(zipHandle.read())

    except IOError:

        raise("Unzip file failed!")  

def mkdir(path):

    #Created uncompress path folder

    isExists=os.path.exists(path)

    if not isExists:

        os.makedirs(path)

        print(path + " Created folder sucessful!")

        return True

    else:

        #print ("this path is exist")

        return False   

if __name__ == "__main__": 

    #rootdir = "/home/BiodataTest/updb"

    #rootdir = "/home/BiodataTest/pdb"

    rootdir = "/home/BiodataTest/membraneprotein"

    count = 0

    countcon = 0

    start = datetime.datetime.now()

    saveFilePath = "/home/BiodataTest/picale_zyh_1000/picale.pic"

    #Storage = DS.DataStorage('PDB_20180410')

    mempt_path = "/home/BiodataTest/test_picale/Membrane_Filename_list.pic"

    with open(mempt_path, 'rb') as f:

        memprtn_file = pickle.load(f)

    print(memprtn_file)

    count = 0

    for parent,dirnames,filenames in os.walk(rootdir):

        for filename in filenames:

            count = count + 1

            #start unzip,get the target name and make a files

            dirname = filename[4:6]

            filename_for_membrane_1 = filename[3:7]

            filename_for_membrane = filename_for_membrane_1.upper()

            print(filename_for_membrane)

            if(filename_for_membrane in memprtn_file):

                print(filename_for_membrane)

                filename_with_rootdir = "/home/BiodataTest/pdb/pdb/" + str(dirname)+"/"+str(filename)

                unzipFileName = "/home/BiodataTest/membraneprotein/" + str(dirname)+"/"+str(filename)

                mkdir("/home/BiodataTest/membraneprotein/" + str(dirname)+"/")

                try:

                    UnzipFile(filename_with_rootdir,unzipFileName[:-3])

                    pass

                except:

                    continue

            else:

                #print(filename)

                pass

            pass

            #print(filename)

    end = datetime.datetime.now()

    print("alltime = ")

    print (end-start)

    print(count)

print("Done")

于是就获得了这些膜蛋白。

第三十行“MEMBRANCE”的多字母变体英文就大剌剌的显示在那里。。。总结工作也是需要做好的，不然可能之后会忘记吧。

巴特西

20180711-统计PDB中的蛋白质种类、膜蛋白文件个数及信息等

最新文章

热门文章