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参数查看示例用法!")
网友评论