【3】数据分析-9-biopython实践

这两天处理数据的时候,需要解析Blast的结果。经管之前写过解析blast的xml结果的程序,用起来还是不够优雅。刚好,借这次机会来接触一下biopython。

一、简介

biopython是基于python语言来帮助生物信息学工作这解决问题的工具。它可以做:

  1. 将生物信息学文件解析为Python可用的数据结构,包含以下支持的格式:
    • Blast输出结果 – standalone和在线Blast
    • Clustalw
    • FASTA
    • GenBank
    • PubMed和Medline
    • ExPASy文件, 如Enzyme和Prosite
    • SCOP, 包括‘dom’和‘lin’文件
    • UniGene
    • SwissProt
  2. 被支持格式的文件可以通过记录来重复或者通过字典界面来索引。
  3. 处理常见的生物信息学在线数据库的代码:

    • NCBI – Blast, Entrez和PubMed服务
    • ExPASy – Swiss-Prot和Prosite条目, 包括Prosite搜索
  4. 常见生物信息学程序的接口,例如:

    • NCBI的Standalone
    • Clustalw比对
    • EMBOSS命令行工具
  5. 一个能处理序列、ID和序列特征的标准序列类。

  6. 对序列实现常规操作的工具,如翻译,转录和权重计算。

  7. 利用k最近邻接、Bayes或SVM对数据进行分类的代码。

  8. 处理比对的代码,包括创建和处理替换矩阵的标准方法。

  9. 分发并行任务到不同进程的代码。

  10. 实现序列的基本操作,翻译以及BLAST等功能的GUI程序。

  11. 使用这些模块的详细文档和帮助,包括此文件,在线的wiki文档,网站和邮件列表。

  12. 整合BioSQL,一个也被BioPerl和BioJava支持的数据库架构。 感觉好厉害的样子,还是从解析xml开始学起来吧。

biopython下载地址:http://biopython.org/wiki/Download

biopython官网说明:http://biopython-cn.readthedocs.io/zh_CN/latest/index.html

二、实战

2.1.解析xml文件

blastrecord

from Bio.Blast import NCBIXML
def test_bio_python():
    result_file = "tblastn_100000_2.xml"
    output_file = "tblastn_100000_gene.tsv"

    # with open(output_file,"w") as data1:
    data2 = open(output_file, "w")
    data2.write("seq\tlength\tgen\tmap_len\tmismatch\n")

    result_handle = open(result_file)
    blast_records = NCBIXML.parse(result_handle)
    blast_records = list(blast_records)

    for blast_record in blast_records:
       # blast_record为一条序列的比对结果
        # for description in blast_record.descriptions:
        #     print description.title
        #     print description.score
        #     print description.e
        #     print description.num_alignments
        # break
        # print blast_record.
        # print blast_record.description.title

       # print dir(blast_record) 可以看到每条blast_record的基本属性
        query_id = blast_record.query  # query的名字
        query_length = blast_record.query_length
        # print blast_record.reference
        # print blast_record.query_letters
        # print blast_record.query_id

        for aligment in blast_record.alignments: 
        #aligment为一条序列比对到一个基因,其中一个基因有可能被多次Hits
            for hsp in aligment.hsps:
                # print hsp.__getattr__

                # print aligment.title #被Hit的序列的名字,相当于xml中的
                #<Hit_def>Homo sapiens ER membrane protein complex subunit 1 (EMC1), transcript variant 2, mRNA</Hit_def>

                #获得()中的基因名字
                if "(" in aligment.title:
                    if aligment.title.count("(") == 1:
                        gene = aligment.title.split("(")[1].split(")")[0]
                    elif aligment.title.count("(") == 2:
                        gene = aligment.title.split(")")[1].split("(")[1]
                else:
                    gene = aligment.title

                # print aligment.length   #比对序列的长度


                # print hsp.score
                # print hsp.bits
                # print hsp.expect
                # print hsp.num_alignments
                # print hsp.identities
                identity_len = hsp.identities  #比对到的碱基个数,相当于<Hsp_identity>4</Hsp_identity>

                # print hsp.positives  #相当于<Hsp_positive>6</Hsp_positive>
                # print hsp.gaps  #相当于<Hsp_gaps>0</Hsp_gaps>
                # print hsp.strand
                # print hsp.frame
                # print hsp.query #相当于<Hsp_qseq>LLVDDQ</Hsp_qseq>
                # print hsp.query_start
                # print hsp.match  #相当于<Hsp_midline>LL+DD+</Hsp_midline>
                # print hsp.sbjct  #相当于<Hsp_hseq>LLIDDE</Hsp_hseq>
                # print hsp.sbjct_start

                gaps = hsp.gaps
                identity_len = hsp.identities

                align_len = len(query_seq) - int(gaps)
                identity_1 = int(identity_len) / int(align_len) * 100 #比对上的序列的相似性
                identity_2 = int(identity_len) / int(query_length) * 100 #全长比对的百分比
                mismatch = int(query_length) - int(identity_len)
                result = [query_id, str(query_length), gene, str(identity_len),
                          str(mismatch)]
                data2.write("%s\n" % "\t".join(result))
                # break
    data2.close()

后面用到,再来接着补充吧

2.蛋白pdb文件的分析

下载蛋白序列文件

