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参数查看示例用法！")