Biopython - 快速指南

Biopython - 简介

Biopython 是最大且最流行的 Python 生物信息学软件包。它包含许多用于常见生物信息学任务的不同子模块。它由 Chapman 和 Chang 开发，主要用 Python 编写。它还包含用于优化软件复杂计算部分的 C 代码。它运行在 Windows、Linux、Mac OS X 等上。

基本上，Biopython 是 Python 模块的集合，提供处理 DNA、RNA 和蛋白质序列操作的功能，例如 DNA 字符串的反向互补、在蛋白质序列中查找基序等。它提供了大量解析器来读取所有主要遗传数据库如 GenBank、SwissPort、FASTA 等，以及在 python 环境中运行其他流行生物信息学软件/工具（如 NCBI BLASTN、Entrez 等）的包装器/接口。它有 BioPerl、BioJava 和 BioRuby 等兄弟项目。

特征

Biopython 是可移植的、清晰的并且具有易于学习的语法。下面列出了一些显着特征 -

  • 解释性、交互性和面向对象。

  • 支持 FASTA、PDB、GenBank、Blast、SCOP、PubMed/Medline、ExPASy 相关格式。

  • 处理序列格式的选项。

  • 管理蛋白质结构的工具。

  • BioSQL - 用于存储序列以及特征和注释的标准 SQL 表集。

  • 访问在线服务和数据库，包括 NCBI 服务（Blast、Entrez、PubMed）和 ExPASY 服务（SwissProt、Prosite）。

  • 访问本地服务，包括 Blast、Clustalw、EMBOSS。

目标

Biopython 的目标是通过 python 语言提供简单、标准和广泛的生物信息学访问。Biopython 的具体目标如下：

  • 提供对生物信息学资源的标准化访问。

  • 高质量、可重用的模块和脚本。

  • 可在集群代码、PDB、NaiveBayes 和马尔可夫模型中使用的快速数组操作。

  • 基因组数据分析。

优点

Biopython 需要的代码非常少，并具有以下优点 -

  • 提供聚类中使用的微阵列数据类型。

  • 读取和写入树视图类型文件。

  • 支持用于PDB解析、表示和分析的结构数据。

  • 支持 Medline 应用程序中使用的期刊数据。

  • 支持BioSQL数据库，这是所有生物信息学项目中广泛使用的标准数据库。

  • 通过提供将生物信息学文件解析为格式特定记录对象或序列加特征的通用类的模块来支持解析器开发。

  • 基于食谱风格的清晰文档。

案例研究样本

让我们检查一些用例（群体遗传学、RNA 结构等），并尝试了解 Biopython 如何在该领域发挥重要作用 -

群体遗传学

群体遗传学是对群体内遗传变异的研究，涉及对群体中基因和等位基因频率随空间和时间变化的检查和建模。

Biopython 提供了用于群体遗传学的 Bio.PopGen 模块。该模块包含收集有关经典群体遗传学信息的所有必要功能。

RNA结构

我们生命所必需的三大生物大分子是DNA、RNA和蛋白质。蛋白质是细胞的主力，并像酶一样发挥着重要作用。DNA（脱氧核糖核酸）被认为是细胞的“蓝图”。它携带细胞生长、吸收营养和繁殖所需的所有遗传信息。RNA（核糖核酸）在细胞中充当“DNA 复印件”。

Biopython 提供代表核苷酸、DNA 和 RNA 构建块的 Bio.Sequence 对象。

Biopython - 安装

本节介绍如何在您的计算机上安装 Biopython。安装非常简单，不会超过五分钟。

第 1 步- 验证 Python 安装

Biopython 设计用于与 Python 2.5 或更高版本配合使用。所以，必须先安装python。在命令提示符中运行以下命令 -

> python --version

它的定义如下 -

验证 Python 安装

如果安装正确，它会显示 python 的版本。否则，请下载最新版本的 python，安装它，然后再次运行该命令。

步骤 2 - 使用 pip 安装 Biopython

在所有平台上使用 pip 从命令行轻松安装 Biopython。输入以下命令 -

> pip install biopython

您的屏幕上将显示以下响应 -

使用 pip 安装 Biopython

用于更新旧版本的 Biopython -

> pip install biopython –-upgrade

您的屏幕上将显示以下响应 -

更新旧版本

执行此命令后，旧版本的 Biopython 和 NumPy（Biopython 依赖于它）将被删除，然后再安装最新版本。

步骤 3 - 验证 Biopython 安装

现在，您已经在计算机上成功安装了 Biopython。要验证 Biopython 是否已正确安装，请在 python 控制台上键入以下命令 -

验证 Biopython 安装

它显示了 Biopython 的版本。

替代方法 - 使用源代码安装 Biopython

要使用源代码安装 Biopython，请按照以下说明操作 -

从以下链接下载最新版本的 Biopython - https://biopython.org/wiki/Download

截至目前，最新版本是biopython-1.72。

下载文件并解压压缩存档文件，进入源代码文件夹并输入以下命令 -

> python setup.py build

这将从源代码构建 Biopython，如下所示 -

使用源安装 Biopython

现在，使用以下命令测试代码 -

> python setup.py test

测试代码

最后，使用以下命令安装 -

> python setup.py install

最后安装

Biopython - 创建简单的应用程序

让我们创建一个简单的 Biopython 应用程序来解析生物信息学文件并打印内容。这将帮助我们理解 Biopython 的一般概念以及它如何在生物信息学领域提供帮助。

步骤 1 - 首先，创建一个示例序列文件“example.fasta”并将以下内容放入其中。

>sp|P25730|FMS1_ECOLI CS1 fimbrial subunit A precursor (CS1 pilin) 
MKLKKTIGAMALATLFATMGASAVEKTISVTASVDPTVDLLQSDGSALPNSVALTYSPAV
NNFEAHTINTVVHTNDSDKGVVVKLSADPVLSNVLNPTLQIPVSVNFAGKPLSTTGITID 
SNDLNFASSGVNKVSSTQKLSIHADATRVTGGALTAGQYQGLVSIILTKSTTTTTTTKGT 

>sp|P15488|FMS3_ECOLI CS3 fimbrial subunit A precursor (CS3 pilin) 
MLKIKYLLIGLSLSAMSSYSLAAAGPTLTKELALNVLSPAALDATWAPQDNLTLSNTGVS 
NTLVGVLTLSNTSIDTVSIASTNVSDTSKNGTVTFAHETNNSASFATTISTDNANITLDK 
NAGNTIVKTTNGSQLPTNLPLKFITTEGNEHLVSGNYRANITITSTIKGGGTKKGTTDKK

扩展名fasta指的是序列文件的文件格式。FASTA源于生物信息学软件FASTA，因而得名。FASTA格式有多个序列一一排列，每个序列都有自己的id、名称、描述和实际的序列数据。

步骤 2 - 创建一个新的 python 脚本 *simple_example.py" 并输入以下代码并保存。

from Bio.SeqIO import parse 
from Bio.SeqRecord import SeqRecord 
from Bio.Seq import Seq 

file = open("example.fasta") 

records = parse(file, "fasta") for record in records:    
   print("Id: %s" % record.id) 
   print("Name: %s" % record.name) 
   print("Description: %s" % record.description) 
   print("Annotations: %s" % record.annotations) 
   print("Sequence Data: %s" % record.seq) 
   print("Sequence Alphabet: %s" % record.seq.alphabet)

