美文网首页
Biopython:为生物信息而设计的python

Biopython:为生物信息而设计的python

作者: wo_monic | 来源:发表于2020-12-30 15:12 被阅读0次

    Biopython文档中文版 https://biopython-cn.readthedocs.io/zh_CN/latest/cn/chr02.html

    biopython安装

    pip3 install biopython

    引入Bio模块

    import Bio
    from Bio.Seq import Seq

    新建序列

    seq1 = Seq("ATTCCGTTTGACTGCGTGTTGAAAAACCTGGGGCCCATGCTGCATGGC")

    GC含量

    from Bio.SeqUtils import GC
    GC(seq1)

    序列处理

    注意:下表是从0开始,区间是左闭右开。

    seq1[0:10:1] #获取第1到10个序列,step是1.

    反向

    seq1.complement()

    反向互补

    seq1.reverse_complement()

    转录DNA成mRNA

    mRNA_seq1 = seq1.reverse_complement().transcribe()

    转录翻译DNA成蛋白序列

    protein_seq1 = seq1.translate() #默认是常用table id 1

    biopython的翻译表是基于NCBI.可以根据实际设定table

    protein_seq1 = seq1.translate(table=2,to_stop=True) #设置为线粒体模式,遇到终止密码子则停止
    protein_seq1 = seq1.translate(to_stop=True) #遇到终止密码子则停止
    protein_seq1 = seq1.translate(cds=True) #说明你输入的是完整的cds序列,不含中止密码子,如果出现终止密码子则会抛出异常

    翻译表

    from Bio.Data import CodonTable
    standard_table = CodonTable.unambiguous_dna_by_name["Standard"]
    print(standard_table)

    MutableSeq对象

    Seq对象的序列是只读的。使用MutableSeq可以转换成可变的序列

    from Bio.Seq import MutableSeq
    mutable_seq = MutableSeq("GCCATTGTAATGGGCCGCTGGCCATTTGCCCGACCATTGTAATGGGCCGCTGAAAGGGTGCCCGAGAAAGGGTGCCCGA")
    mutable_seq[5] = "G"
    mutable_seq.remove("T")

    当编辑完成时,可以用toseq()转变为Seq序列。

    new_Seq = mutable_seq.toseq()

    UnknownSeq对象 使用UnknownSeq可以节省内存

    from Bio.Seq import UnknownSeq
    unk1 = UnknownSeq(20) #20个未知序列,默认是用?
    unk2 = UnknownSeq(30,character="N") #指定字符为N

    第4章 序列注释对象

    from Bio.SeqRecord import SeqRecord

    第5章 序列输入和输出

    读取单条序列

    from Bio import SeqIO
    single_seq_record = Bio.SeqIO.read("single.fa","fasta")

    读取多条序列

    from Bio import SeqIO
    for seq_record in SeqIO.parse("Sua.5LTR3.fa","fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

    把结果读取到list

    records = list(SeqIO.parse("unknown.lib.fa","fasta"))

    使用文件句柄

    with open("test.fa") as handle:
    print(sum(len(r) for r in SeqIO.parse(handle,"fasta")))

    打开gzip压缩文件

    import gzip
    from Bio import SeqIO
    with gzip.open("chr.fa.gz","rt") as handle:
    print(sum(len(r) for r in SeqIO.Parse(handle,"fasta")))

    打开bzip2压缩文件

    import bz2
    from Bio import SeqIO
    with bz2.open("cha.fa.bz2","rt") as handle:
    print(sum(len(r) for r in SeqIO.Parse(handle,"fasta")))

    写入序列文件

    from Bio.Seq import Seq
    from Bio.SeqRecord import SeqRecord
    from Bio import SeqIO
    SeqIO.write(proteins,"translations.fasta","fasta")

    模块化1:把DNA序列转换为蛋白质序列

    修改最后一行的输入文件和输出文件即可运行

    #转换DNA序列为protein序列
    def dna2protein(in_fasta,out_fasta,table=1):
        ##定义转换函数
        from Bio.SeqRecord import SeqRecord
        def make_protein_record(dna_record):
            """return a new SeqRecord with the translate sequence"""
            return SeqRecord(seq = dna_record.seq.translate(table=table),id = dna_record.id,description = "translation of protein ")
        ##读取fasta序列
        from Bio import SeqIO
        proteins = (make_protein_record(dna_record) for dna_record in SeqIO.parse(in_fasta,"fasta"))
        SeqIO.write(proteins,out_fasta,"fasta")
    dna2protein("unknown.fa","unknown.protein.fa")
    

    模块化2:转换处理fasta文件的id

    farename.py脚本代码如下,需要根据实际修改split_head函数里的命令

    #!/usr/bin/env python3
    # -*- coding: utf-8 -*-
    
    """
    用于修改fasta文件的头部,使用前先根据实际修改split_str函数的切割命令
    """
    __version__     = "0.1"
    __date__        = "2021-02-26"
    import sys
    import os
    import gzip
    import Bio
    
    
    #读取用户的输入
    args = sys.argv
    
    #定义对description切片的语法(每次运行前,先根据实际修改切片语法)
    def split_head(description):
        #指定description的第4个,用=分割1次,取第1个
        split_str=description.split()[3].split("=",1)[1]
        return split_str
    
    #读取普通fasta文件,修改文件头部
    def re_header():
        '''
        #读取fasta,并更改header
        '''
        from Bio import SeqIO
        #读取fasta
        for seq_record in SeqIO.parse(InputFa,"fasta"):
            #print(seq_record.id)原始的fasta的id
            #print(seq_record.seq)
            #print(len(seq_record))
            #print(seq_record.description) 原始的fasta的id后的内容
            #指定description的第4个,用=分割1次,取第1个
            #此处的分割需要根据实际进行修改
            description_id=split_head(seq_record.description)
            liststr= ['>',description_id]
            new_id= ''.join(liststr)
            print(new_id)
            print(seq_record.seq)
    
    #读取fasta.gz的文件,修改文件头部
    def re_gz_header():
        '''
        read file.fasta.gz  
        '''
        from Bio import SeqIO
        with gzip.open(InputFa, "rt") as handle:
            for seq_record in SeqIO.parse(handle, "fasta"):
                #print(seq_record.id)
                #print(seq_record.seq)
                #print(len(seq_record))
                #指定description的第4个,用=分割1次,取第1个
                #此处的分割需要根据实际进行修改
                description_id=split_head(seq_record.description)
                liststr= ['>',description_id]
                new_id= ''.join(liststr)
                print(new_id)
                print(seq_record.seq)
    
    
    #判断用户输入的是fasta还是fasta.gz
    def judge_Input():
        if InputFa.endswith(".gz"):
            re_gz_header()
        else:
            re_header()
    
    #判断用户输入
    if len(args) >1:
        if args[1] == "-h" or args[1] == "--help":
            print("Usage:(输入文件支持fasta或fasta.gz)\n"
            "用于修改fasta文件的头部\n"
            "python3 farename.py test.fa #使用切片输出新的fasta.\n"
            "python3 farename.py test.fa.gz #使用切片输出新的fasta.\n"
            "Notes:切片前先检测你的fasta文件的头部的格式,之后修改函数split_head里的切片参数")
        elif len(args) == 2:
            InputFa = args[1]
            judge_Input()
        else:
            print("输入参数有误,请使用-h参数查看示例用法!")
    else:
        print("输入参数有误,请使用-h参数查看示例用法!")
    

    相关文章

      网友评论

          本文标题:Biopython:为生物信息而设计的python

          本文链接:https://www.haomeiwen.com/subject/clsfoktx.html