Biopython 简明教程
Biopython - Introduction
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 等姊妹项目。
Features
Biopython 可移植、清晰且具有易于学习的语法。下面列出了一些显著的特性：
-
解释型、交互式且面向对象。
-
支持 FASTA、PDB、GenBank、Blast、SCOP、PubMed/Medline、ExPASy 相关格式。
-
处理序列格式的选项。
-
管理蛋白质结构的工具。
-
BioSQL - 用于存储序列及特性和注释的标准 SQL 表集。
-
访问在线服务和数据库，包括 NCBI 服务（Blast、Entrez、PubMed）和 ExPASY 服务（SwissProt、Prosite）。
-
访问本地服务，包括 Blast、Clustalw、EMBOSS。
Goals
Biopython 的目标是通过 Python 语言为生物信息学提供简单、标准且广泛的访问途径。以下列出 Biopython 的具体目标：
-
提供对生物信息学资源的标准化访问。
-
高质量、可重复使用的模块和脚本。
-
可在集群代码、PDB、朴素贝叶斯和马尔可夫模型中使用的快速数组操纵。
-
Genomic data analysis.
Advantages
Biopython 所需的代码极少，并具有以下优势：
-
提供用于聚类中的微阵列数据类型。
-
读取并写入 Tree-View 类型文件。
-
支持用于 PDB 解析、表示和分析的结构数据。
-
支持在 Medline 应用程序中使用的日志数据。
-
支持 BioSQL 数据库，这是所有生物信息学项目中广泛使用的标准数据库。
-
通过提供用于将生物信息学文件解析为特定记录对象格式或序列加上特性的通用类别的模块来支持解析器开发。
-
基于食谱风格的清晰文档。
Biopython - Installation
本节说明如何在机器上安装 Biopython。安装它非常容易，不超过五分钟。
Step 1 − 验证 Python 安装
Biopython 被设计为与 Python 2.5 或更高版本一起使用。因此，必须先安装 python。在命令提示符中运行以下命令 −
> python --version它在下面定义 −
如果安装正确，它将显示 python 的版本。否则，下载最新版本的 python，将其安装，然后重新运行该命令。
Step 2 − 使用 pip 安装 Biopython
在所有平台上，使用 pip 从命令行安装 Biopython 非常简单。键入以下命令 −
> pip install biopython将在屏幕上看到以下响应 −
要更新较早版本的 Biopython −
> pip install biopython –-upgrade将在屏幕上看到以下响应 −
执行此命令后，在安装最新版本的 Biopython 和 NumPy（Biopython 依赖于它）之前，将删除旧版本的 Biopython 和 NumPy。
Step 3 − 验证 Biopython 安装
现在, 你已经成功地在你的机器上安装了 Biopython。为了验证 Biopython 是否正确地安装, 在你的 python 控制台上输入一下命令 −
它显示 Biopython 的版本。
Alternate Way − Installing Biopython using Source
要使用源代码安装 Biopython, 请遵循下面的说明 −
从以下链接下载最新版的 Biopython − https://biopython.org/wiki/Download
到目前为止, 最新版本是 biopython-1.72 。
下载文件并解压压缩存档文件, 移到源代码文件夹并输入以下命令 −
> python setup.py build这将从源代码构建 Biopython, 如下所示 −
现在, 使用以下命令测试代码 −
> python setup.py test
最后, 使用以下命令安装 −
> python setup.py install
Biopython - Creating Simple Application
让我们创建一个简单的 Biopython 应用程序来解析生物信息学文件并打印其内容。这将帮助我们理解 Biopython 的一般概念以及它如何在生物信息学领域发挥作用。
Step 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、name、description 和实际序列数据。
Step 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)让我们对代码进行更深入的了解 -
Line 1 导入 Bio.SeqIO 模块中可用的 parse 类。Bio.SeqIO 模块用于读写不同格式的序列文件，`parse' 类用于解析序列文件的内容。
Line 2 导入 Bio.SeqRecord 模块中可用的 SeqRecord 类。此模块用于处理序列记录，SeqRecord 类用于表示序列文件中可用的特定序列。
*Line 3" 导入 Seq 类，该类在 Bio.Seq 模块中可用。此模块用于处理序列数据，Seq 类用于表示序列文件中可用特定序列记录的序列数据。
Line 5 使用常规 python 函数 open 打开 "example.fasta" 文件。
Line 7 解析序列文件的内容，并以 SeqRecord 对象的列表形式返回内容。
Line 9-15 使用 python for 循环遍历记录，并打印序列记录（SqlRecord）的属性，例如 id、name、description、sequence data 等。
Line 15 使用 Alphabet 类打印序列的类型。
Step 3 - 打开命令提示符并转到包含序列文件 “example.fasta” 的文件夹并运行以下命令 -
> python simple_example.pyStep 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 - Sequence
一个序列是一系列用于表示生物体的蛋白质、DNA 或 RNA 的字母。由 Seq 类表示。Seq 类在 Bio.Seq 模块中定义。
让我们像下面那样在 Biopython 中创建一个简单的序列 −
>>> from Bio.Seq import Seq
>>> seq = Seq("AGCT")
>>> seq
Seq('AGCT')
>>> print(seq)
AGCT在这里，我们创建了一个简单的蛋白质序列 AGCT ，每个字母分别表示*A*兰ine、*G*lycine、*C*ysteine 和 *T*hreonine。
每个 Seq 对象都有两个重要的属性 −
-
data − 实际的序列字符串 (AGCT)
-
alphabet − 用于表示序列类型。例如，DNA 序列、RNA 序列等。默认情况下，它不表示任何序列，并且本质上是通用的。
Alphabet Module
Seq 对象包含 Alphabet 属性以指定序列类型、字母和可能的操作。它在 Bio.Alphabet 模块中定义。可以如下定义字母表 −
>>> from Bio.Seq import Seq
>>> myseq = Seq("AGCT")
>>> myseq
Seq('AGCT')
>>> myseq.alphabet
Alphabet()Alphabet 模块提供了以下类来表示不同类型的序列。Alphabet - 所有类型字母表的基本类。
SingleLetterAlphabet - 通用字母表，字母大小为一。它源自 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 (protein) − IUPAC 20 个标准氨基酸蛋白质字母表。
-
ExtendedIUPACProtein (extended_protein) − 扩展的大写 IUPAC 蛋白质单字母字母表，包括 X。
-
IUPACAmbiguousDNA (ambiguous_dna) − 大写 IUPAC 含混 DNA。
-
IUPACUnambiguousDNA (unambiguous_dna) − 大写 IUPAC 明确 DNA (GATC)。
-
ExtendedIUPACDNA (extended_dna) − 扩展 IUPAC DNA 字母表。
-
IUPACAmbiguousRNA (ambiguous_rna) − 大写 IUPAC 含混 RNA。
-
IUPACUnambiguousRNA (unambiguous_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'Basic Operations
本节简要地阐述了 Seq 类中可用的所有基本操作。序列类似于 python 字符串。我们可以在序列中执行 python 字符串操作，如切片、计数、连接、查找、拆分和去除。
使用以下代码可获得各种输出。
To get the first value in sequence.
>>> seq_string = Seq("AGCTAGCT")
>>> seq_string[0]
'A'To print the first two values.
>>> seq_string[0:2]
Seq('AG')To print all the values.
>>> seq_string[ : ]
Seq('AGCTAGCT')To perform length and count operations.
>>> len(seq_string)
8
>>> seq_string.count('A')
2To add two sequences.
>>> 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())在下文中给出了各种代码以根据要求获得输出。
To change the case of sequence.
>>> from Bio.Alphabet import generic_rna
>>> rna = Seq("agct", generic_rna)
>>> rna.upper()
Seq('AGCT', RNAAlphabet())To check python membership and identity operator.
>>> rna = Seq("agct", generic_rna)
>>> 'a' in rna
True
>>> 'A' in rna
False
>>> rna1 = Seq("AGCT", generic_dna)
>>> rna is rna1
FalseTo find single letter or sequence of letter inside the given sequence.
>>> protein_seq = Seq('AGUACACUGGU', generic_protein)
>>> protein_seq.find('G')
1
>>> protein_seq.find('GG')
8To perform splitting operation.
>>> protein_seq = Seq('AGUACACUGGU', generic_protein)
>>> protein_seq.split('A')
[Seq('', ProteinAlphabet()), Seq('GU', ProteinAlphabet()),
Seq('C', ProteinAlphabet()), Seq('CUGGU', ProteinAlphabet())]To perform strip operations in the sequence.
>>> strip_seq = Seq(" AGCT ")
>>> strip_seq
Seq(' AGCT ')
>>> strip_seq.strip()
Seq('AGCT')Biopython - Advanced Sequence Operations
在本章中，我们将讨论 Biopython 提供的某些高级序列特征。
Complement and Reverse Complement
可以对核苷酸序列进行反向互补以获取新序列。此外，互补序列可以进行反向互补以获取原始序列。Biopython 提供了两种方法来执行此功能 − complement 和 reverse_complement 。代码如下 −
>>> from Bio.Alphabet import IUPAC
>>> nucleotide = Seq('TCGAAGTCAGTC', IUPAC.ambiguous_dna)
>>> nucleotide.complement()
Seq('AGCTTCAGTCAG', IUPACAmbiguousDNA())
>>>在此，complement() 方法允许对 DNA 或 RNA 序列进行互补。reverse_complement() 方法对结果序列进行互补，并从左到右对结果序列进行反转。如下所示 −
>>> nucleotide.reverse_complement()
Seq('GACTGACTTCGA', IUPACAmbiguousDNA())Biopython 使用 Bio.Data.IUPACData 提供的 ambiguous_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 Content
预测基因组 DNA 碱基组分（GC 含量）会显著影响基因组功能和物种生态。GC 含量是 GC 核苷酸数除以总核苷酸数。
要获取 GC 核苷酸含量，请导入以下模块并执行以下步骤：
>>> from Bio.SeqUtils import GC
>>> nucleotide = Seq("GACTGACTTCGA",IUPAC.unambiguous_dna)
>>> GC(nucleotide)
50.0Transcription
转录是将 DNA 序列更改为 RNA 序列的过程。实际的生物转录过程是执行反向互补（TCAG → CUGA）以将 DNA 视为模板链来获取 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())Translation
翻译是将 RNA 序列翻译成肽序列的过程。考虑如下所示的 RNA 序列 -
>>> rna_seq = Seq("AUGGCCAUUGUAAU",IUPAC.unambiguous_rna)
>>> rna_seq
Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())现在，对上述代码应用 translate() 函数 -
>>> rna_seq.translate()
Seq('MAIV', IUPACProtein())上述 RNA 序列很简单。考虑 RNA 序列 AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGA 并应用 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())此处，终止密码子不包含在结果序列中，因为它不包含终止密码子。
Translation Table
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 - Sequence I/O Operations
Biopython 提供了一个模块 Bio.SeqIO, 分别用于从文件（任何流）读取和写入序列。它支持生物信息学中几乎所有可用的文件格式。大多数软件为不同的文件格式提供了不同的方法。但是, Biopython 认真地遵循一种方法, 通过其 SeqRecord 对象向用户显示解析的序列数据。
让我们在以下部分详细了解 SeqRecord。
SeqRecord
Bio.SeqRecord 模块提供 SeqRecord 来保存序列的元信息以及序列数据本身, 如下所示 −
-
seq − 这是实际序列。
-
id − 这是给定序列的主标识符。默认类型为字符串。
-
name − 这是序列的名称。默认类型为字符串。
-
description − 它显示关于该序列的人类可读信息。
-
annotations − 这是一个关于该序列的附加信息的字典。
以下是可以导入 SeqRecord 的指定方式
from Bio.SeqRecord import SeqRecord让我们逐步使用真实的序列文件理解解析序列文件时的微妙差别。
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 parseparse() 方法包含两个参数，第一个是文件句柄，第二个是文件格式。
>>> 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()
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 格式不支持序列注释。
list comprehension
我们可以使用列表解析将可迭代对象转换为列表，如下所示
>>> 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")GenBank
这是一种更为丰富的基因序列格式，并且包括各种注释类的字段。Biopython 提供了一个 GenBank 示例文件，该文件可通过 https://github.com/biopython/biopython/blob/master/Doc/examples/ls_orchid.fasta. 访问
将该文件下载并保存在您的 Biopython 示例目录中，作为 ‘orchid.gbk’ 。
由于 Biopython 提供了一个函数来解析所有生物信息学格式，所以解析 GenBank 格式只需在解析方法中改变 format 选项即可。
以下是相对应的代码：
>>> 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 - Sequence Alignments
Sequence alignment 是按特定顺序布置两个或多个序列（DNA、RNA 或肽序列）的过程，以识别它们之间的相似区域。
识别相似区域使我们能够推断大量信息，例如物种之间有哪些性状得到了保守，不同物种在遗传上有多接近，物种如何进化等。Biopython 为序列比对提供了广泛的支持。
让我们在本章中学习 Biopython 提供的一些重要功能 -
Parsing Sequence Alignment
Biopython 提供了一个模块 Bio.AlignIO 来读写序列比对。在生物信息学中，有很多可用的格式可用于指定序列比对数据，类似于前面学到的序列数据。Bio.AlignIO 提供类似于 Bio.SeqIO 的 API，不同之处在于 Bio.SeqIO 用作序列数据，而 Bio.AlignIO 用作序列比对数据。
在开始学习之前，我们从互联网下载一个序列比对示例文件。
若要下载示例文件，请执行以下步骤 -
Step 1 - 打开您最喜欢的浏览器并转到 http://pfam.xfam.org/family/browse 网站。它会按字母顺序列出所有 Pfam 家族。
Step 2 - 选择任何种子值较少的家族。它包含最少的数据，使我们能够轻松使用比对。在此，我们选择了/点击了 PF18225，然后它打开转到 http://pfam.xfam.org/family/PF18225 并显示有关它的完整详细信息，包括序列比对。
Step 3 - 转到比对部分并下载斯德哥尔摩格式的序列比对文件 (PF18225_seed.txt)。
让我们尝试使用 Bio.AlignIO 读取下载的序列比对文件，如下所示 -
Import Bio.AlignIO module
>>> from Bio import AlignIO使用read方法进行读取比对。read方法用于读取给定文件中可用的单个比对数据。如果给定文件包含多个比对，我们可以使用parse方法。parse方法返回可迭代比对对象，类似于Bio.SeqIO模块中的parse方法。
>>> alignment = AlignIO.read(open("PF18225_seed.txt"), "stockholm")Print the alignment object.
>>> 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
>>>Multiple Alignments
通常情况下，大多数序列比对文件包含单个比对数据，只要使用 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 Sequence Alignment
Pairwise sequence alignment 一次只比较两个序列，并提供尽可能好的序列比对。 Pairwise 易于理解，并且可以从得到的序列比对中推断出例外。
Biopython提供了一个特殊模块 Bio.pairwise2 ，以使用成对方法识别比对序列。Biopython应用最佳算法找到比对序列，并且它与其他软件相当。
让我们通过使用成对模块编写一个示例来查找两个简单假设序列的序列比对。这将帮助我们理解序列比对的概念以及如何使用Biopython对其进行编程。
Step 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是指缺口扩展罚分。
Step 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)Step 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
>>>Support for Sequence Alignment Tools
Biopython 通过 Bio.Align.Applications 模块提供大量序列比对工具的接口。下面列出了一些工具 -
-
ClustalW
-
MUSCLE
-
EMBOSS needle and water
让我们编写一个简单的 Biopython 示例，通过最流行的对齐工具 ClustalW 创建序列对齐。
Step 1 − 从 http://www.clustal.org/download/current/ 下载 Clustalw 程序并将其安装。此外，使用“clustal”安装路径更新系统 PATH。
Step 2 − 从 Bio.Align.Applications 模块导入 ClustalwCommanLine。
>>> from Bio.Align.Applications import ClustalwCommandlineStep 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.fastaStep 4 − 调用 cmd() 将运行 clustalw 命令，并给出结果对齐文件 opuntia.aln 的输出。
>>> stdout, stderr = cmd()Step 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 - Overview of BLAST
BLAST 代表 Basic Local Alignment Search Tool 。它在生物序列之间找到相似区域。Biopython 提供了 Bio.Blast 模块来处理 NCBI BLAST 操作。你可以在本地连接或互联网连接上运行 BLAST。
让我们在以下部分简要了解这两个连接 -
Running over Internet
Biopython 提供 Bio.Blast.NCBIWWW 模块来调用 BLAST 的在线版本。要做到这一点，我们需要导入以下模块 -
>>> from Bio.Blast import NCBIWWWNCBIWW 模块提供 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 函数易于理解，并且减少了使用它的学习曲线。
Connecting and Searching
为了理解连接和搜索 BLAST 在线版本的过程，让我们通过 Biopython 对本地序列文件中的简单序列搜索进行在线 BLAST 服务器搜索。
Step 1 − 在 Biopython 目录中创建一个名为 blast_example.fasta 的文件，并输入以下序列信息作为输入
Example of a single sequence in FASTA/Pearson format:
>sequence A ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc
>sequence B ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattca
tattctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtcStep 2 − 导入 NCBIWWW 模块。
>>> from Bio.Blast import NCBIWWWStep 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'Step 4 − 现在，调用 qblast 函数，将序列数据作为主要参数传递。其他参数表示数据库 (nt) 和内部程序 (blastn)。
>>> result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data)
>>> result_handle
<_io.StringIO object at 0x000001EC9FAA4558>blast_results 保存了我们的搜索结果。可以将其存储到文件中以便稍后使用，还可以进行解析以获取详细信息。我们将在下一部分学习如何执行此操作。
Step 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 将自动为您的序列分配一个标识符。
Step 6 − result_handle 对象将拥有整个结果，可以将其保存到文件中以便稍后使用。
>>> with open('results.xml', 'w') as save_file:
>>> blast_results = result_handle.read()
>>> save_file.write(blast_results)我们将在后面的部分中了解如何解析结果文件。
Running Standalone BLAST
此部分说明如何在本地系统上运行 BLAST。如果你在本地系统上运行 BLAST，可能会更快，还可以让你创建自己的数据库来搜索序列。
Connecting BLAST
通常，不建议在本地运行 BLAST，因为它体积过大，需要额外的精力来运行软件，并且涉及成本。在线 BLAST 足以满足基本和高级目的。当然，有时候你可能需要在本地安装它。
考虑到你正在进行频繁的在线搜索，这可能需要大量时间和高网络容量，如果你有专有序列数据或 IP 相关问题，那么建议在本地安装它。
要做到这一点，我们需要执行以下步骤 −
Step 1 − 使用给定的链接下载并安装最新的 blast 二进制文件 − ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
Step 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 来对照其进行查询。
Step 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 文件收集的；因此，它与我们的数据库匹配。
Step 4 − BLAST 软件提供了许多应用程序来搜索数据库，我们使用 blastn。 blastn application requires minimum of three arguments, db, query and 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>
<Statistics_entropy>0.85</Statistics_entropy>
</Statistics>
</Iteration_stat>
</Iteration>
</BlastOutput_iterations>
</BlastOutput>可以使用以下代码在 Python 内运行上述命令：
>>> from Bio.Blast.Applications import NcbiblastnCommandline
>>> blastn_cline = NcbiblastnCommandline(query = "search.fasta", db = "alun",
outfmt = 5, out = "results.xml")
>>> stdout, stderr = blastn_cline()其中，第一个是 BLAST 输出的处理程序，而第二个是 BLAST 命令生成的可能的错误输出。
由于我们已提供输出文件作为命令行参数（out = “results.xml”）并将输出格式设置为 XML（outfmt = 5），因此输出文件将保存于当前的工作目录。
Parsing BLAST Result
通常，BLAST 输出使用 NCBIXML 模块解析为 XML 格式。为此，我们需要导入以下模块：
>>> from Bio.Blast import NCBIXML现在， open the file directly using python open method 和 use NCBIXML parse method 如下所示：
>>> E_VALUE_THRESH = 1e-20
>>> for record in NCBIXML.parse(open("results.xml")):
>>> if record.alignments:
>>> print("\n")
>>> print("query: %s" % record.query[:100])
>>> for align in record.alignments:
>>> for hsp in align.hsps:
>>> if hsp.expect < E_VALUE_THRESH:
>>> print("match: %s " % align.title[:100])这将产生以下输出：
query: gnl|alu|Z15030_HSAL001056 (Alu-J)
match: gnl|alu|Z15030_HSAL001056 (Alu-J)
match: gnl|alu|L12964_HSAL003860 (Alu-J)
match: gnl|alu|L13042_HSAL003863 (Alu-FLA?)
match: gnl|alu|M86249_HSAL001462 (Alu-FLA?)
match: gnl|alu|M29484_HSAL002265 (Alu-J)
query: gnl|alu|D00596_HSAL003180 (Alu-Sx)
match: gnl|alu|D00596_HSAL003180 (Alu-Sx)
match: gnl|alu|J03071_HSAL001860 (Alu-J)
match: gnl|alu|X72409_HSAL005025 (Alu-Sx)
query: gnl|alu|X55502_HSAL000745 (Alu-J)
match: gnl|alu|X55502_HSAL000745 (Alu-J)Biopython - Entrez Database
Entrez 是 NCBI 提供的在线搜索系统。它提供对几乎所有已知的分子生物学数据库的访问，并支持布尔运算符和字段搜索的集成全局查询。它返回来自所有数据库的结果，其中包含信息，例如来自每个数据库的点击数、带有指向原始数据库的链接的记录等。
下面列出了可以通过 Entrez 访问的一些流行数据库 -
-
Pubmed
-
Pubmed Central
-
Nucleotide (GenBank Sequence Database)
-
Protein (Sequence Database)
-
Genome (Whole Genome Database)
-
结构（三维大分子结构）
-
Taxonomy (Organisms in GenBank)
-
SNP (Single Nucleotide Polymorphism)
-
UniGene（以基因为导向的转录序列簇）
-
CDD（保守蛋白结构域数据库）
-
3D 域（来自 Entrez 结构的域）
除了以上数据库之外，Entrez 还提供了更多数据库来执行字段搜索。
Biopython 提供了一个 Entrez 特定的模块 Bio.Entrez 来访问 Entrez 数据库。让我们在本章学习如何使用 Biopython 访问 Entrez -
Database Connection Steps
要添加 Entrez 的功能，请导入以下模块 -
>>> from Bio import Entrez接下来，设置您的电子邮件以识别与下面给出的代码相连的是谁 -
>>> Entrez.email = '<youremail>'然后，设置 Entrez 工具参数，默认情况下，它为 Biopython。
>>> Entrez.tool = 'Demoscript'现在， call einfo function to find index term counts, last update, and available links for each database 如下所示 -
>>> info = Entrez.einfo()einfo 方法返回一个对象，它可以通过其 read 方法访问信息，如下所示 -
>>> data = info.read()
>>> print(data)
<?xml version = "1.0" encoding = "UTF-8" ?>
<!DOCTYPE eInfoResult PUBLIC "-//NLM//DTD einfo 20130322//EN"
"https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20130322/einfo.dtd">
<eInfoResult>
<DbList>
<DbName>pubmed</DbName>
<DbName>protein</DbName>
<DbName>nuccore</DbName>
<DbName>ipg</DbName>
<DbName>nucleotide</DbName>
<DbName>nucgss</DbName>
<DbName>nucest</DbName>
<DbName>structure</DbName>
<DbName>sparcle</DbName>
<DbName>genome</DbName>
<DbName>annotinfo</DbName>
<DbName>assembly</DbName>
<DbName>bioproject</DbName>
<DbName>biosample</DbName>
<DbName>blastdbinfo</DbName>
<DbName>books</DbName>
<DbName>cdd</DbName>
<DbName>clinvar</DbName>
<DbName>clone</DbName>
<DbName>gap</DbName>
<DbName>gapplus</DbName>
<DbName>grasp</DbName>
<DbName>dbvar</DbName>
<DbName>gene</DbName>
<DbName>gds</DbName>
<DbName>geoprofiles</DbName>
<DbName>homologene</DbName>
<DbName>medgen</DbName>
<DbName>mesh</DbName>
<DbName>ncbisearch</DbName>
<DbName>nlmcatalog</DbName>
<DbName>omim</DbName>
<DbName>orgtrack</DbName>
<DbName>pmc</DbName>
<DbName>popset</DbName>
<DbName>probe</DbName>
<DbName>proteinclusters</DbName>
<DbName>pcassay</DbName>
<DbName>biosystems</DbName>
<DbName>pccompound</DbName>
<DbName>pcsubstance</DbName>
<DbName>pubmedhealth</DbName>
<DbName>seqannot</DbName>
<DbName>snp</DbName>
<DbName>sra</DbName>
<DbName>taxonomy</DbName>
<DbName>biocollections</DbName>
<DbName>unigene</DbName>
<DbName>gencoll</DbName>
<DbName>gtr</DbName>
</DbList>
</eInfoResult>数据以 XML 格式存在，要将数据作为 python 对象获得，请使用 Entrez.read 方法，只要调用 Entrez.einfo() 方法 -
>>> info = Entrez.einfo()
>>> record = Entrez.read(info)此处，record 是一个字典，它有一个键 DbList，如下所示 -
>>> record.keys()
[u'DbList']访问 DbList 键返回下面显示的数据库名称列表 -
>>> record[u'DbList']
['pubmed', 'protein', 'nuccore', 'ipg', 'nucleotide', 'nucgss',
'nucest', 'structure', 'sparcle', 'genome', 'annotinfo', 'assembly',
'bioproject', 'biosample', 'blastdbinfo', 'books', 'cdd', 'clinvar',
'clone', 'gap', 'gapplus', 'grasp', 'dbvar', 'gene', 'gds', 'geoprofiles',
'homologene', 'medgen', 'mesh', 'ncbisearch', 'nlmcatalog', 'omim',
'orgtrack', 'pmc', 'popset', 'probe', 'proteinclusters', 'pcassay',
'biosystems', 'pccompound', 'pcsubstance', 'pubmedhealth', 'seqannot',
'snp', 'sra', 'taxonomy', 'biocollections', 'unigene', 'gencoll', 'gtr']
>>>基本上，Entrez 模块解析 Entrez 搜索系统返回的 XML 并将其作为 python 字典和列表提供。
Search Database
要搜索 Entrez 数据库的任何一个，我们可使用 Bio.Entrez.esearch() 模块。它在下面进行了定义 −
>>> info = Entrez.einfo()
>>> info = Entrez.esearch(db = "pubmed",term = "genome")
>>> record = Entrez.read(info)
>>>print(record)
DictElement({u'Count': '1146113', u'RetMax': '20', u'IdList':
['30347444', '30347404', '30347317', '30347292',
'30347286', '30347249', '30347194', '30347187',
'30347172', '30347088', '30347075', '30346992',
'30346990', '30346982', '30346980', '30346969',
'30346962', '30346954', '30346941', '30346939'],
u'TranslationStack': [DictElement({u'Count':
'927819', u'Field': 'MeSH Terms', u'Term': '"genome"[MeSH Terms]',
u'Explode': 'Y'}, attributes = {})
, DictElement({u'Count': '422712', u'Field':
'All Fields', u'Term': '"genome"[All Fields]', u'Explode': 'N'}, attributes = {}),
'OR', 'GROUP'], u'TranslationSet': [DictElement({u'To': '"genome"[MeSH Terms]
OR "genome"[All Fields]', u'From': 'genome'}, attributes = {})], u'RetStart': '0',
u'QueryTranslation': '"genome"[MeSH Terms] OR "genome"[All Fields]'},
attributes = {})
>>>如果您分配了不正确的 db，则它会返回
>>> info = Entrez.esearch(db = "blastdbinfo",term = "books")
>>> record = Entrez.read(info)
>>> print(record)
DictElement({u'Count': '0', u'RetMax': '0', u'IdList': [],
u'WarningList': DictElement({u'OutputMessage': ['No items found.'],
u'PhraseIgnored': [], u'QuotedPhraseNotFound': []}, attributes = {}),
u'ErrorList': DictElement({u'FieldNotFound': [], u'PhraseNotFound':
['books']}, attributes = {}), u'TranslationSet': [], u'RetStart': '0',
u'QueryTranslation': '(books[All Fields])'}, attributes = {})如果您想跨数据库搜索，则可以使用 Entrez.egquery 。这与 Entrez.esearch 类似，只不过它只需要指定关键字并跳过数据库参数即可。
>>>info = Entrez.egquery(term = "entrez")
>>> record = Entrez.read(info)
>>> for row in record["eGQueryResult"]:
... print(row["DbName"], row["Count"])
...
pubmed 458
pmc 12779 mesh 1
...
...
...
biosample 7
biocollections 0Fetch Records
Entrez 提供了一种特殊方法 efetch，用于从 Entrez 中搜索和下载记录的详细信息。请考虑以下简单示例 −
>>> handle = Entrez.efetch(
db = "nucleotide", id = "EU490707", rettype = "fasta")现在，我们可以使用 SeqIO 对象简单地读取记录
>>> record = SeqIO.read( handle, "fasta" )
>>> record
SeqRecord(seq = Seq('ATTTTTTACGAACCTGTGGAAATTTTTGGTTATGACAATAAATCTAGTTTAGTA...GAA',
SingleLetterAlphabet()), id = 'EU490707.1', name = 'EU490707.1',
description = 'EU490707.1
Selenipedium aequinoctiale maturase K (matK) gene, partial cds; chloroplast',
dbxrefs = [])Biopython - PDB Module
Biopython 提供 Bio.PDB 模块来控制多肽结构。PDB（蛋白质数据库）是当前在线可用的最大的蛋白质结构资源。它承载了很多不同的蛋白质结构，包括蛋白质-蛋白质、蛋白质 DNA、蛋白质 RNA 复合物。
为加载 PDB，键入以下命令：
from Bio.PDB import *Protein Structure File Formats
PDB 将蛋白质结构分发为三种不同的格式：
-
XML 为基础的文件格式，不被 Biopython 支持
-
pdb 文件格式，为经过特殊格式化的文本文件
-
PDBx/mmCIF files format
蛋白质数据库分发的 PDB 文件可能包含导致它们变得不明确或难以解析的格式错误。Bio.PDB 模块试图自动处理这些错误。
Bio.PDB 模块实现了两个不同的解析器，一个是 mmCIF 格式，另一个是 pdb 格式。
我们来详细了解如何解析每个格式：
mmCIF Parser
让我们使用以下命令从 PDB 服务器下载 mmCIF 格式的示例数据库：
>>> pdbl = PDBList()
>>> pdbl.retrieve_pdb_file('2FAT', pdir = '.', file_format = 'mmCif')这会从服务器下载指定文件（2fat.cif）并将其存储于当前的工作目录中。
PDBList 在此处提供从在线 PDB FTP 服务器列出和下载文件的选项。retrieve_pdb_file 方法需要不带扩展名的待下载文件的名称。retrieve_pdb_file 还具有指定下载目录、pdir 和文件格式、file_format 的选项。文件格式的可能值如下：
-
“mmCif” (default, PDBx/mmCif file)
-
“pdb” (format PDB)
-
“xml” (PMDML/XML format)
-
“mmtf” (highly compressed)
-
“bundle”（大结构的 PDB 格式化存档）
若要加载 cif 文件，请按如下所示使用 Bio.MMCIF.MMCIFParser −
>>> parser = MMCIFParser(QUIET = True)
>>> data = parser.get_structure("2FAT", "2FAT.cif")此处，QUIET 在解析文件时禁止生成警告 get_structure will parse the file and return the structure with id as 2FAT （第一个参数）。
运行以上命令后，它将解析文件并打印可能存在的警告（如果存在）。
现在，使用下列命令检查结构 −
>>> data
<Structure id = 2FAT>
To get the type, use type method as specified below,
>>> print(type(data))
<class 'Bio.PDB.Structure.Structure'>我们已成功解析文件并获得了蛋白质结构。我们将在后面的章节了解蛋白质结构的详细信息以及如何获取蛋白质结构。
PDB Parser
让我们使用以下命令从 pdb 服务器中下载 PDB 格式的示例数据库 −
>>> pdbl = PDBList()
>>> pdbl.retrieve_pdb_file('2FAT', pdir = '.', file_format = 'pdb')这将从服务器下载指定的文件（pdb2fat.ent），并将其存储在当前工作目录中。
若要加载 pdb 文件，请按如下所示使用 Bio.PDB.PDBParser −
>>> parser = PDBParser(PERMISSIVE = True, QUIET = True)
>>> data = parser.get_structure("2fat","pdb2fat.ent")在此，get_structure 与 MMCIFParser 类似。PERMISSIVE 选项尝试尽可能灵活地解析蛋白质数据。
现在，使用给出的代码片段检查结构及其类型 −
>>> data
<Structure id = 2fat>
>>> print(type(data))
<class 'Bio.PDB.Structure.Structure'>那么，头部结构存储字典信息。为执行此操作，请键入以下命令 −
>>> print(data.header.keys()) dict_keys([
'name', 'head', 'deposition_date', 'release_date', 'structure_method', 'resolution',
'structure_reference', 'journal_reference', 'author', 'compound', 'source',
'keywords', 'journal'])
>>>为获取名称，请使用以下代码 −
>>> print(data.header["name"])
an anti-urokinase plasminogen activator receptor (upar) antibody: crystal
structure and binding epitope
>>>您还可以使用以下代码检查日期和分辨率 −
>>> print(data.header["release_date"]) 2006-11-14
>>> print(data.header["resolution"]) 1.77PDB Structure
PDB 结构由一个模型组成，其中包含两个链。
-
L 链，包含许多残基
-
H 链，包含许多残基
每个残基由多个原子组成，每个原子具有由 (x，y，z) 坐标表示的三维位置。
让我们在下面的小节中了解如何详细获取原子结构 −
Model
Structure.get_models() 方法返回模型上的迭代器。其定义如下 −
>>> model = data.get_models()
>>> model
<generator object get_models at 0x103fa1c80>
>>> models = list(model)
>>> models [<Model id = 0>]
>>> type(models[0])
<class 'Bio.PDB.Model.Model'>此处，一个模型描述了恰好一个 3D 构象。它包含一个或多个链。
Chain
Model.get_chain() 方法返回一个链的迭代器。它在下面定义 −
>>> chains = list(models[0].get_chains())
>>> chains
[<Chain id = L>, <Chain id = H>]
>>> type(chains[0])
<class 'Bio.PDB.Chain.Chain'>此处，Chain 描述了一个适当的多肽结构，即一个连续的结合残留序列。
Biopython - Motif Objects
序列基序是一个核苷酸或氨基酸序列模式。序列基序是由可能不相邻的氨基酸的三维排列形成的。Biopython 提供了一个单独的模块 Bio.motifs，用于访问如下指定的序列基序的功能 −
from Bio import motifsCreating Simple DNA Motif
让我们使用以下命令创建一个简单的 DNA 基序序列 −
>>> from Bio import motifs
>>> from Bio.Seq import Seq
>>> DNA_motif = [ Seq("AGCT"),
... Seq("TCGA"),
... Seq("AACT"),
... ]
>>> seq = motifs.create(DNA_motif)
>>> print(seq) AGCT TCGA AACT要计数序列值，请使用以下命令 −
>>> print(seq.counts)
0 1 2 3
A: 2.00 1.00 0.00 1.00
C: 0.00 1.00 2.00 0.00
G: 0.00 1.00 1.00 0.00
T: 1.00 0.00 0.00 2.00使用以下代码来计数序列中的“A” −
>>> seq.counts["A", :]
(2, 1, 0, 1)如果您想访问计数列，请使用以下命令 −
>>> seq.counts[:, 3]
{'A': 1, 'C': 0, 'T': 2, 'G': 0}Creating a Sequence Logo
我们现在将讨论如何创建序列标识。
考虑以下序列 −
AGCTTACG
ATCGTACC
TTCCGAAT
GGTACGTA
AAGCTTGG您可以使用以下链接创建自己的标识 − http://weblogo.berkeley.edu/
添加上面序列并创建一个新标识，并将图像命名为 seq.png 保存到您的 Biopython 文件夹。
seq.png
在创建图像后，现在运行以下命令 −
>>> seq.weblogo("seq.png")此 DNA 序列基序被表示为 LexA 结合基序的序列标识。
JASPAR Database
JASPAR 是最流行的数据库之一。它为阅读、写入和扫描序列提供了任何基序格式的设施。它为每个基序存储元信息。 The module Bio.motifs contains a specialized class jaspar.Motif to represent meta-information attributes 。
它拥有以下显著的属性类型 −
-
matrix_id − 唯一的 JASPAR 基序 ID
-
名称 − 基序名称
-
tf_family − 基序家族，例如“螺旋-环状-螺旋”
-
data_type − 基序中使用的数据类型
让我们在 biopython 文件夹中创建一个名为 sample.sites 的 JASPAR 位点格式。它在下面定义 −
sample.sites
>MA0001 ARNT 1
AACGTGatgtccta
>MA0001 ARNT 2
CAGGTGggatgtac
>MA0001 ARNT 3
TACGTAgctcatgc
>MA0001 ARNT 4
AACGTGacagcgct
>MA0001 ARNT 5
CACGTGcacgtcgt
>MA0001 ARNT 6
cggcctCGCGTGc在上面的文件中，我们创建了基序实例。现在，让我们从上述实例创建一个基序对象 −
>>> from Bio import motifs
>>> with open("sample.sites") as handle:
... data = motifs.read(handle,"sites")
...
>>> print(data)
TF name None
Matrix ID None
Matrix:
0 1 2 3 4 5
A: 2.00 5.00 0.00 0.00 0.00 1.00
C: 3.00 0.00 5.00 0.00 0.00 0.00
G: 0.00 1.00 1.00 6.00 0.00 5.00
T: 1.00 0.00 0.00 0.00 6.00 0.00此处，data 从 sample.sites 文件读取所有基序实例。
要打印 data 中的所有实例，请使用以下命令 −
>>> for instance in data.instances:
... print(instance)
...
AACGTG
CAGGTG
TACGTA
AACGTG
CACGTG
CGCGTG使用以下命令统计所有值 −
>>> print(data.counts)
0 1 2 3 4 5
A: 2.00 5.00 0.00 0.00 0.00 1.00
C: 3.00 0.00 5.00 0.00 0.00 0.00
G: 0.00 1.00 1.00 6.00 0.00 5.00
T: 1.00 0.00 0.00 0.00 6.00 0.00
>>>Biopython - BioSQL Module
BioSQL 是一种通用数据库架构，主要设计为存储所有 RDBMS 引擎的序列及其相关数据。它的设计方式是它保存来自所有流行的生物信息学数据库（如 GenBank、Swissport 等）的数据。它还可以用于存储内部数据。
BioSQL 目前为以下数据库提供特定架构 −
-
MySQL (biosqldb-mysql.sql)
-
PostgreSQL (biosqldb-pg.sql)
-
Oracle (biosqldb-ora/*.sql)
-
SQLite (biosqldb-sqlite.sql)
它还为基于 Java 的 HSQLDB 和 Derby 数据库提供最小的支持。
BioPython 提供非常简单、容易和高级的 ORM 功能，用于处理基于 BioSQL 的数据库。 BioPython provides a module, BioSQL 以执行以下功能 −
-
Create/remove a BioSQL database
-
连接到 BioSQL 数据库
-
解析序列数据库，如 GenBank、Swisport、BLAST 结果、Entrez 结果等，并直接将其加载到 BioSQL 数据库
-
从 BioSQL 数据库中获取序列数据
-
从 NCBI BLAST 中获取分类学数据并将其存储在 BioSQL 数据库中
-
对 BioSQL 数据库运行任何 SQL 查询
Overview of BioSQL Database Schema
在深入了解 BioSQL 之前，让我们了解 BioSQL 架构的基础知识。BioSQL 架构提供了 25+ 表来保存序列数据、序列特征、序列类别/本体和分类信息。一些重要的表如下所示 −
-
biodatabase
-
bioentry
-
biosequence
-
seqfeature
-
taxon
-
taxon_name
-
antology
-
term
-
dxref
Creating a BioSQL Database
在本节中，让我们使用由 BioSQL 团队提供的架构创建一个示例 BioSQL 数据库 biosql。我们将使用 SQLite 数据库，因为它很容易上手，并且没有复杂的设置。
这里，我们将按照以下步骤创建基于 SQLite 的 BioSQL 数据库。
Step 1 - 下载 SQLite 数据库引擎并安装它。
Step 2 - 从 GitHub URL 下载 BioSQL 项目。 https://github.com/biosql/biosql
Step 3 - 打开控制台，使用 mkdir 创建一个目录并进入其中。
cd /path/to/your/biopython/sample
mkdir sqlite-biosql
cd sqlite-biosqlStep 4 - 运行以下命令以创建一个新的 SQLite 数据库。
> sqlite3.exe mybiosql.db
SQLite version 3.25.2 2018-09-25 19:08:10
Enter ".help" for usage hints.
sqlite>Step 5 - 从 BioSQL 项目（/sql/biosqldb-sqlite.sql`）复制 biosqldb-sqlite.sql 文件，并将其存储在当前目录中。
Step 6 - 运行以下命令以创建所有表。
sqlite> .read biosqldb-sqlite.sql现在，所有表均已在我们的新数据库中创建。
Step 7 - 运行以下命令以查看数据库中的所有新表。
sqlite> .headers on
sqlite> .mode column
sqlite> .separator ROW "\n"
sqlite> SELECT name FROM sqlite_master WHERE type = 'table';
biodatabase
taxon
taxon_name
ontology
term
term_synonym
term_dbxref
term_relationship
term_relationship_term
term_path
bioentry
bioentry_relationship
bioentry_path
biosequence
dbxref
dbxref_qualifier_value
bioentry_dbxref
reference
bioentry_reference
comment
bioentry_qualifier_value
seqfeature
seqfeature_relationship
seqfeature_path
seqfeature_qualifier_value
seqfeature_dbxref
location
location_qualifier_value
sqlite>前三个命令是配置命令，用于配置 SQLite 以格式化方式显示结果。
Step 8 - 复制 BioPython 团队提供的样例 GenBank 文件，ls_orchid.gbk https://raw.githubusercontent.com/biopython/biopython/master/Doc/examples/ls_orchid.gbk 到当前目录，并将其另存为 orchid.gbk。
Step 9 - 使用以下代码创建 Python 脚本，load_orchid.py 并执行它。
from Bio import SeqIO
from BioSQL import BioSeqDatabase
import os
server = BioSeqDatabase.open_database(driver = 'sqlite3', db = "orchid.db")
db = server.new_database("orchid")
count = db.load(SeqIO.parse("orchid.gbk", "gb"), True) server.commit()
server.close()以上代码解析了文件中的记录，并将其转换为 Python 对象并插入到 BioSQL 数据库中。我们将在后面的章节中分析代码。
最后，我们创建了一个新的 BioSQL 数据库，并将一些示例数据加载到其中。我们将在下一章讨论重要表格。
Simple ER Diagram
biodatabase 表位于层次结构的顶部，其主要目的是将一组序列数据组织成一个组/虚拟数据库。 Every entry in the biodatabase refers to a separate database and it does not mingle with another database. BioSQL 数据库中的所有相关表都引用 biodatabase 条目。
bioentry 表保存有关序列的所有详细信息，但序列数据除外。特定 bioentry 的序列数据将存储在 biosequence 表中。
taxon 和 taxon_name 是分类信息，并且每个条目都引用该表以指定其分类信息。
在理解了模式之后，让我们在下一节中研究一些查询。
BioSQL Queries
让我们深入了解一些 SQL 查询，以更好地了解数据的组织方式以及表彼此之间的关系。在继续之前，让我们使用以下命令打开数据库并设置一些格式化命令 -
> sqlite3 orchid.db
SQLite version 3.25.2 2018-09-25 19:08:10
Enter ".help" for usage hints.
sqlite> .header on
sqlite> .mode columns.header and .mode are formatting options to better visualize the data 。您还可以使用任何 SQLite 编辑器来运行查询。
如下所示，列出系统中可用的虚拟序列数据库 −
select
*
from
biodatabase;
*** Result ***
sqlite> .width 15 15 15 15
sqlite> select * from biodatabase;
biodatabase_id name authority description
--------------- --------------- --------------- ---------------
1 orchid
sqlite>这里，我们只有一个数据库， orchid 。
使用以下给出的代码列出数据库 orchid 中可用的条目（前 3 项）
select
be.*,
bd.name
from
bioentry be
inner join
biodatabase bd
on bd.biodatabase_id = be.biodatabase_id
where
bd.name = 'orchid' Limit 1,
3;
*** Result ***
sqlite> .width 15 15 10 10 10 10 10 50 10 10
sqlite> select be.*, bd.name from bioentry be inner join biodatabase bd on
bd.biodatabase_id = be.biodatabase_id where bd.name = 'orchid' Limit 1,3;
bioentry_id biodatabase_id taxon_id name accession identifier division description version name
--------------- --------------- ---------- ---------- ---------- ---------- ----------
---------- ---------- ----------- ---------- --------- ---------- ----------
2 1 19 Z78532 Z78532 2765657 PLN
C.californicum 5.8S rRNA gene and ITS1 and ITS2 DN 1
orchid
3 1 20 Z78531 Z78531 2765656 PLN
C.fasciculatum 5.8S rRNA gene and ITS1 and ITS2 DN 1
orchid
4 1 21 Z78530 Z78530 2765655 PLN
C.margaritaceum 5.8S rRNA gene and ITS1 and ITS2 D 1
orchid
sqlite>使用给定的代码列出与条目（登录号 − Z78530、名称 − C. fasciculatum 5.8S rRNA 基因和 ITS1 及 ITS2 DNA）相关的序列详细信息 −
select
substr(cast(bs.seq as varchar), 0, 10) || '...' as seq,
bs.length,
be.accession,
be.description,
bd.name
from
biosequence bs
inner join
bioentry be
on be.bioentry_id = bs.bioentry_id
inner join
biodatabase bd
on bd.biodatabase_id = be.biodatabase_id
where
bd.name = 'orchid'
and be.accession = 'Z78532';
*** Result ***
sqlite> .width 15 5 10 50 10
sqlite> select substr(cast(bs.seq as varchar), 0, 10) || '...' as seq,
bs.length, be.accession, be.description, bd.name from biosequence bs inner
join bioentry be on be.bioentry_id = bs.bioentry_id inner join biodatabase bd
on bd.biodatabase_id = be.biodatabase_id where bd.name = 'orchid' and
be.accession = 'Z78532';
seq length accession description name
------------ ---------- ---------- ------------ ------------ ---------- ---------- -----------------
CGTAACAAG... 753 Z78532 C.californicum 5.8S rRNA gene and ITS1 and ITS2 DNA orchid
sqlite>使用以下代码获取与条目（登录号 − Z78530、名称 − C. fasciculatum 5.8S rRNA 基因和 ITS1 及 ITS2 DNA）相关的完整序列 −
select
bs.seq
from
biosequence bs
inner join
bioentry be
on be.bioentry_id = bs.bioentry_id
inner join
biodatabase bd
on bd.biodatabase_id = be.biodatabase_id
where
bd.name = 'orchid'
and be.accession = 'Z78532';
*** Result ***
sqlite> .width 1000
sqlite> select bs.seq from biosequence bs inner join bioentry be on
be.bioentry_id = bs.bioentry_id inner join biodatabase bd on bd.biodatabase_id =
be.biodatabase_id where bd.name = 'orchid' and be.accession = 'Z78532';
seq
----------------------------------------------------------------------------------------
----------------------------
CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAACAGAATATATGATCGAGTGAATCT
GGAGGACCTGTGGTAACTCAGCTCGTCGTGGCACTGCTTTTGTCGTGACCCTGCTTTGTTGTTGGGCCTCC
TCAAGAGCTTTCATGGCAGGTTTGAACTTTAGTACGGTGCAGTTTGCGCCAAGTCATATAAAGCATCACTGATGAATGACATTATTGT
CAGAAAAAATCAGAGGGGCAGTATGCTACTGAGCATGCCAGTGAATTTTTATGACTCTCGCAACGGATATCTTGGCTC
TAACATCGATGAAGAACGCAG
sqlite>列出与生物数据库相关的分类单元，兰花
select distinct
tn.name
from
biodatabase d
inner join
bioentry e
on e.biodatabase_id = d.biodatabase_id
inner join
taxon t
on t.taxon_id = e.taxon_id
inner join
taxon_name tn
on tn.taxon_id = t.taxon_id
where
d.name = 'orchid' limit 10;
*** Result ***
sqlite> select distinct tn.name from biodatabase d inner join bioentry e on
e.biodatabase_id = d.biodatabase_id inner join taxon t on t.taxon_id =
e.taxon_id inner join taxon_name tn on tn.taxon_id = t.taxon_id where d.name =
'orchid' limit 10;
name
------------------------------
Cypripedium irapeanum
Cypripedium californicum
Cypripedium fasciculatum
Cypripedium margaritaceum
Cypripedium lichiangense
Cypripedium yatabeanum
Cypripedium guttatum
Cypripedium acaule
pink lady's slipper
Cypripedium formosanum
sqlite>Load Data into BioSQL Database
让我们学习如何在本章中将序列数据加载到 BioSQL 数据库中。我们在前一节中已经有了将数据加载到数据库中的代码，代码如下 −
from Bio import SeqIO
from BioSQL import BioSeqDatabase
import os
server = BioSeqDatabase.open_database(driver = 'sqlite3', db = "orchid.db")
DBSCHEMA = "biosqldb-sqlite.sql"
SQL_FILE = os.path.join(os.getcwd(), DBSCHEMA)
server.load_database_sql(SQL_FILE)
server.commit()
db = server.new_database("orchid")
count = db.load(SeqIO.parse("orchid.gbk", "gb"), True) server.commit()
server.close()我们将深入了解代码的每一行及其目的 −
Line 1 − 加载 SeqIO 模块。
Line 2 − 加载 BioSeqDatabase 模块。此模块提供与 BioSQL 数据库交互的所有功能。
Line 3 − 加载 os 模块。
Line 5 − open_database 使用配置的驱动程序 (driver) 打开指定的数据库 (db) 并返回对 BioSQL 数据库 (server) 的句柄。Biopython 支持 sqlite、mysql、postgresql 和 oracle 数据库。
Line 6-10 − load_database_sql 方法从外部文件加载 sql 并执行它。commit 方法提交事务。我们可以跳过此步骤，因为我们已经使用 schema 创建了数据库。
Line 12 − new_database 方法创建新的虚拟数据库 orchid 并返回句柄 db 以对 orchid 数据库执行命令。
Line 13 − load 方法将序列条目（可迭代 SeqRecord）加载到 orchid 数据库中。SqlIO.parse 解析 GenBank 数据库并以可迭代 SeqRecord 的形式返回其中的所有序列。load 方法的第二个参数 (True) 指示它从 NCBI blast 网站获取序列数据的分类单元详细信息，如果系统中尚未提供这些详细信息。
Line 14 − commit 提交事务。
Line 15 − close 关闭数据库连接并销毁服务器句柄。
Fetch the Sequence Data
让我们按如下所示从 orchid 数据库中获取标识符为 2765658 的序列 −
from BioSQL import BioSeqDatabase
server = BioSeqDatabase.open_database(driver = 'sqlite3', db = "orchid.db")
db = server["orchid"]
seq_record = db.lookup(gi = 2765658)
print(seq_record.id, seq_record.description[:50] + "...")
print("Sequence length %i," % len(seq_record.seq))此处，server["orchid"] 返回从虚拟数据库 orchid 中提取数据的句柄。 lookup 方法提供了一个根据标准选择序列的选项，我们已经使用标识符 2765658 选择了序列。 lookup 以 SeqRecordObject 的形式返回序列信息。由于我们已经知道如何使用 SeqRecord，因此很容易从中获取数据。
Biopython - Population Genetics
种群遗传学在进化论中扮演着重要的角色。它分析物种之间的遗传差异以及同一物种内的两个或更多个体之间的遗传差异。
Biopython 为种群遗传学提供 Bio.PopGen 模块，并主要支持 `GenePop，一个由 Michel Raymond 和 Francois Rousset 开发的流行遗传学包。
A simple parser
让我们编写一个简单的应用程序来解析 GenePop 格式并理解该概念。
在下面给出的链接中下载 Biopython 团队提供的 genePop 文件 − https://raw.githubusercontent.com/biopython/biopython/master/Tests/PopGen/c3line.gen
使用下面的代码段加载 GenePop 模块 −
from Bio.PopGen import GenePop按照下面使用 GenePop.read 方法解析文件 −
record = GenePop.read(open("c3line.gen"))显示下面给出的基因座和种群信息 −
>>> record.loci_list
['136255903', '136257048', '136257636']
>>> record.pop_list
['4', 'b3', '5']
>>> record.populations
[[('1', [(3, 3), (4, 4), (2, 2)]), ('2', [(3, 3), (3, 4), (2, 2)]),
('3', [(3, 3), (4, 4), (2, 2)]), ('4', [(3, 3), (4, 3), (None, None)])],
[('b1', [(None, None), (4, 4), (2, 2)]), ('b2', [(None, None), (4, 4), (2, 2)]),
('b3', [(None, None), (4, 4), (2, 2)])],
[('1', [(3, 3), (4, 4), (2, 2)]), ('2', [(3, 3), (1, 4), (2, 2)]),
('3', [(3, 2), (1, 1), (2, 2)]), ('4',
[(None, None), (4, 4), (2, 2)]), ('5', [(3, 3), (4, 4), (2, 2)])]]
>>>此处，文件中有三个基因座和三组种群：第一组种群有 4 条记录，第二组种群有 3 条记录，第三组种群有 5 条记录。record.populations 显示所有种群集以及每个基因座的等位基因数据。
Manipulate the GenePop file
Biopython 提供移除基因座和种群数据的选项。
Remove a population set by position,
>>> record.remove_population(0)
>>> record.populations
[[('b1', [(None, None), (4, 4), (2, 2)]),
('b2', [(None, None), (4, 4), (2, 2)]),
('b3', [(None, None), (4, 4), (2, 2)])],
[('1', [(3, 3), (4, 4), (2, 2)]),
('2', [(3, 3), (1, 4), (2, 2)]),
('3', [(3, 2), (1, 1), (2, 2)]),
('4', [(None, None), (4, 4), (2, 2)]),
('5', [(3, 3), (4, 4), (2, 2)])]]
>>>Remove a locus by position,
>>> record.remove_locus_by_position(0)
>>> record.loci_list
['136257048', '136257636']
>>> record.populations
[[('b1', [(4, 4), (2, 2)]), ('b2', [(4, 4), (2, 2)]), ('b3', [(4, 4), (2, 2)])],
[('1', [(4, 4), (2, 2)]), ('2', [(1, 4), (2, 2)]),
('3', [(1, 1), (2, 2)]), ('4', [(4, 4), (2, 2)]), ('5', [(4, 4), (2, 2)])]]
>>>Remove a locus by name,
>>> record.remove_locus_by_name('136257636') >>> record.loci_list
['136257048']
>>> record.populations
[[('b1', [(4, 4)]), ('b2', [(4, 4)]), ('b3', [(4, 4)])],
[('1', [(4, 4)]), ('2', [(1, 4)]),
('3', [(1, 1)]), ('4', [(4, 4)]), ('5', [(4, 4)])]]
>>>Interface with GenePop Software
Biopython 提供了与 GenePop 软件交互的接口，由此公开了该软件的许多功能。为此，使用了 Bio.PopGen.GenePop 模块。EasyController 是一个易于使用的接口。让我们了解如何解析 GenePop 文件以及使用 EasyController 执行一些分析。
首先，安装 GenePop 软件并将安装文件夹放入系统路径中。要获取关于 GenePop 文件的基本信息，请创建一个 EasyController 对象，然后按照下面指定的调用 get_basic_info 方法−
>>> from Bio.PopGen.GenePop.EasyController import EasyController
>>> ec = EasyController('c3line.gen')
>>> print(ec.get_basic_info())
(['4', 'b3', '5'], ['136255903', '136257048', '136257636'])
>>>此处，第一项是种群列表，第二项是基因座列表。
要获取特定基因座的所有等位基因列表，请按如下指定传递基因座名称来调用 get_alleles_all_pops 方法−
>>> allele_list = ec.get_alleles_all_pops("136255903")
>>> print(allele_list)
[2, 3]要按特定种群和基因座获取等位基因列表，请按如下指定传递基因座名称和种群位置来调用 get_alleles −
>>> allele_list = ec.get_alleles(0, "136255903")
>>> print(allele_list)
[]
>>> allele_list = ec.get_alleles(1, "136255903")
>>> print(allele_list)
[]
>>> allele_list = ec.get_alleles(2, "136255903")
>>> print(allele_list)
[2, 3]
>>>同样，EasyController 公开了很多功能：等位基因频率、基因型频率、多基因座 F 统计量、哈代-温伯格平衡、连锁不平衡等。
Biopython - Genome Analysis
基因组是 DNA 的完整集合，包括其所有基因。基因组分析是指研究个体基因及其在遗传中的作用。
Genome Diagram
基因组图以图表的形式展示遗传信息。Biopython 使用 Bio.Graphics.GenomeDiagram 模块来展示 GenomeDiagram。GenomeDiagram 模块需要安装 ReportLab。
Steps for creating a diagram
创建图表的流程通常遵循以下简单模式：
-
为想要展示的每个单独的特性集创建一个 FeatureSet，并向它们添加 Bio.SeqFeature 对象。
-
为想要展示的每个图表创建一个 GraphSet，并向它们添加图表数据。
-
为图表中想要的每个轨道创建一个 Track，并向需要的轨道添加 GraphSets 和 FeatureSets。
-
创建一个 Diagram，并向其中添加 Track。
-
告诉 Diagram 绘制图像。
-
将图像写入文件。
让我们以输入 GenBank 文件为例：
https://raw.githubusercontent.com/biopython/biopython/master/Doc/examples/ls_orchid.gbk 并从 SeqRecord 对象中读取记录，然后最终绘制基因组图表。它在下面进行了说明：
我们首先需要导入所有模块，如下所示：
>>> from reportlab.lib import colors
>>> from reportlab.lib.units import cm
>>> from Bio.Graphics import GenomeDiagram现在，导入 SeqIO 模块来读取数据：
>>> from Bio import SeqIO
record = SeqIO.read("example.gb", "genbank")在这里，记录从 genbank 文件中读取序列。
现在，创建一个空的图表来添加轨道和特性集：
>>> diagram = GenomeDiagram.Diagram(
"Yersinia pestis biovar Microtus plasmid pPCP1")
>>> track = diagram.new_track(1, name="Annotated Features")
>>> feature = track.new_set()现在，我们可以使用下面定义的绿色到灰色的交替颜色应用颜色主题更改：
>>> for feature in record.features:
>>> if feature.type != "gene":
>>> continue
>>> if len(feature) % 2 == 0:
>>> color = colors.blue
>>> else:
>>> color = colors.red
>>>
>>> feature.add_feature(feature, color=color, label=True)现在您可以在屏幕上看到以下响应：
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d3dc90>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d3dfd0>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x1007627d0>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d57290>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d57050>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d57390>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d57590>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d57410>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d57490>
<Bio.Graphics.GenomeDiagram._Feature.Feature object at 0x105d574d0>让我们为上述输入记录绘制一张图表 -
>>> diagram.draw(
format = "linear", orientation = "landscape", pagesize = 'A4',
... fragments = 4, start = 0, end = len(record))
>>> diagram.write("orchid.pdf", "PDF")
>>> diagram.write("orchid.eps", "EPS")
>>> diagram.write("orchid.svg", "SVG")
>>> diagram.write("orchid.png", "PNG")执行上述命令后，您可以在 Biopython 目录中看到已保存的以下图像。
** Result **
genome.png
您还可以通过进行以下更改来以圆形格式绘制图像 -
>>> diagram.draw(
format = "circular", circular = True, pagesize = (20*cm,20*cm),
... start = 0, end = len(record), circle_core = 0.7)
>>> diagram.write("circular.pdf", "PDF")Biopython - Phenotype Microarray
表型被定义为由生物体对特定化学物质或环境表现出的可观察到的特征或性状。表型微阵列同时测量生物体对大量化学物质和环境的反应，并分析数据以了解基因突变、基因特性等。
Biopython 提供了一个优秀的模块 Bio.Phenotype 来分析表型数据。让我们在本章学习如何在该表型数据中解析、内插、抽取和分析表型微阵列数据。
Parsing
表型微阵列数据可以采用两种格式：CSV 和 JSON。Biopython 支持这两种格式。Biopython 解析器解析表型微阵列数据并返回为 PlateRecord 对象的集合。每个 PlateRecord 对象都包含一个 WellRecord 对象的集合。每个 WellRecord 对象以 8 行和 12 列格式保存数据。8 行表示为 A 到 H，12 列表示为 01 到 12。例如，第 4 行和第 6 列表示为 D06。
让我们通过以下示例了解解析的格式和概念−
Step 1 − 下载由 Biopython 团队提供的 Plates.csv 文件− https://raw.githubusercontent.com/biopython/biopython/master/Doc/examples/Plates.csv
Step 2 − 如下加载 phenotpe 模块− Step 3
>>> from Bio import phenotypeStep 3 − 调用 phenotype.parse 方法，传递数据文件和格式选项（“pm-csv”）。它返回可迭代的 PlateRecord 如下所示：
>>> plates = list(phenotype.parse('Plates.csv', "pm-csv"))
>>> plates
[PlateRecord('WellRecord['A01'], WellRecord['A02'], WellRecord['A03'], ..., WellRecord['H12']'),
PlateRecord('WellRecord['A01'], WellRecord['A02'], WellRecord['A03'], ..., WellRecord['H12']'),
PlateRecord('WellRecord['A01'], WellRecord['A02'], WellRecord['A03'], ..., WellRecord['H12']'),
PlateRecord('WellRecord['A01'], WellRecord['A02'],WellRecord['A03'], ..., WellRecord['H12']')]
>>>Step 4 − 从列表中访问第一个平板，如下所示：
>>> plate = plates[0]
>>> plate
PlateRecord('WellRecord['A01'], WellRecord['A02'], WellRecord['A03'], ...,
WellRecord['H12']')
>>>Step 5 − 如前所述，一个平板包含 8 行，每行有 12 个项目。WellRecord 可通过如下指定的两种方式访问：
>>> well = plate["A04"]
>>> well = plate[0, 4]
>>> well WellRecord('(0.0, 0.0), (0.25, 0.0), (0.5, 0.0), (0.75, 0.0),
(1.0, 0.0), ..., (71.75, 388.0)')
>>>Step 6 − 每个孔都有在不同时间点的一系列测量，并且可以使用循环访问，如下所示：
>>> for v1, v2 in well:
... print(v1, v2)
...
0.0 0.0
0.25 0.0
0.5 0.0
0.75 0.0
1.0 0.0
...
71.25 388.0
71.5 388.0
71.75 388.0
>>>Interpolation
插值可以更深入地了解数据。Biopython 提供了内插 WellRecord 数据的方法以获取中间时间点的信息。语法类似于列表索引，因此易于学习。
要获取 20.1 小时的数据，只需按如下指定传递为索引值即可−
>>> well[20.10]
69.40000000000003
>>>我们可以传递开始时间点和结束时间点，如下所述：
>>> well[20:30]
[67.0, 84.0, 102.0, 119.0, 135.0, 147.0, 158.0, 168.0, 179.0, 186.0]
>>>上述命令使用 1 小时间隔以内插 20 小时至 30 小时的数据。默认情况下，间隔为 1 小时，我们可以将其更改为任何值。例如，让我们按照下方指定的方式赋予 15 分钟（0.25 小时）的间隔：
>>> well[20:21:0.25]
[67.0, 73.0, 75.0, 81.0]
>>>Analyze and Extract
Biopython 提供了一种方法，可以使用 Gompertz、Logistic 和 Richards sigmoid 函数分析 WellRecord 数据。默认情况下，拟合方法使用 Gompertz 函数。我们需要调用 WellRecord 对象的拟合方法来完成任务。编码如下：
>>> well.fit()
Traceback (most recent call last):
...
Bio.MissingPythonDependencyError: Install scipy to extract curve parameters.
>>> well.model
>>> getattr(well, 'min') 0.0
>>> getattr(well, 'max') 388.0
>>> getattr(well, 'average_height')
205.42708333333334
>>>Biopython 依赖 scipy 模块来执行高级分析。它将计算最小值、最大值和平均值详情，而无需使用 scipy 模块。
Biopython - Plotting
本章解释了如何绘制序列。在进入该主题之前，我们先了解一下绘图的基础知识。
Plotting
Matplotlib 是一个 Python 绘图库，它以各种格式生成高质量的图形。我们可以创建各种类型的图形，如折线图、直方图、条形图、饼图、散点图等。
*pyLab 是属于 matplotlib 的一个模块，它将数字模块 numpy 与绘图模块 pyplot 结合在一起。*Biopython 使用 pylab 模块绘制序列。为此，我们需要导入以下代码：
import pylab在导入之前，我们需要使用以下给出的命令使用 pip 命令安装 matplotlib 包：
pip install matplotlibSample Input File
在你的 Biopython 目录中创建一个名为 plot.fasta 的示例文件，并添加以下更改：
>seq0 FQTWEEFSRAAEKLYLADPMKVRVVLKYRHVDGNLCIKVTDDLVCLVYRTDQAQDVKKIEKF
>seq1 KYRTWEEFTRAAEKLYQADPMKVRVVLKYRHCDGNLCIKVTDDVVCLLYRTDQAQDVKKIEKFHSQLMRLME
>seq2 EEYQTWEEFARAAEKLYLTDPMKVRVVLKYRHCDGNLCMKVTDDAVCLQYKTDQAQDVKKVEKLHGK
>seq3 MYQVWEEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVCLQYKTDQAQDV
>seq4 EEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVVSYEMRLFGVQKDNFALEHSLL
>seq5 SWEEFAKAAEVLYLEDPMKCRMCTKYRHVDHKLVVKLTDNHTVLKYVTDMAQDVKKIEKLTTLLMR
>seq6 FTNWEEFAKAAERLHSANPEKCRFVTKYNHTKGELVLKLTDDVVCLQYSTNQLQDVKKLEKLSSTLLRSI
>seq7 SWEEFVERSVQLFRGDPNATRYVMKYRHCEGKLVLKVTDDRECLKFKTDQAQDAKKMEKLNNIFF
>seq8 SWDEFVDRSVQLFRADPESTRYVMKYRHCDGKLVLKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM
>seq9 KNWEDFEIAAENMYMANPQNCRYTMKYVHSKGHILLKMSDNVKCVQYRAENMPDLKK
>seq10 FDSWDEFVSKSVELFRNHPDTTRYVVKYRHCEGKLVLKVTDNHECLKFKTDQAQDAKKMEKLine Plot
现在，让我们为上述 fasta 文件创建一个简单的折线图。
Step 1 − 导入 SeqIO 模块来读取 fasta 文件。
>>> from Bio import SeqIOStep 2 − 解析输入文件。
>>> records = [len(rec) for rec in SeqIO.parse("plot.fasta", "fasta")]
>>> len(records)
11
>>> max(records)
72
>>> min(records)
57Step 3 − 让我们导入 pylab 模块。
>>> import pylabStep 4 − 通过分配 x 和 y 轴标签来配置折线图。
>>> pylab.xlabel("sequence length")
Text(0.5, 0, 'sequence length')
>>> pylab.ylabel("count")
Text(0, 0.5, 'count')
>>>Step 5 − 通过设置网格显示来配置折线图。
>>> pylab.grid()Step 6 − 通过调用绘图方法并提供记录作为输入来绘制简单的折线图。
>>> pylab.plot(records)
[<matplotlib.lines.Line2D object at 0x10b6869d 0>]Step 7 − 最后，使用以下命令保存图表。
>>> pylab.savefig("lines.png")
Histogram Chart
直方图用于连续数据，其中柱表示数据范围。绘制直方图与折线图相同，但使用 pylab.plot 除外。相反，使用记录和一些自定义的柱值（5）调用 pylab 模块的 hist 方法。完整的编码如下：
Step 1 − 导入 SeqIO 模块来读取 fasta 文件。
>>> from Bio import SeqIOStep 2 − 解析输入文件。
>>> records = [len(rec) for rec in SeqIO.parse("plot.fasta", "fasta")]
>>> len(records)
11
>>> max(records)
72
>>> min(records)
57Step 3 − 让我们导入 pylab 模块。
>>> import pylabStep 4 − 通过分配 x 和 y 轴标签来配置折线图。
>>> pylab.xlabel("sequence length")
Text(0.5, 0, 'sequence length')
>>> pylab.ylabel("count")
Text(0, 0.5, 'count')
>>>Step 5 − 通过设置网格显示来配置折线图。
>>> pylab.grid()Step 6 − 通过调用绘图方法并提供记录作为输入来绘制简单的折线图。
>>> pylab.hist(records,bins=5)
(array([2., 3., 1., 3., 2.]), array([57., 60., 63., 66., 69., 72.]), <a list
of 5 Patch objects>)
>>>Step 7 − 最后，使用以下命令保存图表。
>>> pylab.savefig("hist.png")
GC Percentage in Sequence
GC 百分比是最常用的分析数据之一，用于比较不同的序列。我们可以使用一组序列的 GC 百分比绘制一个简单的折线图，并立即对其进行比较。在此，我们可以将数据从序列长度更改为 GC 百分比。完整的编码如下：
Step 1 − 导入 SeqIO 模块来读取 fasta 文件。
>>> from Bio import SeqIOStep 2 − 解析输入文件。
>>> from Bio.SeqUtils import GC
>>> gc = sorted(GC(rec.seq) for rec in SeqIO.parse("plot.fasta", "fasta"))Step 3 − 让我们导入 pylab 模块。
>>> import pylabStep 4 − 通过分配 x 和 y 轴标签来配置折线图。
>>> pylab.xlabel("Genes")
Text(0.5, 0, 'Genes')
>>> pylab.ylabel("GC Percentage")
Text(0, 0.5, 'GC Percentage')
>>>Step 5 − 通过设置网格显示来配置折线图。
>>> pylab.grid()Step 6 − 通过调用绘图方法并提供记录作为输入来绘制简单的折线图。
>>> pylab.plot(gc)
[<matplotlib.lines.Line2D object at 0x10b6869d 0>]Step 7 − 最后，使用以下命令保存图表。
>>> pylab.savefig("gc.png")
Biopython - Cluster Analysis
一般来讲，聚类分析会将一堆物体归为同一组。这个概念主要应用于数据挖掘、统计数据分析、机器学习、模式识别、图像分析、生物信息学等领域。我们可以通过各种算法来理解聚类在不同分析中的广泛应用。
根据生物信息学，聚类分析主要用于基因表达数据分析，以找出具有相似基因表达的一组基因。
在本章中，我们将了解 Biopython 中的重要算法，以理解真实数据集上聚类的基本原理。
Biopython 使用 Bio.Cluster 模块来实现所有算法。它支持以下算法 −
-
Hierarchical Clustering
-
K - Clustering
-
Self-Organizing Maps
-
Principal Component Analysis
让我们简单介绍一下上述算法。
Hierarchical Clustering
层级聚类用于通过距离测量将每个节点与其最近邻节点相连，并创建一个簇。Bio.Cluster 节点具有三个属性：左、右和距离。让我们创建一个简单的簇，如下所示 −
>>> from Bio.Cluster import Node
>>> n = Node(1,10)
>>> n.left = 11
>>> n.right = 0
>>> n.distance = 1
>>> print(n)
(11, 0): 1如果你想构建基于树的聚类，请使用以下命令 −
>>> n1 = [Node(1, 2, 0.2), Node(0, -1, 0.5)] >>> n1_tree = Tree(n1)
>>> print(n1_tree)
(1, 2): 0.2
(0, -1): 0.5
>>> print(n1_tree[0])
(1, 2): 0.2让我们使用 Bio.Cluster 模块执行层级聚类。
考虑距离是在一个数组中定义的。
>>> import numpy as np
>>> distance = array([[1,2,3],[4,5,6],[3,5,7]])现在将距离数组添加到树簇中。
>>> from Bio.Cluster import treecluster
>>> cluster = treecluster(distance)
>>> print(cluster)
(2, 1): 0.666667
(-1, 0): 9.66667上述函数返回一个树簇对象。该对象包含节点，其中项目数按行或列聚类。
K - Clustering
它是一种分区算法，分为 k 均值、中位数和 medoids 聚类。让我们简要了解一下每种聚类。
K-means Clustering
这种方法在数据挖掘中很流行。此算法的目标是找出数据中的组，其中组的数量由变量 K 表示。
该算法迭代运行，根据提供 的特征将每个数据点分配到 K 个组中的一个。数据点根据特征相似性进行聚类。
>>> from Bio.Cluster import kcluster
>>> from numpy import array
>>> data = array([[1, 2], [3, 4], [5, 6]])
>>> clusterid, error,found = kcluster(data)
>>> print(clusterid) [0 0 1]
>>> print(found)
1K-medoids Clustering
此方法基于给定的一组项目，使用用户传递的距离矩阵和簇数。
考虑如下定义的距离矩阵 −
>>> distance = array([[1,2,3],[4,5,6],[3,5,7]])我们可以使用以下命令计算 k 均值聚类 −
>>> from Bio.Cluster import kmedoids
>>> clusterid, error, found = kmedoids(distance)让我们来看一个示例。
kcluster 函数以数据矩阵作为输入，而不是 Seq 实例。您需要将序列转换为矩阵并提供给 kcluster 函数。
将数据转换为仅包含数值元素的矩阵的一种方法是使用 numpy.fromstring 函数。它基本上将序列中的每个字母翻译为其 ASCII 对应字符。
这将创建 kcluster 函数识别并用于对序列进行聚类的已编码序列的二维数组。
>>> from Bio.Cluster import kcluster
>>> import numpy as np
>>> sequence = [ 'AGCT','CGTA','AAGT','TCCG']
>>> matrix = np.asarray([np.fromstring(s, dtype=np.uint8) for s in sequence])
>>> clusterid,error,found = kcluster(matrix)
>>> print(clusterid) [1 0 0 1]Self-Organizing Maps
这种方法是一种人工神经网络。它是由 Kohonen 开发的，通常称为 Kohonen 映射。它根据矩形拓扑将项目组织成群集。
让我们使用如下所示的相同数组距离来创建简单的群集 −
>>> from Bio.Cluster import somcluster
>>> from numpy import array
>>> data = array([[1, 2], [3, 4], [5, 6]])
>>> clusterid,map = somcluster(data)
>>> print(map)
[[[-1.36032469 0.38667395]]
[[-0.41170578 1.35295911]]]
>>> print(clusterid)
[[1 0]
[1 0]
[1 0]]此处， clusterid 是一个有两列的数组，其中，行数等于已聚类的项目数，而 data 是一个具有行或列维度的数组。
Principal Component Analysis
主成分分析对于可视化高维数据非常有用。它是一种使用线性代数和统计学中的简单矩阵操作来计算原始数据到相同数量或更少维度的投影的方法。
主成分分析返回一个包含列均值、坐标、组件和特征值的元组。让我们了解一下这一概念的基础知识。
>>> from numpy import array
>>> from numpy import mean
>>> from numpy import cov
>>> from numpy.linalg import eig
# define a matrix
>>> A = array([[1, 2], [3, 4], [5, 6]])
>>> print(A)
[[1 2]
[3 4]
[5 6]]
# calculate the mean of each column
>>> M = mean(A.T, axis = 1)
>>> print(M)
[ 3. 4.]
# center columns by subtracting column means
>>> C = A - M
>>> print(C)
[[-2. -2.]
[ 0. 0.]
[ 2. 2.]]
# calculate covariance matrix of centered matrix
>>> V = cov(C.T)
>>> print(V)
[[ 4. 4.]
[ 4. 4.]]
# eigendecomposition of covariance matrix
>>> values, vectors = eig(V)
>>> print(vectors)
[[ 0.70710678 -0.70710678]
[ 0.70710678 0.70710678]]
>>> print(values)
[ 8. 0.]让我们将相同的 rectangular 矩阵数据应用于如下定义的 Bio.Cluster 模块 −
>>> from Bio.Cluster import pca
>>> from numpy import array
>>> data = array([[1, 2], [3, 4], [5, 6]])
>>> columnmean, coordinates, components, eigenvalues = pca(data)
>>> print(columnmean)
[ 3. 4.]
>>> print(coordinates)
[[-2.82842712 0. ]
[ 0. 0. ]
[ 2.82842712 0. ]]
>>> print(components)
[[ 0.70710678 0.70710678]
[ 0.70710678 -0.70710678]]
>>> print(eigenvalues)
[ 4. 0.]Biopython - Machine Learning
生物信息学是应用机器学习算法的一个绝佳领域。在这里，我们拥有大量生物体的遗传信息，并且不可能手动分析所有这些信息。如果使用适当的机器学习算法，我们可以从这些数据中提取大量有用的信息。Biopython 提供了一组有用的算法来进行监督式机器学习。
监督式学习基于输入变量 (X) 和输出变量 (Y)。它使用一种算法来了解从输入到输出的映射函数。它在下面定义 −
Y = f(X)这种方法的主要目的是逼近映射函数，当您有新的输入数据 (x) 时，您可以预测该数据的输出变量 (Y)。
Logistic Regression Model
逻辑回归是一种监督式机器学习算法。它用于使用预测变量的加权和来找出 K 个类之间的差异。它计算事件发生的概率，可用于癌症检测。
Biopython 提供 Bio.LogisticRegression 模块，以根据逻辑回归算法预测变量。目前，Biopython 仅针对两个类实现了逻辑回归算法 (K = 2)。
Biopython - Testing Techniques
Biopython 具有大量的测试脚本，可以在不同的条件下测试软件以确保软件没有错误。要运行测试脚本，请下载 Biopython 的源代码，然后运行以下命令 −
python run_tests.py这将运行所有测试脚本，并给出以下输出 −
Python version: 2.7.12 (v2.7.12:d33e0cf91556, Jun 26 2016, 12:10:39)
[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)]
Operating system: posix darwin
test_Ace ... ok
test_Affy ... ok
test_AlignIO ... ok
test_AlignIO_ClustalIO ... ok
test_AlignIO_EmbossIO ... ok
test_AlignIO_FastaIO ... ok
test_AlignIO_MauveIO ... ok
test_AlignIO_PhylipIO ... ok
test_AlignIO_convert ... ok
...........................................
...........................................我们也可以运行专门的测试脚本，如下面指定的 −
python test_AlignIO.py