让我们更深入地研究一下代码 -

第 1行导入 Bio.SeqIO 模块中可用的解析类。Bio.SeqIO模块用于读写不同格式的序列文件，“parse”类用于解析序列文件的内容。

第 2 行导入 Bio.SeqRecord 模块中可用的 SeqRecord 类。该模块用于操作序列记录，SeqRecord 类用于表示序列文件中可用的特定序列。

*第3行导入Bio.Seq模块中可用的Seq类。该模块用于操作序列数据，Seq类用于表示序列文件中可用的特定序列记录的序列数据。

第 5 行使用常规 python 函数 open 打开“example.fasta”文件。

第7行解析序列文件的内容并将内容作为SeqRecord对象的列表返回。

第 9-15 行使用 python for 循环遍历记录，并打印序列记录 (SqlRecord) 的属性，例如 id、名称、描述、序列数据等。

第 15 行使用 Alphabet 类打印序列的类型。

步骤 3 - 打开命令提示符并转到包含序列文件“example.fasta”的文件夹，然后运行以下命令 -

> python simple_example.py

步骤 4 - Python 运行脚本并打印示例文件“example.fasta”中可用的所有序列数据。输出将类似于以下内容。

Id: sp|P25730|FMS1_ECOLI 
Name: sp|P25730|FMS1_ECOLI 
Decription: sp|P25730|FMS1_ECOLI CS1 fimbrial subunit A precursor (CS1 pilin) 
Annotations: {} 
Sequence Data: MKLKKTIGAMALATLFATMGASAVEKTISVTASVDPTVDLLQSDGSALPNSVALTYSPAVNNFEAHTINTVVHTNDSD
KGVVVKLSADPVLSNVLNPTLQIPVSVNFAGKPLSTTGITIDSNDLNFASSGVNKVSSTQKLSIHADATRVTGGALTA
GQYQGLVSIILTKSTTTTTTTKGT 
Sequence Alphabet: SingleLetterAlphabet() 
Id: sp|P15488|FMS3_ECOLI 
Name: sp|P15488|FMS3_ECOLI 
Decription: sp|P15488|FMS3_ECOLI CS3 fimbrial subunit A precursor (CS3 pilin) 
Annotations: {} 
Sequence Data: MLKIKYLLIGLSLSAMSSYSLAAAGPTLTKELALNVLSPAALDATWAPQDNLTLSNTGVSNTLVGVLTLSNTSIDTVS
IASTNVSDTSKNGTVTFAHETNNSASFATTISTDNANITLDKNAGNTIVKTTNGSQLPTNLPLKFITTEGNEHLVSGN
YRANITITSTIKGGGTKKGTTDKK 
Sequence Alphabet: SingleLetterAlphabet()

在这个例子中我们看到了三个类：parse、SeqRecord 和 Seq。这三个类提供了大部分功能，我们将在下一节中学习这些类。

Biopython - 序列

序列是用于表示生物体的蛋白质、DNA 或 RNA 的一系列字母。它由 Seq 类表示。Seq 类在 Bio.Seq 模块中定义。

让我们在 Biopython 中创建一个简单的序列，如下所示 -

>>> from Bio.Seq import Seq 
>>> seq = Seq("AGCT") 
>>> seq 
Seq('AGCT') 
>>> print(seq) 
AGCT

在这里，我们创建了一个简单的蛋白质序列AGCT，每个字母代表A丙氨酸、G赖氨酸、半胱氨酸和T苏氨酸。

每个 Seq 对象都有两个重要的属性 -

  • 数据 - 实际序列字符串 (AGCT)

  • 字母表- 用于表示序列的类型。例如DNA序列、RNA序列等。默认情况下，它不代表任何序列，本质上是通用的。

字母模块

Seq 对象包含 Alphabet 属性来指定序列类型、字母和可能的操作。它在 Bio.Alphabet 模块中定义。字母表可以定义如下 -

>>> from Bio.Seq import Seq 
>>> myseq = Seq("AGCT") 
>>> myseq 
Seq('AGCT') 
>>> myseq.alphabet 
Alphabet()

Alphabet 模块提供以下类来表示不同类型的序列。Alphabet - 所有类型字母的基类。

SingleLetterAlphabet - 字母大小为 1 的通用字母表。它源自 Alphabet，所有其他字母类型均源自它。

>>> from Bio.Seq import Seq 
>>> from Bio.Alphabet import single_letter_alphabet 
>>> test_seq = Seq('AGTACACTGGT', single_letter_alphabet) 
>>> test_seq 
Seq('AGTACACTGGT', SingleLetterAlphabet())

ProteinAlphabet - 通用单字母蛋白质字母表。

>>> from Bio.Seq import Seq 
>>> from Bio.Alphabet import generic_protein 
>>> test_seq = Seq('AGTACACTGGT', generic_protein) 
>>> test_seq 
Seq('AGTACACTGGT', ProteinAlphabet())

NucleotideAlphabet - 通用单字母核苷酸字母表。

>>> from Bio.Seq import Seq 
>>> from Bio.Alphabet import generic_nucleotide 
>>> test_seq = Seq('AGTACACTGGT', generic_nucleotide) >>> test_seq 
Seq('AGTACACTGGT', NucleotideAlphabet())

DNAAlphabet - 通用单字母 DNA 字母表。

>>> from Bio.Seq import Seq 
>>> from Bio.Alphabet import generic_dna 
>>> test_seq = Seq('AGTACACTGGT', generic_dna) 
>>> test_seq 
Seq('AGTACACTGGT', DNAAlphabet())

RNAAlphabet - 通用单字母 RNA 字母表。

>>> from Bio.Seq import Seq 
>>> from Bio.Alphabet import generic_rna 
>>> test_seq = Seq('AGTACACTGGT', generic_rna) 
>>> test_seq 
Seq('AGTACACTGGT', RNAAlphabet())

Biopython 模块 Bio.Alphabet.IUPAC 提供了 IUPAC 社区定义的基本序列类型。它包含以下类 -

  • IUPACProtein（蛋白质） - 20 个标准氨基酸的 IUPAC 蛋白质字母表。

  • Extended IUPACProtein (extended_ Protein) - 扩展大写 IUPAC 蛋白质单字母字母表，包括 X。

  • IUPACAmbigeousDNA (ambiguously_dna) - 大写 IUPAC 模糊 DNA。

  • IUPACUnambiguousDNA (unambigously_dna) - 大写的 IUPAC 明确 DNA (GATC)。

  • ExtendedIUPACDNA (extended_dna) - 扩展的 IUPAC DNA 字母表。

  • IUPACAmbigouslyRNA (ambiguously_rna) - 大写 IUPAC 模糊 RNA。

  • IUPACUnambiguousRNA (unambigously_rna) - 大写 IUPAC 明确 RNA (GAUC)。

考虑 IUPACProtein 类的一个简单示例，如下所示 -

>>> from Bio.Alphabet import IUPAC 
>>> protein_seq = Seq("AGCT", IUPAC.protein) 
>>> protein_seq 
Seq('AGCT', IUPACProtein()) 
>>> protein_seq.alphabet

此外，Biopython 通过 Bio.Data 模块公开所有生物信息学相关的配置数据。例如，IUPACData. Protein_letters 具有 IUPACProtein 字母表中可能的字母。

