用 Python 玩转常用生物序列

发布于 2021-04-29 07:33

一、准备工作

1、获取感兴趣的基因,蛋白质,转录本等生物序列 FASTA 或 GenBank

这里举例,进入 NCBI 获取的GeneBank / FASTA 的数据格式

比如查看 POU5F1 基因:https://www.ncbi.nlm.nih.gov/gene/5460

2、搭建 Python 环境与项目目录

现在我们的目录结构是这样的

搭建目录结构及Python环境参考:https://blog.csdn.net/u011262253/article/details/105902060

二、操作生物序列

1、读取常见的序列文件格式(fasta,gb)
 1from Bio import SeqIO
2
3# 读取包含单个序列 Fasta 格式文件
4fa_seq = SeqIO.read("res/sequence1.fasta""fasta")
5# print fa_seq
6# 读取包含多个序列的 fasta 格式文件
7for fa in SeqIO.parse("res/multi.fasta""fasta"):
8    print (fa.seq)
9# 一个多序列文件中的所有序列
10seqs = [fa.seq for fa in SeqIO.parse("res/multi.fasta",  "fasta")]
11print (seqs)
12# 如果不想要seq对象中的字母表,可以用str()来强制类型转换
13seqs = [str(fa.seq) for fa in SeqIO.parse("res/multi.fasta",  "fasta")]
14print (seqs)
15
16# 读取包含单个序列的 gb 格式文件
17gb_seq = SeqIO.read("res/sequence1.gb""genbank")
18print (gb_seq)
2、浏览 fasta 序列文件内容
 1from Bio import SeqIO
2
3# 读取 Fasta 文件详细信息
4fa_seq = SeqIO.read("res/sequence1.fasta""fasta")
5
6# =====获取详细的信息=====
7# 提取基因ID,name
8# Fasta 文件中序列名所在行的第一个词被作为 id 和 name
9print ("id: ", fa_seq.id)
10print ("name: ", fa_seq.name)
11# 基因 Description 是fasta文件格式中的第一行
12print ("description: ", fa_seq.description)
13# 序列
14print ("seq: ", fa_seq.seq)
15# 序列来源库信息(NCBI的数据库信息会包括数据库交叉引用)
16print ("dbxrefs: ", fa_seq.dbxrefs)
17# 全部序列的注释信息
18print ("annotations: ", fa_seq.annotations)
19# 序列中每个字母的注释信息
20print ("letter_annotations: ", fa_seq.letter_annotations)
21# 部分序列的注释信息
22print ("features: ", fa_seq.features)
3、读取 genebank 文件详细信息
 1from Bio import SeqIO
2
3# 读取包含单个序列的 gb 格式文件
4gb_seq = SeqIO.read("res/sequence1.gb""genbank")    
5print (gb_seq)
6
7# =====获取详细的信息=====
8# 提取基因ID,name
9# gb文件中序列名包含比fasta更加详细的序列信息,下面分别是 id 和 name
10print ("id: ", gb_seq.id)
11print ("name: ", gb_seq.name)
12# 基因 Description 是fasta文件格式中的第一行
13print ("description: ",  gb_seq.description)
14# 序列信息, 这里的序列信息是以 bioPython 中的seq对象存储
15print ("seq: ", gb_seq.seq)
16# 序列来源库信息(NCBI的数据库信息会包括数据库交叉引用)
17print ("dbxrefs: ", gb_seq.dbxrefs)
18# 全部序列的注释信息
19print ("annotations: ", gb_seq.annotations)
20# 序列中每个字母的注释信息
21print ("letter_annotations: ", gb_seq.letter_annotations)
22# 部分序列的注释信息,SeqFeature 对象的形式保存了features table中的所有entries(如genes和CDS等)
23print ("features: ", gb_seq.features)
24# 该基因的物种信息
25print ("organism: ", gb_seq.annotations["organism"])
26# 关于序列的注释信息,相关数据库的交叉引用号
27print ("comment: ", gb_seq.annotations["comment"])
28# 序列来源的物种名
29print ("source: ", gb_seq.annotations["source"])
30# 该基因的分类学信息
31print ("taxonomy: ", gb_seq.annotations["taxonomy"])
32# 该基因的整理后的注释信息
33print ("structured_comment: ", gb_seq.annotations["structured_comment"])
34# 该基因序列相关的关键词
35print ("keywords: ", gb_seq.annotations["keywords"])
36# 该基因的相关文献编号,或递交序列的注册信息
37print ("references: ", gb_seq.annotations["references"])
38# 该基因的入库时,给的基因编号,以及在染色体上的位点信息
39print ("accessions: ", gb_seq.annotations["accessions"])
40# 该基因的分子类型,一般为 DNA
41print ("molecule_type: ", gb_seq.annotations["molecule_type"])
42# 该基因的数据文件划分方式
43print ("data_file_division: ", gb_seq.annotations["data_file_division"])
44# 基因发布时间
45print ("date: ", gb_seq.annotations["date"])
46# 该基因的更新版本
47print ("sequence_version: ", gb_seq.annotations["sequence_version"])
48# 该基因的拓扑结构
49print ("topology: ", gb_seq.annotations["topology"])

相信大家可以看到 GeneBank 比 fasta 格式更加详细和贴心,但是对于大量待处理的序列来说内存占用和运行时间比这些详细信息更加重要。这就使fasta成为我们一般在序列分析中常用的格式。

4、新建序列文件
1from Bio.Seq import Seq
2
3# 新建一个DNA序列对象
4dna_seq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_dna)
5# 新建一个RNA序列对象
6rna_seq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_rna)
7# # 新建一个蛋白质序列对象
8protein_seq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.protein)

序列对象由一段字符串和其对应的编码表所定义。我们可以从上述的代码中看到,字符串内容一样,唯一不同的就是第二个参数IUPAC值不一样。IUPAC (International Union of Pure and Applied Chemistry ) 是一个制定化学相关标准的组织,Biopython 所使用的编码表就是由它制定的,想了解详细细节可以参考http://www.bioinformatics.org/sms2/iupac.html ,详细定义如下:

名称编码表
ambiguous_dna_lettersGATCRYWSMKHBVDN
unambiguous_dna_lettersGATC
ambiguous_rna_lettersGAUCRYWSMKHBVDN
unambiguous_rna_lettersGAUC
proteinARNDCQEGHILKMFPSTWYV
5、修改序列文件

在生物学意义上,序列是不可以随便更改的,也就是不可变的。如果强行修改,那么就会报错TypeError: 'Seq' object does not support item assignment

1dna_seq[0] = "G"

如果你一定要修改,当然也可以,但不建议这么做

1dna_seq_mutable = dna_seq.tomutable()
2dna_seq_mutable[0] = "G"
6、操作序列文件
 1from Bio.Seq import Seq
2
3# 新建一个DNA序列对象
4dna_seq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_dna)
5# 序列信息
6print ("Sequence: ", dna_seq)
7# 序列长度
8print ("Length : ", len(dna_seq))
9# 单个核苷酸计数
10print (