from Bio.PDB.PDBList import PDBList
def download_pdb(pdb_id, download_dir='pdb'):
    pdbl = PDBList()
    pdf_name = os.path.join(download_dir, 'pdb%s.ent' % pdb_id.lower())
    if not os.path.exists(pdf_name):
        pdbl.retrieve_pdb_file(pdb_id, pdir=download_dir, file_format='pdb')
    return pdf_name

根据指定位置提取残基的向量

from Bio.PDB.PDBParser import PDBParser
def get_vector(pdb_id, model, chain, chain_s, chain_e, download_dir='pdb'):
    structure_info = []
    parser = PDBParser()
    pdf_name = download_pdb(pdb_id, download_dir)
    chain_info = parser.get_structure("test", pdf_name)[model][chain]

    #for ii in range(chain_s, chain_e + 1):

        # 注意这里了,这个地方我一直报错,原因在于hetatm这一行读为residue以后;residue的key变成了('H_MSE', ii, ' '),而不再是(' ', ii, ' '),所以,就不能直接用chain_info[ii]来提取了,当然如果这个残基不是MSE,就需哟啊灵当别论了

    #   try :
    #       atom = chain_info[ii]['CA']
    #   except:
     #      atom = chain_info[('H_MSE', ii, ' ')]['CA']

      #  atom = chain_info[ii]['CA']

       # vector_info = atom.get_vector()
        #structure_info.append(list(vector_info))

     locs_range = range(chain_s, chain_e + 1)

        for one_residue in chain_info:
            residue_id = one_residue.get_id()[1]
            if residue_id in locs_range :
                atom = one_residue['CA']
                vector_info = atom.get_vector()
                structure_info.append(list(vector_info))


    return structure_info

3.修改指定参加的所有原子的向量

这里面涉及到怎么通过biopython写PDB文件,怎么提取指定的atom,并如何修改atom,短短几行代码,做了如此多的事情,精妙呀

from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.PDBList import PDBList
from Bio.PDB.PDBIO import PDBIO,Select

from Bio.PDB import extract

class SelectAtom(Select):
    def __init__(self,chain_s,chain_e,new_vector):
        self.chain_s = chain_s
        self.chain_e = chain_e
        self.new_vector = new_vector

    def accept_residue(self,residue):
        id_name = residue.get_id()[1]
        if id_name < self.chain_s or id_name > self.chain_e:
            return False
        else:
            for atom in residue :
                atom.set_coord(self.new_vector)
            return True

def get_atom(pdb_id, model, chain, chain_s, chain_e, download_dir='pdb'):
    parser = PDBParser()
    pdf_name = download_pdb(pdb_id, download_dir)

    chain_info = parser.get_structure("test", pdf_name)[model][chain]
    new_list = [2,2,2]
    io = PDBIO()
    io.set_structure(chain_info)
    io.save('out.pdb',SelectAtom(chain_s,chain_e,new_list))

2.3 读fasta序列

from Bio import SeqIO
fastas = SeqIO.parse(fasta_file, 'fasta')
for fasta in fastas:
    print fasta.id
    print fasta.seq

2.4 核酸序列转氨基酸序列

from Bio.Seq import Seq
seq = Seq(test_seq)
print seq.translate()

2.5 PDB提取某些链,并修改链的名字

# 提取链,并修改名字
class SelectChain(Select):
    def __init__(self,chain_hash):
        self.chain_hash = chain_hash

    def accept_chain(self,chain):
        id_name = chain.get_id()
        if id_name not in self.chain_hash:
            return False
        else:
            chain.id = self.chain_hash[id_name]  # 修改链的名字
            return  True

# 只提取链
class SelectChain_2(Select):
    def __init__(self,chain_hash):
        self.chain_hash = chain_hash

    def accept_chain(self,chain):
        id_name = chain.get_id()
        if id_name not in self.chain_hash:
            return False
        else:
            return  True

# pdb_id:  '4XGZ'
# chain_hash:{'A':'H','B':'L'}
def keep_chain(pdb_id,chain_hash,result_dir='./'):
    pdf_name = download_pdb(pdb_id)
    parser = PDBParser()
    chain_info = parser.get_structure("test", pdf_name)
    io = PDBIO()
    io.set_structure(chain_info)
    pdb_result = '%s/%s.pdb'% (result_dir,pdb_id)
    pdb_result_2 = '%s/%s_2.pdb'% (result_dir,pdb_id)
    try:
        io.save(pdb_result,SelectChain(chain_hash))
    except:
        io.save(pdb_result_2, SelectChain_2(chain_hash))
        parser_2 = PDBParser()
        chain_info_2 = parser_2.get_structure("test", pdb_result_2)
        io_2 = PDBIO()
        io_2.set_structure(chain_info_2)
        pdb_result = '%s/%s.pdb' % (result_dir, pdb_id)
        io_2.save(pdb_result, SelectChain(chain_hash))
        os.system('rm %s' % pdb_result_2)

注:先尝试提取链的时候,也修改链的名字,但是有的PDB修改后的名字已经存在其他链中,就会报错:

ValueError: Cannot change id from `A` to `h`. The id `h` is already used for a sibling of this entity.

所以,先提取,提取以后再修改

参考资料

个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn

Sam avatar
About Sam
专注生物信息 专注转化医学