>>> from Bio.Data import IUPACData 
>>> IUPACData.protein_letters 
'ACDEFGHIKLMNPQRSTVWY'

基本操作

本节简要说明 Seq 类中可用的所有基本操作。序列类似于 python 字符串。我们可以执行Python字符串操作，如序列中的切片、计数、连接、查找、分割和剥离。

使用以下代码获得各种输出。

获取序列中的第一个值。

>>> seq_string = Seq("AGCTAGCT") 
>>> seq_string[0] 
'A'

打印前两个值。

>>> seq_string[0:2] 
Seq('AG')

打印所有值。

>>> seq_string[ : ] 
Seq('AGCTAGCT')

执行长度和计数操作。

>>> len(seq_string) 
8 
>>> seq_string.count('A') 
2

添加两个序列。

>>> from Bio.Alphabet import generic_dna, generic_protein 
>>> seq1 = Seq("AGCT", generic_dna) 
>>> seq2 = Seq("TCGA", generic_dna)
>>> seq1+seq2 
Seq('AGCTTCGA', DNAAlphabet())

这里，上述两个序列对象 seq1、seq2 是通用 DNA 序列，因此您可以添加它们并生成新序列。您不能添加具有不兼容字母的序列，例如蛋白质序列和 DNA 序列，如下所示 -

>>> dna_seq = Seq('AGTACACTGGT', generic_dna) 
>>> protein_seq = Seq('AGUACACUGGU', generic_protein) 
>>> dna_seq + protein_seq 
..... 
..... 
TypeError: Incompatible alphabets DNAAlphabet() and ProteinAlphabet() 
>>>

要添加两个或多个序列，首先将其存储在 python 列表中，然后使用“for 循环”检索它，最后将其添加在一起，如下所示 -

>>> from Bio.Alphabet import generic_dna 
>>> list = [Seq("AGCT",generic_dna),Seq("TCGA",generic_dna),Seq("AAA",generic_dna)] 
>>> for s in list: 
... print(s) 
... 
AGCT 
TCGA 
AAA 
>>> final_seq = Seq(" ",generic_dna) 
>>> for s in list: 
... final_seq = final_seq + s 
... 
>>> final_seq 
Seq('AGCTTCGAAAA', DNAAlphabet())

在下面的部分中，给出了各种代码以根据要求获得输出。

更改顺序的大小写。

>>> from Bio.Alphabet import generic_rna 
>>> rna = Seq("agct", generic_rna) 
>>> rna.upper() 
Seq('AGCT', RNAAlphabet())

检查 python 成员资格和身份运算符。

>>> rna = Seq("agct", generic_rna) 
>>> 'a' in rna 
True 
>>> 'A' in rna 
False 
>>> rna1 = Seq("AGCT", generic_dna) 
>>> rna is rna1 
False

在给定序列中查找单个字母或字母序列。

>>> protein_seq = Seq('AGUACACUGGU', generic_protein) 
>>> protein_seq.find('G') 
1 
>>> protein_seq.find('GG') 
8

进行分裂操作。

>>> protein_seq = Seq('AGUACACUGGU', generic_protein) 
>>> protein_seq.split('A') 
[Seq('', ProteinAlphabet()), Seq('GU', ProteinAlphabet()), 
   Seq('C', ProteinAlphabet()), Seq('CUGGU', ProteinAlphabet())]

按顺序执行剥离操作。

>>> strip_seq = Seq(" AGCT ") 
>>> strip_seq 
Seq(' AGCT ') 
>>> strip_seq.strip() 
Seq('AGCT')

Biopython - 高级序列操作

在本章中，我们将讨论 Biopython 提供的一些高级序列功能。

补体和反向补体

核苷酸序列可以反向互补以获得新的序列。另外，可以将互补序列反向互补以获得原始序列。Biopython 提供了两种方法来执行此功能 -补足和反向补足。其代码如下 -

>>> from Bio.Alphabet import IUPAC 
>>> nucleotide = Seq('TCGAAGTCAGTC', IUPAC.ambiguous_dna) 
>>> nucleotide.complement() 
Seq('AGCTTCAGTCAG', IUPACAmbiguousDNA()) 
>>>

这里，complement() 方法允许补充 DNA 或 RNA 序列。verse_complement() 方法对结果序列进行补足并从左到右反转。如下所示 -

>>> nucleotide.reverse_complement() 
Seq('GACTGACTTCGA', IUPACAmbiguousDNA())

Biopython使用Bio.Data.IUPACData提供的ambigously_dna_complement变量来进行补码操作。

>>> from Bio.Data import IUPACData 
>>> import pprint 
>>> pprint.pprint(IUPACData.ambiguous_dna_complement) {
   'A': 'T',
   'B': 'V',
   'C': 'G',
   'D': 'H',
   'G': 'C',
   'H': 'D',
   'K': 'M',
   'M': 'K',
   'N': 'N',
   'R': 'Y',
   'S': 'S',
   'T': 'A',
   'V': 'B',
   'W': 'W',
   'X': 'X',
   'Y': 'R'} 
>>>

GC含量

预计基因组 DNA 碱基组成（GC 含量）将显着影响基因组功能和物种生态。GC含量是GC核苷酸的数量除以总核苷酸。

要获取 GC 核苷酸含量，请导入以下模块并执行以下步骤 -

>>> from Bio.SeqUtils import GC 
>>> nucleotide = Seq("GACTGACTTCGA",IUPAC.unambiguous_dna) 
>>> GC(nucleotide) 
50.0

转录

转录是将DNA序列转变为RNA序列的过程。实际的生物转录过程是以DNA为模板链进行反向互补（TCAG → CUGA）以获得mRNA。然而，在生物信息学以及 Biopython 中，我们通常直接使用编码链，并且可以通过将字母 T 更改为 U 来获得 mRNA 序列。

上述的简单示例如下 -

>>> from Bio.Seq import Seq 
>>> from Bio.Seq import transcribe 
>>> from Bio.Alphabet import IUPAC 
>>> dna_seq = Seq("ATGCCGATCGTAT",IUPAC.unambiguous_dna) >>> transcribe(dna_seq) 
Seq('AUGCCGAUCGUAU', IUPACUnambiguousRNA()) 
>>>

要反转转录，T 更改为 U，如下面的代码所示 -

>>> rna_seq = transcribe(dna_seq) 
>>> rna_seq.back_transcribe() 
Seq('ATGCCGATCGTAT', IUPACUnambiguousDNA())

要获得 DNA 模板链，请反向补充反转录的 RNA，如下所示 -

>>> rna_seq.back_transcribe().reverse_complement() 
Seq('ATACGATCGGCAT', IUPACUnambiguousDNA())

翻译

翻译是将RNA序列翻译成蛋白质序列的过程。考虑如下所示的 RNA 序列 -

>>> rna_seq = Seq("AUGGCCAUUGUAAU",IUPAC.unambiguous_rna) 
>>> rna_seq 
Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())

现在，将 translate() 函数应用于上面的代码 -

>>> rna_seq.translate() 
Seq('MAIV', IUPACProtein())

上面的RNA序列很简单。考虑 RNA 序列 AUGGCCAUUGUAAUGGGCCCUGAAAGGGUGCCCGA 并应用 translate() -

>>> rna = Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGA', IUPAC.unambiguous_rna) 
>>> rna.translate() 
Seq('MAIVMGR*KGAR', HasStopCodon(IUPACProtein(), '*'))

