使用 Bio.motifs 进行序列基序分析
本章概述了 Biopython 中包含的 Bio.motifs 包的功能。它是为参与序列基序分析的人员准备的，因此我假设您熟悉基序分析的基本概念。如果某些内容不清楚，请参阅第 有用链接 节以获取一些相关链接。
本章主要介绍 Biopython 1.61 及更高版本中包含的新 Bio.motifs 包，该包取代了 Biopython 1.50 中引入的旧 Bio.Motif 包，而旧的 Bio.Motif 包又是基于两个更早的 Biopython 模块 Bio.AlignAce 和 Bio.MEME。它提供了它们的大部分功能，并具有统一的基序对象实现。
说到其他库，如果您正在阅读本文，您可能对 TAMO 感兴趣，它是另一个用于处理序列基序的 Python 库。它支持更多从头基序发现器，但它不是 Biopython 的一部分，并且在商业使用方面有一些限制。
基序对象
由于我们感兴趣的是基序分析，因此我们需要首先了解 Motif 对象。为此，我们需要导入 Bio.motifs 库
>>> from Bio import motifs
我们可以开始创建第一个基序对象。我们可以从基序实例列表中创建一个 Motif 对象，也可以通过解析来自基序数据库或基序发现软件的文件来获取 Motif 对象。
从实例创建基序
假设我们有以下 DNA 基序的实例
>>> from Bio.Seq import Seq
>>> instances = [
... Seq("TACAA"),
... Seq("TACGC"),
... Seq("TACAC"),
... Seq("TACCC"),
... Seq("AACCC"),
... Seq("AATGC"),
... Seq("AATGC"),
... ]
那么我们可以如下创建一个 Motif 对象
>>> m = motifs.create(instances)
创建此基序的实例存储在 .alignment 属性中
>>> print(m.alignment.sequences)
[Seq('TACAA'), Seq('TACGC'), Seq('TACAC'), Seq('TACCC'), Seq('AACCC'), Seq('AATGC'), Seq('AATGC')]
打印 Motif 对象显示了用来构造它的实例
>>> print(m)
TACAA
TACGC
TACAC
TACCC
AACCC
AATGC
AATGC
基序的长度定义为序列长度，所有实例的序列长度应该相同
>>> len(m)
5
Motif 对象具有一个属性 .counts，其中包含每个位置每个核苷酸的计数。打印此计数矩阵以易于阅读的格式显示它
>>> print(m.counts)
0 1 2 3 4
A: 3.00 7.00 0.00 2.00 1.00
C: 0.00 0.00 5.00 2.00 6.00
G: 0.00 0.00 0.00 3.00 0.00
T: 4.00 0.00 2.00 0.00 0.00
您可以将这些计数作为字典访问
>>> m.counts["A"]
[3.0, 7.0, 0.0, 2.0, 1.0]
但您也可以将其视为一个二维数组，其中第一个维度是核苷酸，第二个维度是位置
>>> m.counts["T", 0]
4.0
>>> m.counts["T", 2]
2.0
>>> m.counts["T", 3]
0.0
您还可以直接访问计数矩阵的列
>>> m.counts[:, 3]
{'A': 2.0, 'C': 2.0, 'T': 0.0, 'G': 3.0}
除了核苷酸本身之外，您还可以使用核苷酸在基序字母表中的索引
>>> m.alphabet
'ACGT'
>>> m.counts["A", :]
(3.0, 7.0, 0.0, 2.0, 1.0)
>>> m.counts[0, :]
(3.0, 7.0, 0.0, 2.0, 1.0)
获取一致性序列
基序的一致性序列定义为基序位置的字母序列，这些字母在 .counts 矩阵的相应列中获得最大值
>>> m.consensus
Seq('TACGC')
相反，反一致性序列对应于 .counts 矩阵的列中的最小值
>>> m.anticonsensus
Seq('CCATG')
请注意，如果某些列中多个核苷酸具有最大或最小计数，则一致性和反一致性序列的定义存在一些歧义。
对于 DNA 序列，您还可以要求一个简并一致性序列，其中使用模糊核苷酸来表示计数较高的多个核苷酸存在的位置
>>> m.degenerate_consensus
Seq('WACVC')
这里，W 和 R 遵循 IUPAC 核苷酸模糊代码：W 是 A 或 T，而 V 是 A、C 或 G [Cornish1985]。简并一致性序列是根据 Cavener [Cavener1987] 指定的规则构建的。
该 motif.counts.calculate_consensus 方法允许您详细指定如何计算一致性序列。此方法很大程度上遵循 EMBOSS 程序 cons 的约定，并采用以下参数
- 替换矩阵
比较序列时使用的评分矩阵。默认情况下，它为
None，在这种情况下，我们只需计算每个字母的频率。除了默认值之外，您还可以使用Bio.Align.substitution_matrices中提供的替换矩阵。常见的选择是 BLOSUM62（也称为 EBLOSUM62）用于蛋白质，以及 NUC.4.4（也称为 EDNAFULL）用于核苷酸。注意：当前，此方法尚未针对除默认值None以外的值实现。- 多数
正匹配数量的阈值，除以列中的总计数，要求达到一致。如果
substitution_matrix为None，则此参数也必须为None，并且会被忽略；否则会引发ValueError。如果substitution_matrix不是None，则多数的默认值为 0.5。- 同一性
同一性的数量，除以列中的总计数，要求定义一致性值。如果同一性的数量小于同一性乘以列中的总计数，则在一致性序列中使用未定义字符（对于核苷酸为
N，对于氨基酸序列为X）。如果identity为 1.0，则只有相同字母的列会对一致性做出贡献。默认值为零。- 设置大小写
正匹配的阈值，除以列中的总计数，超过该阈值的，一致性为大写，低于该阈值的，一致性为小写。默认情况下，它等于 0.5。
这是一个例子
>>> m.counts.calculate_consensus(identity=0.5, setcase=0.7)
'tACNC'
反向互补基序
我们可以通过调用 reverse_complement 方法来获得基序的反向互补
>>> r = m.reverse_complement()
>>> r.consensus
Seq('GCGTA')
>>> r.degenerate_consensus
Seq('GBGTW')
>>> print(r)
TTGTA
GCGTA
GTGTA
GGGTA
GGGTT
GCATT
GCATT
反向互补仅针对 DNA 基序定义。
切片基序
您可以切片基序以获得所选位置的新 Motif 对象
>>> m_sub = m[2:-1]
>>> print(m_sub)
CA
CG
CA
CC
CC
TG
TG
>>> m_sub.consensus
Seq('CG')
>>> m_sub.degenerate_consensus
Seq('CV')
相对熵
基序的第 \(j\) 列的相对熵（或 Kullback-Leibler 距离）\(H_j\) 定义如下 [Schneider1986] [Durbin1998]
其中
\(M\) – 字母表中的字母数量（由
len(m.alphabet)给出）；\(p_{ij}\) – 字母 \(i\) 的观察频率，归一化，在第 \(j\) 列（见下文）；
\(b_{i}\) – 字母 \(i\) 的背景概率（由
m.background[i]给出）。
观察频率 \(p_{ij}\) 的计算如下
其中
\(c_{ij}\) – 字母 \(i\) 在比对的第 \(j\) 列中出现的次数（由
m.counts[i, j]给出）；\(C_{j}\) – 第 \(j\) 列中的字母总数：\(C_{j} = \sum_{i=1}^{M} c_{ij}\)（由
sum(m.counts[:, j])给出）。\(k_i\) – 字母 \(i\) 的伪计数（由
m.pseudocounts[i]给出）。\(k\) – 总伪计数：\(k = \sum_{i=1}^{M} k_i\)（由
sum(m.pseudocounts.values())给出）。
根据这些定义，\(p_{ij}\) 和 \(b_{i}\) 都归一化为 1
如果背景分布是均匀的，则相对熵与信息内容相同。
可以使用 relative_entropy 属性获取基序 m 每列的相对熵。
>>> m.relative_entropy
array([1.01477186, 2. , 1.13687943, 0.44334329, 1.40832722])
这些值使用以 2 为底的对数计算，因此单位为比特。第二列（仅由 A 核苷酸组成）具有最高的相对熵；第四列（由 A、C 或 G 核苷酸组成）具有最低的相对熵）。基序的相对熵可以通过对各列求和来计算。
>>> print(f"Relative entropy is {sum(m.relative_entropy):0.5f}")
Relative entropy is 6.00332
创建序列标识
如果我们有互联网访问权限，我们可以创建一个 weblogo
>>> m.weblogo("mymotif.png")
我们应该将我们的标识保存为指定的 PNG 文件。
读取基序
手动从实例创建基序有点枯燥，因此拥有用于读取和写入基序的一些 I/O 函数非常有用。对于存储基序没有真正确立的标准，但有一些格式比其他格式使用更多。
JASPAR
最受欢迎的基序数据库之一是 JASPAR。除了基序序列信息外，JASPAR 数据库还为每个基序存储大量元信息。模块 Bio.motifs 包含一个专门的类 jaspar.Motif，其中此元信息表示为属性
matrix_id- 唯一的 JASPAR 基序 ID，例如 ’MA0004.1’name- TF 的名称，例如 ’Arnt’collection- 基序所属的 JASPAR 集合，例如 ’CORE’tf_class- 此 TF 的结构类别，例如 ’Zipper-Type’tf_family- 此 TF 所属的家族，例如 ’Helix-Loop-Helix’species- 此 TF 所属的物种，可能具有多个值，这些值指定为分类 ID，例如 10090tax_group- 此基序所属的分类学超群，例如 ’vertebrates’acc- TF 蛋白的登录号，例如 ’P53762’data_type- 用于构建此基序的数据类型，例如 ’SELEX’medline- 支持此基序的文献的 Pubmed ID，可能有多个值，例如 7592839pazar_id- PAZAR 数据库中 TF 的外部引用，例如 ’TF0000003’comment- 包含有关基序构建的说明的自由格式文本
jaspar.Motif 类继承自通用 Motif 类，因此提供了任何基序格式的所有功能——读取基序、写入基序、扫描序列以查找基序实例等。
JASPAR 以多种不同方式存储基序，包括三种不同的平面文件格式和 SQL 数据库。所有这些格式都有助于构建计数矩阵。但是，上面描述的元信息数量会因格式而异。
JASPAR sites 格式
这三种平面文件格式中的第一种包含一个实例列表。例如，这些是 JASPAR Arnt.sites 文件的开头和结尾行，显示了小鼠螺旋-环-螺旋转录因子 Arnt 的已知结合位点。
>MA0004 ARNT 1
CACGTGatgtcctc
>MA0004 ARNT 2
CACGTGggaggtac
>MA0004 ARNT 3
CACGTGccgcgcgc
...
>MA0004 ARNT 18
AACGTGacagccctcc
>MA0004 ARNT 19
AACGTGcacatcgtcc
>MA0004 ARNT 20
aggaatCGCGTGc
序列中用大写字母表示的部分是发现彼此对齐的基序实例。
我们可以从这些实例创建 Motif 对象，如下所示
>>> from Bio import motifs
>>> with open("Arnt.sites") as handle:
... arnt = motifs.read(handle, "sites")
...
创建此基序的实例存储在 .alignment 属性中
>>> print(arnt.alignment.sequences[:3])
[Seq('CACGTG'), Seq('CACGTG'), Seq('CACGTG')]
>>> for sequence in arnt.alignment.sequences:
... print(sequence)
...
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
AACGTG
AACGTG
AACGTG
AACGTG
CGCGTG
此基序的计数矩阵会自动从实例中计算出来
>>> print(arnt.counts)
0 1 2 3 4 5
A: 4.00 19.00 0.00 0.00 0.00 0.00
C: 16.00 0.00 20.00 0.00 0.00 0.00
G: 0.00 1.00 0.00 20.00 0.00 20.00
T: 0.00 0.00 0.00 0.00 20.00 0.00
此格式不存储任何元信息。
JASPAR pfm 格式
JASPAR 还直接将基序作为计数矩阵提供，而无需创建它的实例。此 pfm 格式仅存储单个基序的计数矩阵。例如，这是 JASPAR 文件 SRF.pfm，它包含人类 SRF 转录因子的计数矩阵
2 9 0 1 32 3 46 1 43 15 2 2
1 33 45 45 1 1 0 0 0 1 0 1
39 2 1 0 0 0 0 0 0 0 44 43
4 2 0 0 13 42 0 45 3 30 0 0
我们可以为这个计数矩阵创建基序，如下所示
>>> with open("SRF.pfm") as handle:
... srf = motifs.read(handle, "pfm")
...
>>> print(srf.counts)
0 1 2 3 4 5 6 7 8 9 10 11
A: 2.00 9.00 0.00 1.00 32.00 3.00 46.00 1.00 43.00 15.00 2.00 2.00
C: 1.00 33.00 45.00 45.00 1.00 1.00 0.00 0.00 0.00 1.00 0.00 1.00
G: 39.00 2.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 44.00 43.00
T: 4.00 2.00 0.00 0.00 13.00 42.00 0.00 45.00 3.00 30.00 0.00 0.00
由于此基序直接从计数矩阵创建，因此它没有与之关联的实例
>>> print(srf.alignment)
None
现在我们可以要求这两个基序的共有序列
>>> print(arnt.counts.consensus)
CACGTG
>>> print(srf.counts.consensus)
GCCCATATATGG
与实例文件一样，此格式不存储任何元信息。
JASPAR 格式 jaspar
jaspar 文件格式允许在单个文件中指定多个基序。在此格式中，每个基序记录都包含一个标题行，后面跟着四行定义计数矩阵。标题行以 > 字符开头（类似于 Fasta 文件格式），后面跟着唯一的 JASPAR 矩阵 ID 和 TF 名称。以下示例显示了一个 jaspar 格式的文件，其中包含三个基序 Arnt、RUNX1 和 MEF2A
>MA0004.1 Arnt
A [ 4 19 0 0 0 0 ]
C [16 0 20 0 0 0 ]
G [ 0 1 0 20 0 20 ]
T [ 0 0 0 0 20 0 ]
>MA0002.1 RUNX1
A [10 12 4 1 2 2 0 0 0 8 13 ]
C [ 2 2 7 1 0 8 0 0 1 2 2 ]
G [ 3 1 1 0 23 0 26 26 0 0 4 ]
T [11 11 14 24 1 16 0 0 25 16 7 ]
>MA0052.1 MEF2A
A [ 1 0 57 2 9 6 37 2 56 6 ]
C [50 0 1 1 0 0 0 0 0 0 ]
G [ 0 0 0 0 0 0 0 0 2 50 ]
T [ 7 58 0 55 49 52 21 56 0 2 ]
基序的读取方式如下
>>> fh = open("jaspar_motifs.txt")
>>> for m in motifs.parse(fh, "jaspar"):
... print(m)
...
TF name Arnt
Matrix ID MA0004.1
Matrix:
0 1 2 3 4 5
A: 4.00 19.00 0.00 0.00 0.00 0.00
C: 16.00 0.00 20.00 0.00 0.00 0.00
G: 0.00 1.00 0.00 20.00 0.00 20.00
T: 0.00 0.00 0.00 0.00 20.00 0.00
TF name RUNX1
Matrix ID MA0002.1
Matrix:
0 1 2 3 4 5 6 7 8 9 10
A: 10.00 12.00 4.00 1.00 2.00 2.00 0.00 0.00 0.00 8.00 13.00
C: 2.00 2.00 7.00 1.00 0.00 8.00 0.00 0.00 1.00 2.00 2.00
G: 3.00 1.00 1.00 0.00 23.00 0.00 26.00 26.00 0.00 0.00 4.00
T: 11.00 11.00 14.00 24.00 1.00 16.00 0.00 0.00 25.00 16.00 7.00
TF name MEF2A
Matrix ID MA0052.1
Matrix:
0 1 2 3 4 5 6 7 8 9
A: 1.00 0.00 57.00 2.00 9.00 6.00 37.00 2.00 56.00 6.00
C: 50.00 0.00 1.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00
G: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 50.00
T: 7.00 58.00 0.00 55.00 49.00 52.00 21.00 56.00 0.00 2.00
请注意，打印 JASPAR 基序将同时产生计数数据和可用的元信息。
访问 JASPAR 数据库
除了解析这些平面文件格式外，我们还可以从 JASPAR SQL 数据库中检索基序。与平面文件格式不同，JASPAR 数据库允许存储 JASPAR Motif 类中定义的所有可能的元信息。本文档中不涉及如何设置 JASPAR 数据库的说明（请参阅主要的 JASPAR 网站）。使用 Bio.motifs.jaspar.db 模块从 JASPAR 数据库中读取基序。首先使用 JASPAR5 类连接到 JASPAR 数据库，该类模拟最新的 JASPAR 架构
>>> from Bio.motifs.jaspar.db import JASPAR5
>>>
>>> JASPAR_DB_HOST = "yourhostname" # fill in these values
>>> JASPAR_DB_NAME = "yourdatabase"
>>> JASPAR_DB_USER = "yourusername"
>>> JASPAR_DB_PASS = "yourpassword"
>>>
>>> jdb = JASPAR5(
... host=JASPAR_DB_HOST,
... name=JASPAR_DB_NAME,
... user=JASPAR_DB_USER,
... password=JASPAR_DB_PASS,
... )
现在，我们可以使用 fetch_motif_by_id 方法通过其唯一的 JASPAR ID 获取单个基序。请注意，JASPAR ID 由一个基本 ID 和一个小数点分隔的版本号组成，例如 ’MA0004.1’。fetch_motif_by_id 方法允许您使用完全指定的 ID 或仅使用基本 ID。如果只提供基本 ID，则将返回基序的最新版本。
>>> arnt = jdb.fetch_motif_by_id("MA0004")
打印基序表明 JASPAR SQL 数据库存储的元信息远多于平面文件
>>> print(arnt)
TF name Arnt
Matrix ID MA0004.1
Collection CORE
TF class Zipper-Type
TF family Helix-Loop-Helix
Species 10090
Taxonomic group vertebrates
Accession ['P53762']
Data type used SELEX
Medline 7592839
PAZAR ID TF0000003
Comments -
Matrix:
0 1 2 3 4 5
A: 4.00 19.00 0.00 0.00 0.00 0.00
C: 16.00 0.00 20.00 0.00 0.00 0.00
G: 0.00 1.00 0.00 20.00 0.00 20.00
T: 0.00 0.00 0.00 0.00 20.00 0.00
我们还可以按名称获取基序。名称必须完全匹配（目前不支持部分匹配或数据库通配符）。请注意，由于名称不能保证唯一，因此 fetch_motifs_by_name 方法实际上返回一个列表。
>>> motifs = jdb.fetch_motifs_by_name("Arnt")
>>> print(motifs[0])
TF name Arnt
Matrix ID MA0004.1
Collection CORE
TF class Zipper-Type
TF family Helix-Loop-Helix
Species 10090
Taxonomic group vertebrates
Accession ['P53762']
Data type used SELEX
Medline 7592839
PAZAR ID TF0000003
Comments -
Matrix:
0 1 2 3 4 5
A: 4.00 19.00 0.00 0.00 0.00 0.00
C: 16.00 0.00 20.00 0.00 0.00 0.00
G: 0.00 1.00 0.00 20.00 0.00 20.00
T: 0.00 0.00 0.00 0.00 20.00 0.00
fetch_motifs 方法允许您获取与指定条件集匹配的基序。这些条件包括上面描述的任何元信息，以及某些矩阵属性，例如最小信息含量（以下示例中的 min_ic）、矩阵的最小长度或用于构建矩阵的最小站点数量。只有通过所有指定条件的基序才会返回。请注意，与允许使用多个值的元信息相对应的选择条件可以指定为单个值或值列表，例如以下示例中的 tax_group 和 tf_family。
>>> motifs = jdb.fetch_motifs(
... collection="CORE",
... tax_group=["vertebrates", "insects"],
... tf_class="Winged Helix-Turn-Helix",
... tf_family=["Forkhead", "Ets"],
... min_ic=12,
... )
>>> for motif in motifs:
... pass # do something with the motif
...
与 Perl TFBS 模块的兼容性
需要注意的重要一点是，JASPAR Motif 类旨在与流行的 Perl TFBS 模块 兼容。因此，关于默认背景和伪计数选择以及信息含量计算方式和序列搜索实例方式的某些细节基于此兼容性标准。这些选择将在下面的特定部分中说明。
- 背景选择Perl
TFBS模块似乎允许选择自定义背景概率（尽管文档指出假设统一背景）。但是，默认情况下使用统一背景。因此，建议您使用统一背景来计算位置特异性评分矩阵 (PSSM)。这是使用 Biopythonmotifs模块时的默认设置。 - 伪计数的选择默认情况下，Perl
TFBS模块使用等于 \(\sqrt{N} * \textrm{bg}[\textrm{nucleotide}]\) 的伪计数，其中 \(N\) 表示用于构建矩阵的序列总数。要应用相同的伪计数公式，请使用jaspar.calculate\_pseudcounts()函数设置基序pseudocounts属性>>> motif.pseudocounts = motifs.jaspar.calculate_pseudocounts(motif)
请注意，计数矩阵可能具有构成各列的序列数量不相等。伪计数计算使用构成矩阵的序列的平均数量。但是，当对计数矩阵调用
normalize时，列中的每个计数值都会除以构成该特定列的序列总数，而不是除以序列的平均数量。这与 PerlTFBS模块不同，因为规范化不是作为单独的步骤完成的，因此在整个 pssm 计算过程中使用序列的平均数量。因此，对于具有不等列计数的矩阵，motifs模块计算的 PSSM 将与 PerlTFBS模块计算的 pssm 有所不同。 - 矩阵信息含量的计算使用
PositionSpecificScoringMatrix类的mean方法计算矩阵的信息含量 (IC) 或特异性。但是需要注意的是，在 PerlTFBS模块中，默认行为是在不先应用伪计数的情况下计算 IC，即使默认情况下，PSSMs 是使用上面描述的伪计数计算的。 - 搜索实例使用 Perl
TFBS基序搜索实例通常使用相对分数阈值，即 0 到 1 之间的分数。为了计算对应于相对分数的绝对 PSSM 分数，可以使用以下公式：>>> abs_score = (pssm.max - pssm.min) * rel_score + pssm.min
要将实例的绝对分数转换回相对分数，可以使用以下公式：
>>> rel_score = (abs_score - pssm.min) / (pssm.max - pssm.min)
例如，使用之前的 Arnt 基序，让我们搜索一个相对分数阈值为 0.8 的序列。
>>> test_seq = Seq("TAAGCGTGCACGCGCAACACGTGCATTA") >>> arnt.pseudocounts = motifs.jaspar.calculate_pseudocounts(arnt) >>> pssm = arnt.pssm >>> max_score = pssm.max >>> min_score = pssm.min >>> abs_score_threshold = (max_score - min_score) * 0.8 + min_score >>> for pos, score in pssm.search(test_seq, threshold=abs_score_threshold): ... rel_score = (score - min_score) / (max_score - min_score) ... print(f"Position {pos}: score = {score:5.3f}, rel. score = {rel_score:5.3f}") ... Position 2: score = 5.362, rel. score = 0.801 Position 8: score = 6.112, rel. score = 0.831 Position -20: score = 7.103, rel. score = 0.870 Position 17: score = 10.351, rel. score = 1.000 Position -11: score = 10.351, rel. score = 1.000
MEME
MEME [Bailey1994] 是一款用于在相关 DNA 或蛋白质序列组中发现基序的工具。它将一组 DNA 或蛋白质序列作为输入，并输出尽可能多的基序。因此，与 JASPAR 文件相比，MEME 输出文件通常包含多个基序。这是一个例子。
在 MEME 生成的输出文件顶部显示有关 MEME 和所用 MEME 版本的一些背景信息
********************************************************************************
MEME - Motif discovery tool
********************************************************************************
MEME version 3.0 (Release date: 2004/08/18 09:07:01)
...
向下，训练序列的输入集会重新出现
********************************************************************************
TRAINING SET
********************************************************************************
DATAFILE= INO_up800.s
ALPHABET= ACGT
Sequence name Weight Length Sequence name Weight Length
------------- ------ ------ ------------- ------ ------
CHO1 1.0000 800 CHO2 1.0000 800
FAS1 1.0000 800 FAS2 1.0000 800
ACC1 1.0000 800 INO1 1.0000 800
OPI3 1.0000 800
********************************************************************************
以及所用的确切命令行
********************************************************************************
COMMAND LINE SUMMARY
********************************************************************************
This information can also be useful in the event you wish to report a
problem with the MEME software.
command: meme -mod oops -dna -revcomp -nmotifs 2 -bfile yeast.nc.6.freq INO_up800.s
...
接下来是有关找到的每个基序的详细信息
********************************************************************************
MOTIF 1 width = 12 sites = 7 llr = 95 E-value = 2.0e-001
********************************************************************************
--------------------------------------------------------------------------------
Motif 1 Description
--------------------------------------------------------------------------------
Simplified A :::9:a::::3:
pos.-specific C ::a:9:11691a
probability G ::::1::94:4:
matrix T aa:1::9::11:
要解析此文件（存储为 meme.dna.oops.txt），请使用
>>> with open("meme.INO_up800.classic.oops.xml") as handle:
... record = motifs.parse(handle, "meme")
...
命令 motifs.parse 会直接读取完整文件，因此您可以在调用 motifs.parse 后关闭文件。头信息存储在属性中
>>> record.version
'5.0.1'
>>> record.datafile
'common/INO_up800.s'
>>> record.command
'meme common/INO_up800.s -oc results/meme10 -mod oops -dna -revcomp -bfile common/yeast.nc.6.freq -nmotifs 2 -objfun classic -minw 8 -nostatus '
>>> record.alphabet
'ACGT'
>>> record.sequences
['sequence_0', 'sequence_1', 'sequence_2', 'sequence_3', 'sequence_4', 'sequence_5', 'sequence_6']
记录是 Bio.motifs.meme.Record 类的一个对象。该类继承自 list，您可以将 record 视为一个 Motif 对象列表
>>> len(record)
2
>>> motif = record[0]
>>> print(motif.consensus)
GCGGCATGTGAAA
>>> print(motif.degenerate_consensus)
GSKGCATGTGAAA
除了这些通用基序属性之外，每个基序还存储由 MEME 计算的特定信息。例如，
>>> motif.num_occurrences
7
>>> motif.length
13
>>> evalue = motif.evalue
>>> print("%3.1g" % evalue)
0.2
>>> motif.name
'GSKGCATGTGAAA'
>>> motif.id
'motif_1'
除了像上面那样使用记录中的索引之外，您还可以通过其名称找到它
>>> motif = record["GSKGCATGTGAAA"]
每个基序都有一个属性 .alignment，其中包含找到基序的序列比对，提供有关每个序列的一些信息
>>> len(motif.alignment)
7
>>> motif.alignment.sequences[0]
Instance('GCGGCATGTGAAA')
>>> motif.alignment.sequences[0].motif_name
'GSKGCATGTGAAA'
>>> motif.alignment.sequences[0].sequence_name
'INO1'
>>> motif.alignment.sequences[0].sequence_id
'sequence_5'
>>> motif.alignment.sequences[0].start
620
>>> motif.alignment.sequences[0].strand
'+'
>>> motif.alignment.sequences[0].length
13
>>> pvalue = motif.alignment.sequences[0].pvalue
>>> print("%5.3g" % pvalue)
1.21e-08
MAST
TRANSFAC
TRANSFAC 是一个手动整理的转录因子数据库，其中包含它们的基因组结合位点和 DNA 结合谱 [Matys2003]。虽然 TRANSFAC 数据库中使用的文件格式现在也被其他人使用，但我们将将其称为 TRANSFAC 文件格式。
TRANSFAC 格式的最小文件如下所示
ID motif1
P0 A C G T
01 1 2 2 0 S
02 2 1 2 0 R
03 3 0 1 1 A
04 0 5 0 0 C
05 5 0 0 0 A
06 0 0 4 1 G
07 0 1 4 0 G
08 0 0 0 5 T
09 0 0 5 0 G
10 0 1 2 2 K
11 0 2 0 3 Y
12 1 0 3 1 G
//
此文件显示了 12 个核苷酸的基序 motif1 的频率矩阵。通常，一个 TRANSFAC 格式的文件可以包含多个基序。例如，这是示例 TRANSFAC 文件 transfac.dat 的内容
VV EXAMPLE January 15, 2013
XX
//
ID motif1
P0 A C G T
01 1 2 2 0 S
02 2 1 2 0 R
03 3 0 1 1 A
...
11 0 2 0 3 Y
12 1 0 3 1 G
//
ID motif2
P0 A C G T
01 2 1 2 0 R
02 1 2 2 0 S
...
09 0 0 0 5 T
10 0 2 0 3 Y
//
要解析 TRANSFAC 文件，请使用
>>> with open("transfac.dat") as handle:
... record = motifs.parse(handle, "TRANSFAC")
...
如果检测到文件内容与 TRANSFAC 文件格式之间的任何差异，则会引发 ValueError。请注意，您可能会遇到不严格遵循 TRANSFAC 格式的文件。例如，列之间的空格数可能不同，或者可以使用制表符而不是空格。使用 strict=False 启用解析此类文件，而不会引发 ValueError
>>> record = motifs.parse(handle, "TRANSFAC", strict=False)
解析不符合规范的文件时，建议检查 motif.parse 返回的记录，以确保它与文件内容一致。
总体版本号（如果有）存储为 record.version
>>> record.version
'EXAMPLE January 15, 2013'
中的每个基序 record 都是 Bio.motifs.transfac.Motif 类的实例，它既继承自 Bio.motifs.Motif 类，也继承自 Python 字典。字典使用两位字母键来存储有关基序的任何附加信息
>>> motif = record[0]
>>> motif.degenerate_consensus # Using the Bio.motifs.Motif property
Seq('SRACAGGTGKYG')
>>> motif["ID"] # Using motif as a dictionary
'motif1'
TRANSFAC 文件通常比这个例子复杂得多，包含大量有关基序的附加信息。表 TRANSFAC 文件中常见的字段 列出了 TRANSFAC 文件中常见的两位字母字段代码
|
登录号 |
|
登录号，次要 |
|
统计基础 |
|
结合因子 |
|
矩阵背后的因子结合位点 |
|
评论 |
|
版权声明 |
|
简短的因子描述 |
|
外部数据库 |
|
创建/更新日期 |
|
亚家族 |
|
超家族 |
|
标识符 |
|
结合因子的名称 |
|
分类学分类 |
|
物种/分类单元 |
|
旧版本 |
|
首选版本 |
|
类型 |
|
空行；这些不会存储在记录中。 |
每个基序也都有一个属性 .references，其中包含与基序关联的引用，使用以下两位字母键
|
参考文献号 |
|
参考文献作者 |
|
参考文献数据 |
|
参考文献标题 |
|
PubMed ID |
打印基序会以其本机 TRANSFAC 格式输出它们
>>> print(record)
VV EXAMPLE January 15, 2013
XX
//
ID motif1
XX
P0 A C G T
01 1 2 2 0 S
02 2 1 2 0 R
03 3 0 1 1 A
04 0 5 0 0 C
05 5 0 0 0 A
06 0 0 4 1 G
07 0 1 4 0 G
08 0 0 0 5 T
09 0 0 5 0 G
10 0 1 2 2 K
11 0 2 0 3 Y
12 1 0 3 1 G
XX
//
ID motif2
XX
P0 A C G T
01 2 1 2 0 R
02 1 2 2 0 S
03 0 5 0 0 C
04 3 0 1 1 A
05 0 0 4 1 G
06 5 0 0 0 A
07 0 1 4 0 G
08 0 0 5 0 G
09 0 0 0 5 T
10 0 2 0 3 Y
XX
//
您可以通过将此输出捕获到一个字符串中并将其保存到一个文件中来以 TRANSFAC 格式导出基序
>>> text = str(record)
>>> with open("mytransfacfile.dat", "w") as out_handle:
... out_handle.write(text)
...
写入基序
说到导出，让我们看一下导出函数的总体情况。我们可以使用内置函数 format 以简单的 JASPAR pfm 格式写入基序
>>> print(format(arnt, "pfm"))
4.00 19.00 0.00 0.00 0.00 0.00
16.00 0.00 20.00 0.00 0.00 0.00
0.00 1.00 0.00 20.00 0.00 20.00
0.00 0.00 0.00 0.00 20.00 0.00
同样，我们可以使用 format 以 JASPAR jaspar 格式写入基序
>>> print(format(arnt, "jaspar"))
>MA0004.1 Arnt
A [ 4.00 19.00 0.00 0.00 0.00 0.00]
C [ 16.00 0.00 20.00 0.00 0.00 0.00]
G [ 0.00 1.00 0.00 20.00 0.00 20.00]
T [ 0.00 0.00 0.00 0.00 20.00 0.00]
要以类似 TRANSFAC 的矩阵格式写入基序，请使用
>>> print(format(m, "transfac"))
P0 A C G T
01 3 0 0 4 W
02 7 0 0 0 A
03 0 5 0 2 C
04 2 2 3 0 V
05 1 6 0 0 C
XX
//
要写入多个基序，可以使用 motifs.write。无论基序来自 TRANSFAC 文件还是其他来源，都可以使用此函数。例如，
>>> two_motifs = [arnt, srf]
>>> print(motifs.write(two_motifs, "transfac"))
P0 A C G T
01 4 16 0 0 C
02 19 0 1 0 A
03 0 20 0 0 C
04 0 0 20 0 G
05 0 0 0 20 T
06 0 0 20 0 G
XX
//
P0 A C G T
01 2 1 39 4 G
02 9 33 2 2 C
03 0 45 1 0 C
04 1 45 0 0 C
05 32 1 0 13 A
06 3 1 0 42 T
07 46 0 0 0 A
08 1 0 0 45 T
09 43 0 0 3 A
10 15 1 0 30 W
11 2 0 44 0 G
12 2 1 43 0 G
XX
//
或者，要以 jaspar 格式写入多个基序
>>> two_motifs = [arnt, mef2a]
>>> print(motifs.write(two_motifs, "jaspar"))
>MA0004.1 Arnt
A [ 4.00 19.00 0.00 0.00 0.00 0.00]
C [ 16.00 0.00 20.00 0.00 0.00 0.00]
G [ 0.00 1.00 0.00 20.00 0.00 20.00]
T [ 0.00 0.00 0.00 0.00 20.00 0.00]
>MA0052.1 MEF2A
A [ 1.00 0.00 57.00 2.00 9.00 6.00 37.00 2.00 56.00 6.00]
C [ 50.00 0.00 1.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00]
G [ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 50.00]
T [ 7.00 58.00 0.00 55.00 49.00 52.00 21.00 56.00 0.00 2.00]
位置权重矩阵
Motif 对象的 .counts 属性显示每个核苷酸在比对中每个位置出现的频率。我们可以通过除以比对中的实例数量来对该矩阵进行归一化，从而得到每个核苷酸在比对中每个位置的概率。我们将这些概率称为位置权重矩阵。但是，请注意，在文献中，此术语也可能用于指代位置特定评分矩阵，我们将在下面讨论。
通常，在归一化之前会将伪计数添加到每个位置。这避免了位置权重矩阵过度拟合比对中数量有限的基序实例，并且还可以防止概率变为零。要将固定伪计数添加到所有位置的所有核苷酸，请为 pseudocounts 参数指定一个数字
>>> pwm = m.counts.normalize(pseudocounts=0.5)
>>> print(pwm)
0 1 2 3 4
A: 0.39 0.83 0.06 0.28 0.17
C: 0.06 0.06 0.61 0.28 0.72
G: 0.06 0.06 0.06 0.39 0.06
T: 0.50 0.06 0.28 0.06 0.06
或者，pseudocounts 可以是一个字典，指定每个核苷酸的伪计数。例如，由于人类基因组的 GC 含量约为 40%，您可能希望相应地选择伪计数
>>> pwm = m.counts.normalize(pseudocounts={"A": 0.6, "C": 0.4, "G": 0.4, "T": 0.6})
>>> print(pwm)
0 1 2 3 4
A: 0.40 0.84 0.07 0.29 0.18
C: 0.04 0.04 0.60 0.27 0.71
G: 0.04 0.04 0.04 0.38 0.04
T: 0.51 0.07 0.29 0.07 0.07
位置权重矩阵有其自己的方法来计算共识、反共识和简并共识序列
>>> pwm.consensus
Seq('TACGC')
>>> pwm.anticonsensus
Seq('CCGTG')
>>> pwm.degenerate_consensus
Seq('WACNC')
请注意，由于伪计数，从位置权重矩阵计算出的简并共识序列与从基序中的实例计算出的简并共识序列略有不同
>>> m.degenerate_consensus
Seq('WACVC')
位置权重矩阵的反向互补可以从 pwm 直接计算
>>> rpwm = pwm.reverse_complement()
>>> print(rpwm)
0 1 2 3 4
A: 0.07 0.07 0.29 0.07 0.51
C: 0.04 0.38 0.04 0.04 0.04
G: 0.71 0.27 0.60 0.04 0.04
T: 0.18 0.29 0.07 0.84 0.40
位置特定评分矩阵
使用背景分布和添加了伪计数的 PWM，很容易计算对数几率比，告诉我们特定符号来自基序而不是来自背景的对数几率。我们可以使用位置权重矩阵上的 .log_odds() 方法
>>> pssm = pwm.log_odds()
>>> print(pssm)
0 1 2 3 4
A: 0.68 1.76 -1.91 0.21 -0.49
C: -2.49 -2.49 1.26 0.09 1.51
G: -2.49 -2.49 -2.49 0.60 -2.49
T: 1.03 -1.91 0.21 -1.91 -1.91
这里我们可以看到基序中比背景中更频繁的符号的正值，以及背景中比基序中更频繁的符号的负值。\(0.0\) 表示在背景和基序中看到符号的可能性相同。
这假设 A、C、G 和 T 在背景中出现的可能性相同。要针对具有 A、C、G、T 不相等概率的背景计算位置特定评分矩阵，请使用 background 参数。例如，针对具有 40% GC 含量的背景，请使用
>>> background = {"A": 0.3, "C": 0.2, "G": 0.2, "T": 0.3}
>>> pssm = pwm.log_odds(background)
>>> print(pssm)
0 1 2 3 4
A: 0.42 1.49 -2.17 -0.05 -0.75
C: -2.17 -2.17 1.58 0.42 1.83
G: -2.17 -2.17 -2.17 0.92 -2.17
T: 0.77 -2.17 -0.05 -2.17 -2.17
PSSM 可获得的最大和最小分数分别存储在 .max 和 .min 属性中
>>> print("%4.2f" % pssm.max)
6.59
>>> print("%4.2f" % pssm.min)
-10.85
通过 .mean 和 .std 方法计算针对特定背景的 PSSM 分数的平均值和标准差。
>>> mean = pssm.mean(background)
>>> std = pssm.std(background)
>>> print("mean = %0.2f, standard deviation = %0.2f" % (mean, std))
mean = 3.21, standard deviation = 2.59
如果没有指定 background，则使用统一背景。平均值等于第 相对熵 节中描述的 Kullback-Leibler 散度或相对熵。
可以将 .reverse_complement、.consensus、.anticonsensus 和 .degenerate_consensus 方法直接应用于 PSSM 对象。
搜索实例
基序最常见的用途是在某些序列中找到其实例。在本节的示例中，我们将使用类似这样的一个人工序列
>>> test_seq = Seq("TACACTGCATTACAACCCAAGCATTA")
>>> len(test_seq)
26
搜索完全匹配
寻找实例的最简单方法是寻找基序真实实例的完全匹配
>>> for pos, seq in test_seq.search(m.alignment):
... print("%i %s" % (pos, seq))
...
0 TACAC
10 TACAA
13 AACCC
我们可以对反向互补执行相同的操作（以查找互补链上的实例）
>>> for pos, seq in test_seq.search(r.alignment):
... print("%i %s" % (pos, seq))
...
6 GCATT
20 GCATT
使用 PSSM 分数搜索匹配
寻找产生高对数几率分数的基序位置也很容易
>>> for position, score in pssm.search(test_seq, threshold=3.0):
... print("Position %d: score = %5.3f" % (position, score))
...
Position 0: score = 5.622
Position -20: score = 4.601
Position 10: score = 3.037
Position 13: score = 5.738
Position -6: score = 4.601
负位置指的是在测试序列的反向链上找到的基序实例，并遵循 Python 中负索引的惯例。因此，位于 pos 的基序实例位于 test_seq[pos:pos+len(m)]，无论 pos 的值是正数还是负数。
您可能会注意到阈值参数，这里任意设置为 \(3.0\)。这是以 \(log_2\) 为单位的，因此我们现在只关注那些在基序模型下出现的概率是背景概率的八倍的词语。默认阈值为 \(0.0\)，它选择所有看起来更像基序而不是背景的词语。
您还可以计算序列中所有位置的分数
>>> pssm.calculate(test_seq)
array([ 5.62230396, -5.6796999 , -3.43177247, 0.93827754,
-6.84962511, -2.04066086, -10.84962463, -3.65614533,
-0.03370807, -3.91102552, 3.03734159, -2.14918518,
-0.6016975 , 5.7381525 , -0.50977498, -3.56422281,
-8.73414803, -0.09919716, -0.6016975 , -2.39429784,
-10.84962463, -3.65614533], dtype=float32)
总的来说，这是计算 PSSM 分数最快的方法。由 pssm.calculate 返回的分数仅针对正向链。要获取反向链上的分数，您可以取 PSSM 的反向互补序列。
>>> rpssm = pssm.reverse_complement()
>>> rpssm.calculate(test_seq)
array([ -9.43458748, -3.06172252, -7.18665981, -7.76216221,
-2.04066086, -4.26466274, 4.60124254, -4.2480607 ,
-8.73414803, -2.26503372, -6.49598789, -5.64668512,
-8.73414803, -10.84962463, -4.82356262, -4.82356262,
-5.64668512, -8.73414803, -4.15613794, -5.6796999 ,
4.60124254, -4.2480607 ], dtype=float32)
选择分数阈值
如果您想要使用一种不太任意的方法来选择阈值，您可以探索 PSSM 分数的分布。由于分数分布的空间随着基序长度呈指数级增长，因此我们使用给定精度的近似值来保持计算成本的可控。
>>> distribution = pssm.distribution(background=background, precision=10**4)
可以使用 distribution 对象来确定许多不同的阈值。我们可以指定所需的假阳性率（在背景生成的序列中“找到”基序实例的概率）
>>> threshold = distribution.threshold_fpr(0.01)
>>> print("%5.3f" % threshold)
4.009
或者假阴性率（“未找到”从基序生成的实例的概率）
>>> threshold = distribution.threshold_fnr(0.1)
>>> print("%5.3f" % threshold)
-0.510
或者一个（近似）满足假阳性率和假阴性率之间某种关系的阈值（\(\frac{\textrm{fnr}}{\textrm{fpr}}\simeq t\)）
>>> threshold = distribution.threshold_balanced(1000)
>>> print("%5.3f" % threshold)
6.241
或者一个（大致）满足假阳性率的 \(-log\) 和信息含量（如 Hertz 和 Stormo 在 patser 软件中使用的那样）相等的阈值。
>>> threshold = distribution.threshold_patser()
>>> print("%5.3f" % threshold)
0.346
例如，在我们的基序的情况下，您可以获得一个阈值，该阈值可以为您提供与搜索具有平衡阈值（速率为 \(1000\)）的实例相同的精确结果（对于此序列）。
>>> threshold = distribution.threshold_fpr(0.01)
>>> print("%5.3f" % threshold)
4.009
>>> for position, score in pssm.search(test_seq, threshold=threshold):
... print("Position %d: score = %5.3f" % (position, score))
...
Position 0: score = 5.622
Position -20: score = 4.601
Position 13: score = 5.738
Position -6: score = 4.601
每个基序对象都与一个相关的特定位置评分矩阵
为了方便使用 PSSM 搜索潜在的 TFBS，位置权重矩阵和特定位置评分矩阵都与每个基序相关联。以 Arnt 基序为例
>>> from Bio import motifs
>>> with open("Arnt.sites") as handle:
... motif = motifs.read(handle, "sites")
...
>>> print(motif.counts)
0 1 2 3 4 5
A: 4.00 19.00 0.00 0.00 0.00 0.00
C: 16.00 0.00 20.00 0.00 0.00 0.00
G: 0.00 1.00 0.00 20.00 0.00 20.00
T: 0.00 0.00 0.00 0.00 20.00 0.00
>>> print(motif.pwm)
0 1 2 3 4 5
A: 0.20 0.95 0.00 0.00 0.00 0.00
C: 0.80 0.00 1.00 0.00 0.00 0.00
G: 0.00 0.05 0.00 1.00 0.00 1.00
T: 0.00 0.00 0.00 0.00 1.00 0.00
>>> print(motif.pssm)
0 1 2 3 4 5
A: -0.32 1.93 -inf -inf -inf -inf
C: 1.68 -inf 2.00 -inf -inf -inf
G: -inf -2.32 -inf 2.00 -inf 2.00
T: -inf -inf -inf -inf 2.00 -inf
负无穷大在这里出现是因为频率矩阵中的对应条目为 0，并且默认情况下我们使用零伪计数。
>>> for letter in "ACGT":
... print("%s: %4.2f" % (letter, motif.pseudocounts[letter]))
...
A: 0.00
C: 0.00
G: 0.00
T: 0.00
如果您更改 .pseudocounts 属性，则位置频率矩阵和特定位置评分矩阵将自动重新计算。
>>> motif.pseudocounts = 3.0
>>> for letter in "ACGT":
... print("%s: %4.2f" % (letter, motif.pseudocounts[letter]))
...
A: 3.00
C: 3.00
G: 3.00
T: 3.00
>>> print(motif.pwm)
0 1 2 3 4 5
A: 0.22 0.69 0.09 0.09 0.09 0.09
C: 0.59 0.09 0.72 0.09 0.09 0.09
G: 0.09 0.12 0.09 0.72 0.09 0.72
T: 0.09 0.09 0.09 0.09 0.72 0.09
>>> print(motif.pssm)
0 1 2 3 4 5
A: -0.19 1.46 -1.42 -1.42 -1.42 -1.42
C: 1.25 -1.42 1.52 -1.42 -1.42 -1.42
G: -1.42 -1.00 -1.42 1.52 -1.42 1.52
T: -1.42 -1.42 -1.42 -1.42 1.52 -1.42
如果您想要对它们使用不同的伪计数，也可以将 .pseudocounts 设置为四个核苷酸的字典。将 motif.pseudocounts 设置为 None 将将其重置为其默认值零。
特定位置评分矩阵取决于背景分布，默认情况下为均匀分布。
>>> for letter in "ACGT":
... print("%s: %4.2f" % (letter, motif.background[letter]))
...
A: 0.25
C: 0.25
G: 0.25
T: 0.25
同样，如果您修改背景分布，则特定位置评分矩阵将重新计算。
>>> motif.background = {"A": 0.2, "C": 0.3, "G": 0.3, "T": 0.2}
>>> print(motif.pssm)
0 1 2 3 4 5
A: 0.13 1.78 -1.09 -1.09 -1.09 -1.09
C: 0.98 -1.68 1.26 -1.68 -1.68 -1.68
G: -1.68 -1.26 -1.68 1.26 -1.68 1.26
T: -1.09 -1.09 -1.09 -1.09 1.85 -1.09
将 motif.background 设置为 None 将将其重置为均匀分布。
>>> motif.background = None
>>> for letter in "ACGT":
... print("%s: %4.2f" % (letter, motif.background[letter]))
...
A: 0.25
C: 0.25
G: 0.25
T: 0.25
如果您将 motif.background 设置为单个值，它将被解释为 GC 含量。
>>> motif.background = 0.8
>>> for letter in "ACGT":
... print("%s: %4.2f" % (letter, motif.background[letter]))
...
A: 0.10
C: 0.40
G: 0.40
T: 0.10
请注意，您现在可以计算 PSSM 分数在计算其所针对的背景的平均值。
>>> print("%f" % motif.pssm.mean(motif.background))
4.703928
以及它的标准差。
>>> print("%f" % motif.pssm.std(motif.background))
3.290900
以及它的分布。
>>> distribution = motif.pssm.distribution(background=motif.background)
>>> threshold = distribution.threshold_fpr(0.01)
>>> print("%f" % threshold)
3.854375
请注意，每次您调用 motif.pwm 或 motif.pssm 时，位置权重矩阵和特定位置评分矩阵都会重新计算。如果速度是一个问题，并且您想重复使用 PWM 或 PSSM，您可以将它们保存为一个变量，如
>>> pssm = motif.pssm
比较基序
一旦我们拥有多个基序，我们可能想要比较它们。
在我们开始比较基序之前，我应该指出基序边界通常是相当任意的。这意味着我们经常需要比较不同长度的基序，因此比较需要包含某种形式的对齐。这意味着我们必须考虑两件事
基序的对齐
用于比较对齐基序的某个函数
为了对齐基序，我们使用 PSSM 的无间隙对齐，并将零替换为矩阵开头和结尾处任何缺失的列。这意味着实际上我们正在使用背景分布来表示 PSSM 中缺失的列。然后，距离函数返回基序之间的最小距离，以及它们对齐中的对应偏移量。
举个例子，让我们首先加载另一个基序，它类似于我们的测试基序 m
>>> with open("REB1.pfm") as handle:
... m_reb1 = motifs.read(handle, "pfm")
...
>>> m_reb1.consensus
Seq('GTTACCCGG')
>>> print(m_reb1.counts)
0 1 2 3 4 5 6 7 8
A: 30.00 0.00 0.00 100.00 0.00 0.00 0.00 0.00 15.00
C: 10.00 0.00 0.00 0.00 100.00 100.00 100.00 0.00 15.00
G: 50.00 0.00 0.00 0.00 0.00 0.00 0.00 60.00 55.00
T: 10.00 100.00 100.00 0.00 0.00 0.00 0.00 40.00 15.00
为了使基序具有可比性，我们选择与我们的基序 m 相同的伪计数和背景分布值。
>>> m_reb1.pseudocounts = {"A": 0.6, "C": 0.4, "G": 0.4, "T": 0.6}
>>> m_reb1.background = {"A": 0.3, "C": 0.2, "G": 0.2, "T": 0.3}
>>> pssm_reb1 = m_reb1.pssm
>>> print(pssm_reb1)
0 1 2 3 4 5 6 7 8
A: 0.00 -5.67 -5.67 1.72 -5.67 -5.67 -5.67 -5.67 -0.97
C: -0.97 -5.67 -5.67 -5.67 2.30 2.30 2.30 -5.67 -0.41
G: 1.30 -5.67 -5.67 -5.67 -5.67 -5.67 -5.67 1.57 1.44
T: -1.53 1.72 1.72 -5.67 -5.67 -5.67 -5.67 0.41 -0.97
我们将使用皮尔逊相关性来比较这些基序。由于我们希望它类似于距离度量，因此我们实际上取 \(1-r\)，其中 \(r\) 是皮尔逊相关系数 (PCC)
>>> distance, offset = pssm.dist_pearson(pssm_reb1)
>>> print("distance = %5.3g" % distance)
distance = 0.239
>>> print(offset)
-2
这意味着基序 m 和 m_reb1 之间的最佳 PCC 是通过以下对齐获得的
m: bbTACGCbb
m_reb1: GTTACCCGG
其中 b 代表背景分布。PCC 本身大约为 \(1-0.239=0.761\)。
从头 基序发现
目前，Biopython 仅对从头 基序发现提供有限的支持。具体来说，我们支持运行 xxmotif，以及解析 MEME。由于基序发现工具的数量正在快速增长，因此欢迎贡献新的解析器。
MEME
假设您已经使用您最喜欢的参数在您选择的序列上运行了 MEME，并将输出保存到文件 meme.out 中。您可以通过运行以下代码段来检索 MEME 报告的基序。
>>> from Bio import motifs
>>> with open("meme.psp_test.classic.zoops.xml") as handle:
... motifsM = motifs.parse(handle, "meme")
...
>>> motifsM
[<Bio.motifs.meme.Motif object at 0xc356b0>]
除了最想要的基序列表外，结果对象还包含更多有用的信息，可以通过具有自解释名称的属性访问。
.alphabet.datafile.sequences.version.command
由 MEME 解析器返回的基序可以像常规的基序对象（带有实例）一样对待，它们还通过添加有关实例的额外信息来提供一些额外的功能。
>>> motifsM[0].consensus
Seq('GCTTATGTAA')
>>> motifsM[0].alignment.sequences[0].sequence_name
'iYFL005W'
>>> motifsM[0].alignment.sequences[0].sequence_id
'sequence_15'
>>> motifsM[0].alignment.sequences[0].start
480
>>> motifsM[0].alignment.sequences[0].strand
'+'
>>> motifsM[0].alignment.sequences[0].pvalue
1.97e-06