这里，终止密码子用星号“*”表示。

在translate()方法中可以在第一个终止密码子处停止。要执行此操作，您可以在 translate() 中指定 to_stop=True ，如下所示 -

>>> rna.translate(to_stop = True) 
Seq('MAIVMGR', IUPACProtein())

这里，终止密码子不包含在结果序列中，因为它不包含终止密码子。

翻译表

NCBI 的遗传密码页面提供了 Biopython 使用的翻译表的完整列表。让我们看一个标准表的示例来可视化代码 -

>>> from Bio.Data import CodonTable 
>>> table = CodonTable.unambiguous_dna_by_name["Standard"] 
>>> print(table) 
Table 1 Standard, SGC0
   | T       | C       | A       | G       | 
 --+---------+---------+---------+---------+-- 
 T | TTT F   | TCT S   | TAT Y   | TGT C   | T
 T | TTC F   | TCC S   | TAC Y   | TGC C   | C
 T | TTA L   | TCA S   | TAA Stop| TGA Stop| A
 T | TTG L(s)| TCG S   | TAG Stop| TGG W   | G 
 --+---------+---------+---------+---------+--
 C | CTT L   | CCT P   | CAT H   | CGT R   | T
 C | CTC L   | CCC P   | CAC H   | CGC R   | C
 C | CTA L   | CCA P   | CAA Q   | CGA R   | A
 C | CTG L(s)| CCG P   | CAG Q   | CGG R   | G 
 --+---------+---------+---------+---------+--
 A | ATT I   | ACT T   | AAT N   | AGT S   | T
 A | ATC I   | ACC T   | AAC N   | AGC S   | C
 A | ATA I   | ACA T   | AAA K   | AGA R   | A
 A | ATG M(s)| ACG T   | AAG K   | AGG R   | G 
 --+---------+---------+---------+---------+--
 G | GTT V   | GCT A   | GAT D   | GGT G   | T
 G | GTC V   | GCC A   | GAC D   | GGC G   | C
 G | GTA V   | GCA A   | GAA E   | GGA G   | A
 G | GTG V   | GCG A   | GAG E   | GGG G   | G 
 --+---------+---------+---------+---------+-- 
>>>

Biopython 使用此表将 DNA 翻译为蛋白质并查找终止密码子。

Biopython - 序列 I/O 操作

Biopython 提供了一个模块 Bio.SeqIO 来分别从文件（任何流）读取序列和向文件（任何流）写入序列。它支持生物信息学中几乎所有可用的文件格式。大多数软件针对不同的文件格式提供不同的方法。但是，Biopython 有意识地遵循单一方法通过其 SeqRecord 对象向用户呈现解析的序列数据。

让我们在下一节中了解有关 SeqRecord 的更多信息。

顺序记录

Bio.SeqRecord 模块提供 SeqRecord 来保存序列的元信息以及序列数据本身，如下所示 -

  • seq - 这是一个实际的序列。

  • id - 它是给定序列的主要标识符。默认类型是字符串。

  • name - 它是序列的名称。默认类型是字符串。

  • 描述- 它显示有关序列的人类可读信息。

  • 注释- 它是有关序列的附加信息的字典。

SeqRecord 可以按如下指定导入

from Bio.SeqRecord import SeqRecord

让我们在接下来的部分中了解使用真实序列文件解析序列文件的细微差别。

解析序列文件格式

本节介绍如何解析两种最流行的序列文件格式：FASTA和GenBank。

FASTA

FASTA是存储序列数据的最基本的文件格式。FASTA最初是生物信息学早期发展过程中开发的DNA和蛋白质序列比对软件包，主要用于搜索序列相似性。

Biopython 提供了一个示例 FASTA 文件，可以通过https://github.com/biopython/biopython/blob/master/Doc/examples/ls_orchid.fasta 访问。

下载此文件并将其保存到您的 Biopython 示例目录中，名称为'orchid.fasta'。

Bio.SeqIO 模块提供了 parse() 方法来处理序列文件，可以按如下方式导入 -

from Bio.SeqIO import parse

parse() 方法包含两个参数，第一个是文件句柄，第二个是文件格式。

>>> file = open('path/to/biopython/sample/orchid.fasta') 
>>> for record in parse(file, "fasta"): 
...    print(record.id) 
... 
gi|2765658|emb|Z78533.1|CIZ78533 
gi|2765657|emb|Z78532.1|CCZ78532 
.......... 
.......... 
gi|2765565|emb|Z78440.1|PPZ78440 
gi|2765564|emb|Z78439.1|PBZ78439 
>>>

这里，parse()方法返回一个可迭代对象，该对象在每次迭代时返回SeqRecord。由于可迭代，它提供了许多复杂且简单的方法，让我们看到一些功能。

下一个（）

next()方法返回可迭代对象中可用的下一个项目，我们可以使用它来获取第一个序列，如下所示 -

>>> first_seq_record = next(SeqIO.parse(open('path/to/biopython/sample/orchid.fasta'),'fasta')) 
>>> first_seq_record.id 'gi|2765658|emb|Z78533.1|CIZ78533' 
>>> first_seq_record.name 'gi|2765658|emb|Z78533.1|CIZ78533' 
>>> first_seq_record.seq Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', SingleLetterAlphabet()) 
>>> first_seq_record.description 'gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA' 
>>> first_seq_record.annotations 
{} 
>>>

这里，seq_record.annotations 为空，因为 FASTA 格式不支持序列注释。

列表理解

我们可以使用列表理解将可迭代对象转换为列表，如下所示

>>> seq_iter = SeqIO.parse(open('path/to/biopython/sample/orchid.fasta'),'fasta') 
>>> all_seq = [seq_record for seq_record in seq_iter] >>> len(all_seq) 
94 
>>>

在这里，我们使用 len 方法来获取总数。我们可以得到最大长度的序列如下 -

>>> seq_iter = SeqIO.parse(open('path/to/biopython/sample/orchid.fasta'),'fasta') 
>>> max_seq = max(len(seq_record.seq) for seq_record in seq_iter) 
>>> max_seq 
789 
>>>

我们也可以使用下面的代码过滤序列 -

>>> seq_iter = SeqIO.parse(open('path/to/biopython/sample/orchid.fasta'),'fasta') 
>>> seq_under_600 = [seq_record for seq_record in seq_iter if len(seq_record.seq) < 600] 
>>> for seq in seq_under_600: 
...    print(seq.id) 
... 
gi|2765606|emb|Z78481.1|PIZ78481 
gi|2765605|emb|Z78480.1|PGZ78480 
gi|2765601|emb|Z78476.1|PGZ78476 
gi|2765595|emb|Z78470.1|PPZ78470 
gi|2765594|emb|Z78469.1|PHZ78469 
gi|2765564|emb|Z78439.1|PBZ78439 
>>>

将 SqlRecord 对象的集合（解析的数据）写入文件就像调用 SeqIO.write 方法一样简单，如下所示 -

file = open("converted.fasta", "w) 
SeqIO.write(seq_record, file, "fasta")

此方法可以有效地用于转换格式，如下所示 -

file = open("converted.gbk", "w) 
SeqIO.write(seq_record, file, "genbank")

基因库

它是一种更丰富的基因序列格式，包含各种注释的字段。Biopython 提供了一个示例 GenBank 文件，可以通过https://github.com/biopython/biopython/blob/master/Doc/examples/ls_orchid.fasta 访问。

下载文件并将其保存到 Biopython 示例目录中，命名为“orchid.gbk”

因为，Biopython 提供了一个单一的函数 parse 来解析所有生物信息学格式。解析 GenBank 格式就像更改解析方法中的格式选项一样简单。

下面给出了相同的代码 -

>>> from Bio import SeqIO 
>>> from Bio.SeqIO import parse 
>>> seq_record = next(parse(open('path/to/biopython/sample/orchid.gbk'),'genbank')) 
>>> seq_record.id 
'Z78533.1' 
>>> seq_record.name 
'Z78533' 
>>> seq_record.seq Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', IUPACAmbiguousDNA()) 
>>> seq_record.description 
'C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA' 
>>> seq_record.annotations {
   'molecule_type': 'DNA', 
   'topology': 'linear', 
   'data_file_division': 'PLN', 
   'date': '30-NOV-2006', 
   'accessions': ['Z78533'], 
   'sequence_version': 1, 
   'gi': '2765658', 
   'keywords': ['5.8S ribosomal RNA', '5.8S rRNA gene', 'internal transcribed spacer', 'ITS1', 'ITS2'], 
   'source': 'Cypripedium irapeanum', 
   'organism': 'Cypripedium irapeanum', 
   'taxonomy': [
      'Eukaryota', 
      'Viridiplantae', 
      'Streptophyta', 
      'Embryophyta', 
      'Tracheophyta', 
      'Spermatophyta', 
      'Magnoliophyta', 
      'Liliopsida', 
      'Asparagales', 
      'Orchidaceae', 
      'Cypripedioideae', 
      'Cypripedium'], 
   'references': [
      Reference(title = 'Phylogenetics of the slipper orchids (Cypripedioideae:
      Orchidaceae): nuclear rDNA ITS sequences', ...), 
      Reference(title = 'Direct Submission', ...)
   ]
}

Biopython - 序列比对

序列比对是以特定顺序排列两个或多个序列（DNA、RNA 或蛋白质序列）以识别它们之间的相似区域的过程。

识别相似区域使我们能够推断出很多信息，例如物种之间保守的性状、不同物种在遗传上的接近程度、物种如何进化等。Biopython 为序列比对提供了广泛的支持。

让我们在本章中学习 Biopython 提供的一些重要功能 -

解析序列对齐

Biopython 提供了一个模块 Bio.AlignIO 来读取和写入序列比对。在生物信息学中，有很多格式可用于指定类似于早期学习的序列数据的序列比对数据。Bio.AlignIO 提供与 Bio.SeqIO 类似的 API，不同之处在于 Bio.SeqIO 作用于序列数据，而 Bio.AlignIO 作用于序列比对数据。

在开始学习之前，让我们从互联网上下载一个示例序列比对文件。

要下载示例文件，请按照以下步骤操作 -

步骤 1 - 打开您最喜欢的浏览器并访问http://pfam.xfam.org/family/browse网站。它将按字母顺序显示所有 Pfam 家族。

步骤 2 - 选择任何一个种子价值较少的家庭。它包含最少的数据，使我们能够轻松地进行对齐工作。在这里，我们选择/单击了 PF18225，它会打开http://pfam.xfam.org/family/PF18225，并显示有关它的完整详细信息，包括序列比对。

步骤 3 - 转到比对部分并下载斯德哥尔摩格式的序列比对文件 (PF18225_seed.txt)。

让我们尝试使用 Bio.AlignIO 读取下载的序列比对文件，如下所示 -

导入Bio.AlignIO模块

>>> from Bio import AlignIO

使用read方法读取对齐。read 方法用于读取给定文件中可用的单一对齐数据。如果给定的文件包含很多对齐方式，我们可以使用解析方法。parse 方法返回可迭代的比对对象，类似于 Bio.SeqIO 模块中的 parse 方法。

>>> alignment = AlignIO.read(open("PF18225_seed.txt"), "stockholm")

打印对齐对象。

>>> print(alignment)
SingleLetterAlphabet() alignment with 6 rows and 65 columns
MQNTPAERLPAIIEKAKSKHDINVWLLDRQGRDLLEQRVPAKVA...EGP B7RZ31_9GAMM/59-123 
AKQRGIAGLEEWLHRLDHSEAIPIFLIDEAGKDLLEREVPADIT...KKP A0A0C3NPG9_9PROT/58-119 
ARRHGQEYFQQWLERQPKKVKEQVFAVDQFGRELLGRPLPEDMA...KKP A0A143HL37_9GAMM/57-121 
TRRHGPESFRFWLERQPVEARDRIYAIDRSGAEILDRPIPRGMA...NKP A0A0X3UC67_9GAMM/57-121 
AINRNTQQLTQDLRAMPNWSLRFVYIVDRNNQDLLKRPLPPGIM...NRK B3PFT7_CELJU/62-126 
AVNATEREFTERIRTLPHWARRNVFVLDSQGFEIFDRELPSPVA...NRT K4KEM7_SIMAS/61-125
>>>

我们还可以检查比对中可用的序列（SeqRecord），如下所示 -

>>> for align in alignment: 
... print(align.seq) 
... 
MQNTPAERLPAIIEKAKSKHDINVWLLDRQGRDLLEQRVPAKVATVANQLRGRKRRAFARHREGP 
AKQRGIAGLEEWLHRLDHSEAIPIFLIDEAGKDLLEREVPADITA---RLDRRREHGEHGVRKKP 
ARRHGQEYFQQWLERQPKKVKEQVFAVDQFGRELLGRPLPEDMAPMLIALNYRNRESHAQVDKKP 
TRRHGPESFRFWLERQPVEARDRIYAIDRSGAEILDRPIPRGMAPLFKVLSFRNREDQGLVNNKP 
AINRNTQQLTQDLRAMPNWSLRFVYIVDRNNQDLLKRPLPPGIMVLAPRLTAKHPYDKVQDRNRK 
AVNATEREFTERIRTLPHWARRNVFVLDSQGFEIFDRELPSPVADLMRKLDLDRPFKKLERKNRT 
>>>

多重对齐

一般来说，大多数序列比对文件都包含单一比对数据，使用read方法来解析它就足够了。在多序列比对概念中，比较两个或多个序列以获得它们之间的最佳子序列匹配，并在单个文件中产生多序列比对。

如果输入序列比对格式包含多个序列比对，那么我们需要使用parse方法而不是read方法，如下所示 -

>>> from Bio import AlignIO 
>>> alignments = AlignIO.parse(open("PF18225_seed.txt"), "stockholm") 
>>> print(alignments) 
<generator object parse at 0x000001CD1C7E0360> 
>>> for alignment in alignments: 
... print(alignment) 
... 
SingleLetterAlphabet() alignment with 6 rows and 65 columns 
MQNTPAERLPAIIEKAKSKHDINVWLLDRQGRDLLEQRVPAKVA...EGP B7RZ31_9GAMM/59-123 
AKQRGIAGLEEWLHRLDHSEAIPIFLIDEAGKDLLEREVPADIT...KKP A0A0C3NPG9_9PROT/58-119 
ARRHGQEYFQQWLERQPKKVKEQVFAVDQFGRELLGRPLPEDMA...KKP A0A143HL37_9GAMM/57-121 
TRRHGPESFRFWLERQPVEARDRIYAIDRSGAEILDRPIPRGMA...NKP A0A0X3UC67_9GAMM/57-121 
AINRNTQQLTQDLRAMPNWSLRFVYIVDRNNQDLLKRPLPPGIM...NRK B3PFT7_CELJU/62-126 
AVNATEREFTERIRTLPHWARRNVFVLDSQGFEIFDRELPSPVA...NRT K4KEM7_SIMAS/61-125
>>>

这里，parse方法返回可迭代的对齐对象，并且可以迭代它以获得实际的对齐。

成对序列比对

成对序列比对一次仅比较两个序列，并提供最佳的序列比对。Pairwise很容易理解，并且可以从结果序列比对中推断出来。

Biopython提供了一个特殊的模块Bio.pairwise2来使用pairwise方法来识别比对序列。Biopython采用最好的算法来查找比对序列，与其他软件不相上下。

让我们编写一个示例，使用pairwise 模块查找两个简单的假设序列的序列比对。这将帮助我们理解序列比对的概念以及如何使用 Biopython 对其进行编程。

步骤1

使用下面给出的命令导入模块pairwise2 -

>>> from Bio import pairwise2

第2步

创建两个序列，seq1 和 seq2 -

>>> from Bio.Seq import Seq 
>>> seq1 = Seq("ACCGGT") 
>>> seq2 = Seq("ACGT")

步骤3

调用方法pairwise2.align.globalxx以及seq1和seq2来使用下面的代码行查找比对 -

>>> alignments = pairwise2.align.globalxx(seq1, seq2)

在这里，globalxx方法执行实际工作并找到给定序列中所有可能的最佳比对。实际上，Bio.pairwise2 提供了相当多的方法，它们遵循以下约定来在不同场景下查找比对。

<sequence alignment type>XY

这里，序列比对类型是指比对类型，可以是全局的，也可以是局部的。全局类型是通过考虑整个序列来查找序列比对。局部类型也是通过查看给定序列的子集来查找序列比对。这将很乏味，但可以更好地了解给定序列之间的相似性。

  • X指的是匹配分数。可能的值为 x（完全匹配）、m（基于相同字符的分数）、d（用户提供的包含字符和匹配分数的字典），最后是 c（用户定义的函数，用于提供自定义评分算法）。

  • Y指间隙罚分。可能的值是 x（无间隙罚分）、s（两个序列的罚分相同）、d（每个序列的罚分不同）以及最后的 c（用户定义的函数，用于提供自定义间隙罚分）

因此，localds 也是一种有效的方法，它使用局部比对技术、用户提供的匹配字典以及用户为两个序列提供的空位罚分来查找序列比对。

>>> test_alignments = pairwise2.align.localds(seq1, seq2, blosum62, -10, -1)

这里，blosum62 指的是pairwise2 模块中可用的字典来提供匹配分数。-10 指空位开放罚分，-1 指空位延伸罚分。

步骤4

循环遍历可迭代的对齐对象并获取每个单独的对齐对象并打印它。

>>> for alignment in alignments: 
... print(alignment) 
... 
('ACCGGT', 'A-C-GT', 4.0, 0, 6) 
('ACCGGT', 'AC--GT', 4.0, 0, 6) 
('ACCGGT', 'A-CG-T', 4.0, 0, 6) 
('ACCGGT', 'AC-G-T', 4.0, 0, 6)

步骤5

Bio.pairwise2 模块提供了一种格式化方法 format_alignment 来更好地可视化结果 -

>>> from Bio.pairwise2 import format_alignment 
>>> alignments = pairwise2.align.globalxx(seq1, seq2) 
>>> for alignment in alignments: 
... print(format_alignment(*alignment)) 
...

ACCGGT 
| | || 
A-C-GT 
   Score=4 
   
ACCGGT 
|| || 
AC--GT 
   Score=4 

ACCGGT 
| || | 
A-CG-T 
   Score=4 

ACCGGT 
|| | | 
AC-G-T 
   Score=4

>>>

Biopython还提供了另一个模块来进行序列比对，Align。该模块提供了一组不同的 API 来简单地设置参数，如算法、模式、匹配分数、差距惩罚等，简单查看 Align 对象如下 -

>>> from Bio import Align
>>> aligner = Align.PairwiseAligner()
>>> print(aligner)
Pairwise sequence aligner with parameters
   match score: 1.000000
   mismatch score: 0.000000
   target open gap score: 0.000000
   target extend gap score: 0.000000
   target left open gap score: 0.000000
   target left extend gap score: 0.000000
   target right open gap score: 0.000000
   target right extend gap score: 0.000000
   query open gap score: 0.000000
   query extend gap score: 0.000000
   query left open gap score: 0.000000
   query left extend gap score: 0.000000
   query right open gap score: 0.000000
   query right extend gap score: 0.000000
   mode: global
>>>

支持序列比对工具

Biopython 通过 Bio.Align.Applications 模块提供了许多序列比对工具的接口。下面列出了一些工具 -

  • 克拉斯塔尔W
  • 肌肉
  • 浮雕针和水

让我们在 Biopython 中编写一个简单的示例，通过最流行的比对工具 ClustalW 创建序列比对。

步骤 1 - 从http://www.clustal.org/download/current/下载 Clustalw 程序并安装。另外，使用“clustal”安装路径更新系统路径。

步骤 2 - 从模块 Bio.Align.Applications 导入 ClustalwCommanLine。

>>> from Bio.Align.Applications import ClustalwCommandline

步骤 3 - 通过使用 Biopython 包中提供的输入文件 opuntia.fasta 调用 ClustalwCommanLine 来设置 cmd。 https://raw.githubusercontent.com/biopython/biopython/master/Doc/examples/opuntia.fasta

>>> cmd = ClustalwCommandline("clustalw2",
infile="/path/to/biopython/sample/opuntia.fasta")
>>> print(cmd)
clustalw2 -infile=fasta/opuntia.fasta

步骤 4 - 调用 cmd() 将运行 clustalw 命令并给出结果对齐文件 opuntia.aln 的输出。

>>> stdout, stderr = cmd()

步骤 5 - 读取并打印对齐文件，如下 -

>>> from Bio import AlignIO
>>> align = AlignIO.read("/path/to/biopython/sample/opuntia.aln", "clustal")
>>> print(align)
SingleLetterAlphabet() alignment with 7 rows and 906 columns
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273285|gb|AF191659.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273284|gb|AF191658.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273287|gb|AF191661.1|AF191
TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273286|gb|AF191660.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273290|gb|AF191664.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273289|gb|AF191663.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273291|gb|AF191665.1|AF191
>>>

Biopython - BLAST 概述

BLAST 代表基本局部对齐搜索工具。它发现生物序列之间的相似区域。Biopython提供Bio.Blast模块来处理NCBI BLAST操作。您可以在本地连接或通过 Internet 连接运行 BLAST。

让我们在下一节中简要了解这两种联系 -

通过互联网运行

Biopython提供了Bio.Blast.NCBIWWW模块来调用在线版本的BLAST。为此，我们需要导入以下模块 -

>>> from Bio.Blast import NCBIWWW

NCBIWW模块提供了qblast功能来查询BLAST在线版本，https://blast.ncbi.nlm.nih.gov/Blast.cgi。qblast 支持在线版本支持的所有参数。

要获得有关此模块的任何帮助，请使用以下命令并了解其功能 -

>>> help(NCBIWWW.qblast) 
Help on function qblast in module Bio.Blast.NCBIWWW: 
qblast(
   program, database, sequence, 
   url_base = 'https://blast.ncbi.nlm.nih.gov/Blast.cgi', 
   auto_format = None, 
   composition_based_statistics = None, 
   db_genetic_code =  None, 
   endpoints = None, 
   entrez_query = '(none)', 
   expect = 10.0, 
   filter = None, 
   gapcosts = None, 
   genetic_code = None, 
   hitlist_size = 50, 
   i_thresh = None, 
   layout = None, 
   lcase_mask = None, 
   matrix_name = None, 
   nucl_penalty = None, 
   nucl_reward = None, 
   other_advanced = None, 
   perc_ident = None, 
   phi_pattern = None, 
   query_file = None, 
   query_believe_defline = None, 
   query_from = None, 
   query_to = None, 
   searchsp_eff = None, 
   service = None, 
   threshold = None, 
   ungapped_alignment = None, 
   word_size = None, 
   alignments = 500, 
   alignment_view = None, 
   descriptions = 500, 
   entrez_links_new_window = None, 
   expect_low = None, 
   expect_high = None, 
   format_entrez_query = None, 
   format_object = None, 
   format_type = 'XML', 
   ncbi_gi = None, 
   results_file = None, 
   show_overview = None, 
   megablast = None, 
   template_type = None, 
   template_length = None
) 
   
   BLAST search using NCBI's QBLAST server or a cloud service provider. 
   
   Supports all parameters of the qblast API for Put and Get. 
   
   Please note that BLAST on the cloud supports the NCBI-BLAST Common 
   URL API (http://ncbi.github.io/blast-cloud/dev/api.html). 
   To use this feature, please set url_base to 'http://host.my.cloud.service.provider.com/cgi-bin/blast.cgi' and 
   format_object = 'Alignment'. For more details, please see 8. Biopython – Overview of BLAST
   
https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE = BlastDocs&DOC_TYPE = CloudBlast 
   
   Some useful parameters: 
   
   - program blastn, blastp, blastx, tblastn, or tblastx (lower case) 
   - database Which database to search against (e.g. "nr"). 
   - sequence The sequence to search. 
   - ncbi_gi TRUE/FALSE whether to give 'gi' identifier. 
   - descriptions Number of descriptions to show. Def 500. 
   - alignments Number of alignments to show. Def 500. 
   - expect An expect value cutoff. Def 10.0. 
   - matrix_name Specify an alt. matrix (PAM30, PAM70, BLOSUM80, BLOSUM45). 
   - filter "none" turns off filtering. Default no filtering 
   - format_type "HTML", "Text", "ASN.1", or "XML". Def. "XML". 
   - entrez_query Entrez query to limit Blast search 
   - hitlist_size Number of hits to return. Default 50 
   - megablast TRUE/FALSE whether to use MEga BLAST algorithm (blastn only) 
   - service plain, psi, phi, rpsblast, megablast (lower case) 
   
   This function does no checking of the validity of the parameters 
   and passes the values to the server as is. More help is available at: 
   https://ncbi.github.io/blast-cloud/dev/api.html

通常，qblast 函数的参数基本上类似于您可以在 BLAST 网页上设置的不同参数。这使得 qblast 函数易于理解，并减少了使用它的学习曲线。

连接和搜索

为了了解连接和搜索 BLAST 在线版本的过程，让我们通过 Biopython 对在线 BLAST 服务器进行简单的序列搜索（在我们的本地序列文件中可用）。

步骤 1 - 在 Biopython 目录中创建一个名为blast_example.fasta的文件，并提供以下序列信息作为输入

Example of a single sequence in FASTA/Pearson format: 
>sequence A ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc 

>sequence B ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattca
tattctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc

步骤 2 - 导入 NCBIWWW 模块。

>>> from Bio.Blast import NCBIWWW

步骤 3 - 使用 python IO 模块打开序列文件blast_example.fasta。

>>> sequence_data = open("blast_example.fasta").read() 
>>> sequence_data 
'Example of a single sequence in FASTA/Pearson format:\n\n\n> sequence 
A\nggtaagtcctctagtacaaacacccccaatattgtgatataattaaaatt 
atattcatat\ntctgttgccagaaaaaacacttttaggctatattagagccatcttctttg aagcgttgtc\n\n'

步骤 4 - 现在，调用 qblast 函数，传递序列数据作为主要参数。另一个参数代表数据库（nt）和内部程序（blastn）。

>>> result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data) 
>>> result_handle 
<_io.StringIO object at 0x000001EC9FAA4558>

blast_results保存我们的搜索结果。它可以保存到文件中以供以后使用，也可以进行解析以获取详细信息。我们将在下一节中学习如何做到这一点。

步骤 5 - 也可以使用 Seq 对象来完成相同的功能，而不是使用整个 fasta 文件，如下所示 -

>>> from Bio import SeqIO 
>>> seq_record = next(SeqIO.parse(open('blast_example.fasta'),'fasta')) 
>>> seq_record.id 
'sequence' 
>>> seq_record.seq 
Seq('ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatat...gtc', 
SingleLetterAlphabet())

现在，调用 qblast 函数，传递 Seq 对象 record.seq 作为主要参数。

>>> result_handle = NCBIWWW.qblast("blastn", "nt", seq_record.seq) 
>>> print(result_handle) 
<_io.StringIO object at 0x000001EC9FAA4558>

BLAST 将自动为您的序列分配一个标识符。

步骤 6 - result_handle 对象将拥有完整的结果，并且可以保存到文件中以供以后使用。

>>> with open('results.xml', 'w') as save_file: 
>>>   blast_results = result_handle.read() 
>>>   save_file.write(blast_results)

我们将在后面的部分中看到如何解析结果文件。

运行独立 BLAST

本节介绍如何在本地系统中运行 BLAST。如果您在本地系统中运行 BLAST，它可能会更快，并且还允许您创建自己的数据库来搜索序列。

连接 BLAST

一般来说，不建议在本地运行 BLAST，因为它的体积较大，需要额外的工作来运行软件，并且涉及成本。Online BLAST 足以满足基本和高级目的。当然，有时可能需要您在本地安装。

考虑到您经常进行在线搜索，这可能需要大量时间和高网络量，并且如果您有专有序列数据或 IP 相关问题，则建议在本地安装。

为此，我们需要执行以下步骤 -

步骤 1 - 使用给定链接下载并安装最新的blast二进制文件 - ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

步骤 2 - 使用以下链接下载并解压最新且必要的数据库 - ftp://ftp.ncbi.nlm.nih.gov/blast/db/

BLAST 软件在其网站上提供了大量数据库。让我们从blast数据库站点下载alu.n.gz文件并将其解压到alu文件夹中。该文件为 FASTA 格式。要在我们的blast应用程序中使用此文件，我们需要首先将该文件从FASTA格式转换为blast数据库格式。BLAST 提供了 makeblastdb 应用程序来执行此转换。

使用下面的代码片段 -

cd /path/to/alu 
makeblastdb -in alu.n -parse_seqids -dbtype nucl -out alun

运行上面的代码将解析输入文件 alu.n 并创建 BLAST 数据库作为多个文件 alun.nsq、alun.nsi 等。现在，我们可以查询该数据库以查找序列。

我们已经在本地服务器中安装了 BLAST，并且还有示例 BLAST 数据库alun来对其进行查询。

步骤 3 - 让我们创建一个示例序列文件来查询数据库。创建文件 search.fsa 并将以下数据放入其中。

>gnl|alu|Z15030_HSAL001056 (Alu-J) 
AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCT 
TGAGCCTAGGAGTTTGCGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAA 
AGAAAAAAAAAATAGCTCTGCTGGTGGTGCATGCCTATAGTCTCAGCTATCAGGAGGCTG 
GGACAGGAGGATCACTTGGGCCCGGGAGTTGAGGCTGTGGTGAGCCACGATCACACCACT 
GCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAACAAATAA 
>gnl|alu|D00596_HSAL003180 (Alu-Sx) 
AGCCAGGTGTGGTGGCTCACGCCTGTAATCCCACCGCTTTGGGAGGCTGAGTCAGATCAC 
CTGAGGTTAGGAATTTGGGACCAGCCTGGCCAACATGGCGACACCCCAGTCTCTACTAAT 
AACACAAAAAATTAGCCAGGTGTGCTGGTGCATGTCTGTAATCCCAGCTACTCAGGAGGC 
TGAGGCATGAGAATTGCTCACGAGGCGGAGGTTGTAGTGAGCTGAGATCGTGGCACTGTA
CTCCAGCCTGGCGACAGAGGGAGAACCCATGTCAAAAACAAAAAAAGACACCACCAAAGG 
TCAAAGCATA 
>gnl|alu|X55502_HSAL000745 (Alu-J) 
TGCCTTCCCCATCTGTAATTCTGGCACTTGGGGAGTCCAAGGCAGGATGATCACTTATGC 
CCAAGGAATTTGAGTACCAAGCCTGGGCAATATAACAAGGCCCTGTTTCTACAAAAACTT 
TAAACAATTAGCCAGGTGTGGTGGTGCGTGCCTGTGTCCAGCTACTCAGGAAGCTGAGGC 
AAGAGCTTGAGGCTACAGTGAGCTGTGTTCCACCATGGTGCTCCAGCCTGGGTGACAGGG 
CAAGACCCTGTCAAAAGAAAGGAAGAAAGAACGGAAGGAAAGAAGGAAAGAAACAAGGAG 
AG

序列数据从 alu.n 文件收集；因此，它与我们的数据库匹配。

步骤 4 - BLAST 软件提供了许多应用程序来搜索数据库，我们使用blastn。blastn 应用程序至少需要三个参数：db、query 和 out。db指要搜索的数据库；query是要匹配的序列，out是存储结果的文件。现在，运行以下命令来执行这个简单的查询 -

blastn -db alun -query search.fsa -out results.xml -outfmt 5

运行上述命令将在results.xml文件中搜索并给出输出，如下所示（部分数据） -

<?xml version = "1.0"?> 
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" 
   "http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">
<BlastOutput> 
   <BlastOutput_program>blastn</BlastOutput_program> 
   <BlastOutput_version>BLASTN 2.7.1+</BlastOutput_version> 
   <BlastOutput_reference>Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb 
      Miller (2000), "A greedy algorithm for aligning DNA sequences", J 
      Comput Biol 2000; 7(1-2):203-14.
   </BlastOutput_reference> 
   
   <BlastOutput_db>alun</BlastOutput_db> 
   <BlastOutput_query-ID>Query_1</BlastOutput_query-ID> 
   <BlastOutput_query-def>gnl|alu|Z15030_HSAL001056 (Alu-J)</BlastOutput_query-def>
   <BlastOutput_query-len>292</BlastOutput_query-len> 
   <BlastOutput_param> 
      <Parameters> 
         <Parameters_expect>10</Parameters_expect> 
         <Parameters_sc-match>1</Parameters_sc-match> 
         <Parameters_sc-mismatch>-2</Parameters_sc-mismatch> 
         <Parameters_gap-open>0</Parameters_gap-open> 
         <Parameters_gap-extend>0</Parameters_gap-extend> 
         <Parameters_filter>L;m;</Parameters_filter> 
      </Parameters> 
   </BlastOutput_param> 
   <BlastOutput_iterations> 
      <Iteration> 
         <Iteration_iter-num>1</Iteration_iter-num><Iteration_query-ID>Query_1</Iteration_query-ID> 
         <Iteration_query-def>gnl|alu|Z15030_HSAL001056 (Alu-J)</Iteration_query-def> 
         <Iteration_query-len>292</Iteration_query-len> 
         <Iteration_hits> 
         <Hit>
            <Hit_num>1</Hit_num> 
            <Hit_id>gnl|alu|Z15030_HSAL001056</Hit_id> 
            <Hit_def>(Alu-J)</Hit_def> 
            <Hit_accession>Z15030_HSAL001056</Hit_accession> 
            <Hit_len>292</Hit_len>
            <Hit_hsps> 
               <Hsp>
                 <Hsp_num>1</Hsp_num> 
                  <Hsp_bit-score>540.342</Hsp_bit-score> 
                  <Hsp_score>292</Hsp_score>
                  <Hsp_evalue>4.55414e-156</Hsp_evalue> 
                  <Hsp_query-from>1</Hsp_query-from>
                  <Hsp_query-to>292</Hsp_query-to> 
                  <Hsp_hit-from>1</Hsp_hit-from> 
                  <Hsp_hit-to>292</Hsp_hit-to> 
                  <Hsp_query-frame>1</Hsp_query-frame>
                  <Hsp_hit-frame>1</Hsp_hit-frame>
                  <Hsp_identity>292</Hsp_identity>
                  <Hsp_positive>292</Hsp_positive> 
                  <Hsp_gaps>0</Hsp_gaps> 
                  <Hsp_align-len>292</Hsp_align-len>
                  
                  <Hsp_qseq>
                     AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCTTGAGCCTAGGAGTTTG
                     CGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAAAGAAAAAAAAAATAGCTCTGCTGGTGGTGCATG
                     CCTATAGTCTCAGCTATCAGGAGGCTGGGACAGGAGGATCACTTGGGCCCGGGAGTTGAGGCTGTGGTGAGCC
                     ACGATCACACCACTGCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAACAAATAA
                  </Hsp_qseq> 

                  <Hsp_hseq>
                     AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCTTGAGCCTAGGA
                     GTTTGCGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAAAGAAAAAAAAAATAGCTCTGCT
                     GGTGGTGCATGCCTATAGTCTCAGCTATCAGGAGGCTGGGACAGGAGGATCACTTGGGCCCGGGAGTTGAGG
                     CTGTGGTGAGCCACGATCACACCACTGCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAAC
                     AAATAA
                  </Hsp_hseq>

                  <Hsp_midline>
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||
                  </Hsp_midline>
               </Hsp> 
            </Hit_hsps>
         </Hit>
         ......................... 
         ......................... 
         ......................... 
         </Iteration_hits> 
         <Iteration_stat> 
            <Statistics> 
               <Statistics_db-num>327</Statistics_db-num> 
               <Statistics_db-len>80506</Statistics_db-len> 
               <Statistics_hsp-lenv16</Statistics_hsp-len> 
               <Statistics_eff-space>21528364</Statistics_eff-space> 
               <Statistics_kappa>0.46</Statistics_kappa> 
               <Statistics_lambda>1.28</Statistics_lambda> 
               <St