使用生物基序进行序列基序分析
本章概述了 Bio.motifs Biopython中包含的包。它适合参与序列基序分析的人，所以我假设您熟悉基序分析的基本概念。如果有什么不清楚，请查看部分 有用的链接 获取一些相关链接。
本章的大部分内容都描述了新的 Bio.motifs Biopython 1.61以后包含的包，该包正在取代旧的 Bio.Motif Biopython 1.50中引入的包，该包又基于两个较旧的前Biopython模块， Bio.AlignAce 和 Bio.MEME .它通过统一的主题对象实现来提供大部分功能。
说到其他图书馆，如果您正在阅读这篇文章，您可能会感兴趣 TAMO ，另一个旨在处理序列图案的Python库。它支持更多 de-novo 主题查找器，但它不是Biopython的一部分，并且对商业用途有一些限制。
Motif对象
由于我们对主题分析感兴趣，所以我们需要看看 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]
但您也可以将其视为一个2D阵列，其中核苷为第一维，位置为第二维：
>>> 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 方法可以让您详细指定如何计算共识序列。这种方法在很大程度上遵循CLASS程序的惯例 cons ，并采取以下论点：
- substitution_matrix
比较序列时使用的评分矩阵。默认为
None，在这种情况下，我们只需计算每个字母的频率。您可以使用中提供的替代矩阵，而不是默认值Bio.Align.substitution\_matrices.常见的选择是蛋白质的BLOSUM62（也称为EBLOSUM62），以及核苷酸的NUC.4.4（也称为EDNAAFSYS）。注：目前，该方法尚未针对默认值以外的值实现None.- 多
达到共识所需的阳性匹配数量的阈值，除以列中的总计数。如果
substitution_matrix是None，那么这个论点也一定是None，并被忽略; aValueError否则提出。如果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')
相对熵
相对熵（或库尔巴克-莱布勒距离） \(H_j\) 柱 \(j\) 主题的定义为 [Schneider1986] [Durbin1998] :
where:
\(M\) - 字母表中的字母数（由
len(m.alphabet));\(p_{ij}\) - 观察到的字母频率 \(i\) ，正常化，在 \(j\) - 第th列（见下文）;
\(b_{i}\) - 字母的背景概率 \(i\) （由
m.background[i]).
观察到的频率 \(p_{ij}\) 计算方法如下：
where:
\(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：
如果背景分布均匀，则相对信息量相同。
每列主题的相对熵 m 可以使用 relative_entropy 属性：
>>> m.relative_entropy
array([1.01477186, 2. , 1.13687943, 0.44334329, 1.40832722])
这些值是使用以2为底的log计算的，因此以位为单位。第二列（由 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 .除了基序序列信息外，JASVAR数据库还存储了每个基序的大量元信息。模块 Bio.motifs 包含一个专门的类 jaspar.Motif 其中该元信息被表示为属性：
matrix_id- 唯一的JASVAR基序ID，例如“MA 0004.1”name- TF的名称，例如“Arnt”collection- 该主题所属的JASVAR系列，例如“CORE”tf_class- 该TF的结构类别，例如“拉链型”tf_family- 该TF所属的家族，例如“Helix-Loop-Hattan”species- 此TF所属的物种可能具有多个值，这些值指定为分类学ID，例如10090tax_group- 该主题所属的分类超群，例如“脊椎动物”acc- TF蛋白的登录号，例如“P53762”data_type- 用于构建此基序的数据类型，例如“SEN”medline- 支持此主题的文献的Pubmed ID可能是多个值，例如7592839pazar_id- PAZAR数据库中TF的外部引用，例如“TF 0000003”comment- 包含有关主题结构注释的自由形式文本
的 jaspar.Motif 类继承自类 Motif 类，因此提供了任何主题格式的所有设施-读取主题、写入主题、扫描主题实例的序列等。
JASBAR以多种不同的方式存储主题，包括三种不同的平面文件格式和SQL数据库。所有这些格式都有助于构建计数矩阵。然而，上述可用的Meta信息的量随格式而变化。
JASVAR sites 格式
三种平面文件格式中的第一种包含实例列表。例如，这些是JASVAR的开始行和结束行 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
此格式不存储任何Meta信息。
JASVAR pfm 格式
JASVAR还使图案直接作为计数矩阵提供，而不需要创建它的实例。这 pfm 格式仅存储单个主题的计数矩阵。例如，这是JASVAR文件 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
与实例文件一样，没有以这种格式存储Meta信息。
JASPAR格式 jaspar
的 jaspar 文件格式允许在单个文件中指定多个图案。在这种格式中，每条主题记录都由标题行和定义计数矩阵的四行组成。标题行以一个开头 > 字符（类似于Fasta文件格式），后面是唯一的JASVAR矩阵ID和TF名称。下面的示例演示一个 jaspar 包含三个主题Arnt、RUNX 1和MEF 2A的格式化文件：
>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
请注意，打印JASVAR主题会产生计数数据和可用的元信息。
初始化JASVAR数据库
除了解析这些平面文件格式之外，我们还可以从JASPARSQL数据库中检索主题。与平面文件格式不同，JASPAR数据库允许存储JASPAR中定义的所有可能的Meta信息 Motif class. It is beyond the scope of this document to describe how to set up a JASPAR database (please see the main JASPAR 网站）。使用 Bio.motifs.jaspar.db module.首先使用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,
... )
现在，我们可以通过其唯一的JASPAR ID获取单个主题， fetch_motif_by_id 法请注意，JASPARID由基本ID和版本号组成，版本号由小数点分隔，例如“MA 0004.1”。的 fetch_motif_by_id 方法允许您使用完全指定的ID或仅使用基本ID。如果仅提供基本ID，则返回主题的最新版本。
>>> arnt = jdb.fetch_motif_by_id("MA0004")
打印主题表明JASPARSQL数据库存储的元信息比平面文件多得多：
>>> 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 方法允许您获取与指定条件集相匹配的图案。这些标准包括任何上述Meta信息以及某些矩阵属性，例如最小信息内容 (min_ic 在下面的例子中）、矩阵的最小长度或用于构建矩阵的最小位点数量。仅返回通过所有指定条件的主题。请注意，与允许多个值的Meta信息相对应的选择标准可以指定为单个值或值列表，例如 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模块兼容
需要注意的一件重要事情是JASVAR Motif class was designed to be compatible with the popular Perl TFBS modules .因此，有关背景和伪计数默认值选择以及如何计算信息内容和搜索实例序列的一些细节都基于此兼容性标准。这些选择在下面的具体小节中注明。
- Choice of background:的Perl
TFBS模块似乎允许选择自定义背景概率（尽管文档声明假设背景均匀）。但默认设置是使用统一的背景。因此，建议您使用统一的背景来计算位置特定的评分矩阵（PSSM）。这是使用Biopython时的默认设置motifsmodule. - Choice of pseudocounts:By default, the Perl
TFBSmodules use a pseudocount equal to \(\sqrt{N} * \textrm{bg}[\textrm{nucleotide}]\), where \(N\) represents the total number of sequences used to construct the matrix. To apply this same pseudocount formula, set the motifpseudocountsattribute using thejaspar.calculate\_pseudcounts()function:>>> motif.pseudocounts = motifs.jaspar.calculate_pseudocounts(motif)
请注意，计数矩阵的组成列的序列数量可能不相等。伪计数计算使用组成矩阵的序列的平均数。然而当
normalize在计数矩阵上调用时，一列中的每个计数值都除以组成该特定列的序列总数，而不是除以序列的平均数。这与Perl不同TFBS模块，因为规范化不是作为单独的步骤进行的，因此在整个PSM的计算中使用平均序列数。因此，对于列数不等的矩阵，由motifs模块将与Perl计算的psms有所不同TFBS模块。 - Computation of matrix information content:矩阵的信息含量（IC）或特异性是使用
mean方法PositionSpecificScoringMatrix课然而值得注意的是，在Perl中TFBS模块的默认行为是在不首先应用伪计数的情况下计算IC，尽管默认情况下PSSM是使用如上所述的伪计数来计算的。 - Searching for instances:使用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或蛋白质序列作为输入，并根据要求输出尽可能多的基序。因此，与JASVAR文件相比，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 课类从列表继承，您可以想到 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]
Seq('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
TRANSAC是一个手动策划的转录因子及其基因组结合位点和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
//
该文件显示了主题的频率矩阵 motif1 12个核苷酸。一般来说，一个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
//
要解析TRANSAC文件，请使用
>>> with open("transfac.dat") as handle:
... record = motifs.parse(handle, "TRANSFAC")
...
如果检测到文件内容与TRANSFAC文件格式之间存在任何差异， ValueError 被提出。请注意，您可能会遇到不严格遵循TRANSAC格式的文件。例如，列之间的空白数量可能不同，或者可以使用选项卡来代替空白。使用 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文件中常见的字段 列出了TRANSAC文件中常见的两字母字段代码：
|
登录号 |
|
入学号，中学 |
|
统计基础 |
|
结合因子 |
|
矩阵下的因子结合位点 |
|
评论 |
|
版权声明 |
|
简短因素描述 |
|
外部数据库 |
|
创建/更新日期 |
|
亚家族 |
|
超家族 |
|
标识符 |
|
结合因子名称 |
|
系统分类 |
|
物种/分类群 |
|
旧版本 |
|
优选版本 |
|
类型 |
|
空行;这些不存储在记录中。 |
每个主题也有一个属性 .references 包含与主题相关的引用，使用以下两个字母的键：
|
附图标记 |
|
参考作者 |
|
参考数据 |
|
参考标题 |
|
PubMed ID |
打印主题会以其原生的TRANSAC格式写出：
>>> 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)
...
通用 pfm-four-columns 格式
如果任何工具特定的主题格式都不适用于您的PLM文件，并且您的PLM文件的值以4列格式组织，则您可以尝试通用 pfm-four-columns 主题解析器：
# CIS-BP
Pos A C G T
1 0.00961538461538462 0.00961538461538462 0.00961538461538462 0.971153846153846
2 0.00961538461538462 0.00961538461538462 0.00961538461538462 0.971153846153846
3 0.971153846153846 0.00961538461538462 0.00961538461538462 0.00961538461538462
4 0.00961538461538462 0.00961538461538462 0.00961538461538462 0.971153846153846
5 0.00961538461538462 0.971153846153846 0.00961538461538462 0.00961538461538462
6 0.971153846153846 0.00961538461538462 0.00961538461538462 0.00961538461538462
7 0.00961538461538462 0.971153846153846 0.00961538461538462 0.00961538461538462
8 0.00961538461538462 0.00961538461538462 0.00961538461538462 0.971153846153846
# C2H2-ZFs
Gene ENSG00000197372
Pos A C G T
1 0.341303 0.132427 0.117054 0.409215
2 0.283785 0.077066 0.364552 0.274597
3 0.491055 0.078208 0.310520 0.120217
4 0.492621 0.076117 0.131007 0.300256
5 0.250645 0.361464 0.176504 0.211387
6 0.276694 0.498070 0.197793 0.027444
7 0.056317 0.014631 0.926202 0.002850
8 0.004470 0.007769 0.983797 0.003964
9 0.936213 0.058787 0.002387 0.002613
10 0.004352 0.004030 0.002418 0.989200
11 0.013277 0.008165 0.001991 0.976567
12 0.968132 0.002263 0.002868 0.026737
13 0.397623 0.052017 0.350783 0.199577
14 0.000000 0.000000 1.000000 0.000000
15 1.000000 0.000000 0.000000 0.000000
16 0.000000 0.000000 1.000000 0.000000
17 0.000000 0.000000 1.000000 0.000000
18 1.000000 0.000000 0.000000 0.000000
19 0.000000 1.000000 0.000000 0.000000
20 1.000000 0.000000 0.000000 0.000000
# C2H2-ZFs
Gene FBgn0000210
Motif M1734_0.90
Pos A C G T
1 0.25 0.0833333 0.0833333 0.583333
2 0.75 0.166667 0.0833333 0
3 0.833333 0 0 0.166667
4 1 0 0 0
5 0 0.833333 0.0833333 0.0833333
6 0.333333 0 0 0.666667
7 0.833333 0 0 0.166667
8 0.5 0 0.333333 0.166667
9 0.5 0.0833333 0.166667 0.25
10 0.333333 0.25 0.166667 0.25
11 0.166667 0.25 0.416667 0.166667
# FlyFactorSurvey (Cluster Buster)
>AbdA_Cell_FBgn0000014
1 3 0 14
0 0 0 18
16 0 0 2
18 0 0 0
1 0 0 17
0 0 6 12
15 1 2 0
# HOMER
>ATGACTCATC AP-1(bZIP)/ThioMac-PU.1-ChIP-Seq(GSE21512)/Homer 6.049537 -1.782996e+03 0 9805.3,5781.0,3085.1,2715.0,0.00e+00
0.419 0.275 0.277 0.028
0.001 0.001 0.001 0.997
0.010 0.002 0.965 0.023
0.984 0.003 0.001 0.012
0.062 0.579 0.305 0.054
0.026 0.001 0.001 0.972
0.043 0.943 0.001 0.012
0.980 0.005 0.001 0.014
0.050 0.172 0.307 0.471
0.149 0.444 0.211 0.195
# HOCOMOCO
> AHR_si
40.51343240527031 18.259112547756697 56.41253757072521 38.77363485291994
10.877470982533044 11.870876719950774 34.66312982331297 96.54723985087516
21.7165707818416 43.883079837598544 20.706746561638717 67.6523201955933
2.5465132509466635 1.3171620263517245 145.8637051322628 4.231336967110781
0.0 150.35847450464382 1.4927836298652875 2.1074592421627525
3.441039751299748 0.7902972158110341 149.37613720253387 0.3512432070271259
0.0 3.441039751299748 0.7024864140542533 149.81519121131782
0.0 0.0 153.95871737667187 0.0
43.07922333291745 66.87558226865211 16.159862546986584 27.844049228115868
# Neph
UW.Motif.0001 atgactca
0.772949 0.089579 0.098612 0.038860
0.026652 0.004653 0.025056 0.943639
0.017663 0.023344 0.918728 0.040264
0.919596 0.025414 0.029759 0.025231
0.060312 0.772259 0.104968 0.062462
0.037406 0.020643 0.006667 0.935284
0.047316 0.899024 0.026928 0.026732
0.948639 0.019497 0.005737 0.026128
# Tiffin
T A G C
30 0 28 40
0 0 0 99
0 55 14 29
0 99 0 0
20 78 0 0
0 52 7 39
19 46 11 22
0 60 38 0
0 33 0 66
73 0 25 0
99 0 0 0
主题如下：
>>> with open("fourcolumns.pfm") as fh:
... for m in motifs.parse(fh, "pfm-four-columns"):
... print(m.name, m.counts, sep="\n")
...
0 1 2 3 4 5 6 7
G: 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
A: 0.01 0.01 0.97 0.01 0.01 0.97 0.01 0.01
T: 0.97 0.97 0.01 0.97 0.01 0.01 0.01 0.97
C: 0.01 0.01 0.01 0.01 0.97 0.01 0.97 0.01
ENSG00000197372
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
G: 0.12 0.36 0.31 0.13 0.18 0.20 0.93 0.98 0.00 0.00 0.00 0.00 0.35 1.00 0.00 1.00 1.00 0.00 0.00 0.00
A: 0.34 0.28 0.49 0.49 0.25 0.28 0.06 0.00 0.94 0.00 0.01 0.97 0.40 0.00 1.00 0.00 0.00 1.00 0.00 1.00
T: 0.41 0.27 0.12 0.30 0.21 0.03 0.00 0.00 0.00 0.99 0.98 0.03 0.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00
C: 0.13 0.08 0.08 0.08 0.36 0.50 0.01 0.01 0.06 0.00 0.01 0.00 0.05 0.00 0.00 0.00 0.00 0.00 1.00 0.00
M1734_0.90
0 1 2 3 4 5 6 7 8 9 10
G: 0.08 0.08 0.00 0.00 0.08 0.00 0.00 0.33 0.17 0.17 0.42
A: 0.25 0.75 0.83 1.00 0.00 0.33 0.83 0.50 0.50 0.33 0.17
T: 0.58 0.00 0.17 0.00 0.08 0.67 0.17 0.17 0.25 0.25 0.17
C: 0.08 0.17 0.00 0.00 0.83 0.00 0.00 0.00 0.08 0.25 0.25
AbdA_Cell_FBgn0000014
0 1 2 3 4 5 6
G: 0.00 0.00 0.00 0.00 0.00 6.00 2.00
A: 1.00 0.00 16.00 18.00 1.00 0.00 15.00
T: 14.00 18.00 2.00 0.00 17.00 12.00 0.00
C: 3.00 0.00 0.00 0.00 0.00 0.00 1.00
ATGACTCATC AP-1(bZIP)/ThioMac-PU.1-ChIP-Seq(GSE21512)/Homer 6.049537 -1.782996e+03 0 9805.3,5781.0,3085.1,2715.0,0.00e+00
0 1 2 3 4 5 6 7 8 9
G: 0.28 0.00 0.96 0.00 0.30 0.00 0.00 0.00 0.31 0.21
A: 0.42 0.00 0.01 0.98 0.06 0.03 0.04 0.98 0.05 0.15
T: 0.03 1.00 0.02 0.01 0.05 0.97 0.01 0.01 0.47 0.20
C: 0.28 0.00 0.00 0.00 0.58 0.00 0.94 0.01 0.17 0.44
AHR_si
0 1 2 3 4 5 6 7 8
G: 56.41 34.66 20.71 145.86 1.49 149.38 0.70 153.96 16.16
A: 40.51 10.88 21.72 2.55 0.00 3.44 0.00 0.00 43.08
T: 38.77 96.55 67.65 4.23 2.11 0.35 149.82 0.00 27.84
C: 18.26 11.87 43.88 1.32 150.36 0.79 3.44 0.00 66.88
0 1 2 3 4 5 6 7
G: 0.10 0.03 0.92 0.03 0.10 0.01 0.03 0.01
A: 0.77 0.03 0.02 0.92 0.06 0.04 0.05 0.95
T: 0.04 0.94 0.04 0.03 0.06 0.94 0.03 0.03
C: 0.09 0.00 0.02 0.03 0.77 0.02 0.90 0.02
0 1 2 3 4 5 6 7 8 9 10
G: 28.00 0.00 14.00 0.00 0.00 7.00 11.00 38.00 0.00 25.00 0.00
A: 0.00 0.00 55.00 99.00 78.00 52.00 46.00 60.00 33.00 0.00 0.00
T: 30.00 0.00 0.00 0.00 20.00 0.00 19.00 0.00 0.00 73.00 99.00
C: 40.00 99.00 29.00 0.00 0.00 39.00 22.00 0.00 66.00 0.00 0.00
通用 pfm-four-rows 格式
如果任何工具特定的主题格式都不适用于您的PLM文件，并且您的PLM文件的值以4行格式组织，则可以尝试通用 pfm-four-rows 主题解析器：
# hDPI
A 0 5 6 5 1 0
C 1 1 0 0 0 4
G 5 0 0 0 3 0
T 0 0 0 1 2 2
# YeTFaSCo
A 0.5 0.0 0.0 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.5 0.0 0.0833333334583333
T 0.0 0.0 0.0 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.0 0.0 0.0833333334583333
G 0.0 1.0 0.0 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.0 1.0 0.249999999875
C 0.5 0.0 1.0 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.5 0.0 0.583333333208333
# FlyFactorSurvey ZFP finger
A | 92 106 231 135 0 1 780 28 0 700 739 94 60 127 130
C | 138 82 129 81 774 1 3 1 0 6 17 49 193 122 148
G | 270 398 54 164 7 659 1 750 755 65 1 41 202 234 205
T | 290 204 375 411 9 127 6 11 36 20 31 605 335 307 308
# ScerTF pcm
A | 9 1 1 97 1 94
T | 80 1 97 1 1 2
C | 9 97 1 1 1 2
G | 2 1 1 1 97 2
# ScerTF pfm
A | 0.090 0.010 0.010 0.970 0.010 0.940
C | 0.090 0.970 0.010 0.010 0.010 0.020
G | 0.020 0.010 0.010 0.010 0.970 0.020
T | 0.800 0.010 0.970 0.010 0.010 0.020
# iDMMPMM
> abd-A
0.218451749734889 0.0230646871686108 0.656680805938494 0.898197242841994 0.040694591728526 0.132953340402969 0.74907211028632 0.628313891834571
0.0896076352067868 0.317338282078473 0.321580063626723 0.0461293743372216 0.0502386002120891 0.040694591728526 0.0284994697773065 0.0339342523860021
0.455991516436904 0.0691940615058324 0.0108695652173913 0.0217391304347826 0.0284994697773065 0.0284994697773065 0.016304347826087 0.160127253446448
0.235949098621421 0.590402969247084 0.0108695652173913 0.0339342523860021 0.880567338282079 0.797852598091198 0.206124072110286 0.17762460233298
# JASPAR
>MA0001.1 AGL3
A [ 0 3 79 40 66 48 65 11 65 0 ]
C [94 75 4 3 1 2 5 2 3 3 ]
G [ 1 0 3 4 1 0 5 3 28 88 ]
T [ 2 19 11 50 29 47 22 81 1 6 ]
# JASPAR
>MA0001.1 AGL3
0 3 79 40 66 48 65 11 65 0
94 75 4 3 1 2 5 2 3 3
1 0 3 4 1 0 5 3 28 88
2 19 11 50 29 47 22 81 1 6
# Cys2His2 Zinc Finger Proteins PWM Predictor
base 1 2 3 4 5 6 7 8 9
a 0.116 0.974 0.444 0.116 0.974 0.444 0.667 0.939 0.068 # noqa: RST301
c 0.718 0.006 0.214 0.718 0.006 0.214 0.143 0.006 0.107 # noqa: RST301
g 0.016 0.020 0.028 0.016 0.020 0.028 0.047 0.045 0.216 # noqa: RST301
t 0.150 0.001 0.314 0.150 0.001 0.314 0.143 0.009 0.609 # noqa: RST301
Ent= 1.210 0.202 1.665 1.210 0.202 1.665 1.399 0.396 1.521
主题如下：
>>> with open("fourrows.pfm") as fh:
... for m in motifs.parse(fh, "pfm-four-rows"):
... print(m.name, m.counts, sep="\n")
...
0 1 2 3 4 5
G: 5.00 0.00 0.00 0.00 3.00 0.00
A: 0.00 5.00 6.00 5.00 1.00 0.00
T: 0.00 0.00 0.00 1.00 2.00 2.00
C: 1.00 1.00 0.00 0.00 0.00 4.00
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
G: 0.00 1.00 0.00 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.00 1.00 0.25
A: 0.50 0.00 0.00 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.50 0.00 0.08
T: 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.00 0.00 0.08
C: 0.50 0.00 1.00 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.50 0.00 0.58
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
G: 270.00 398.00 54.00 164.00 7.00 659.00 1.00 750.00 755.00 65.00 1.00 41.00 202.00 234.00 205.00
A: 92.00 106.00 231.00 135.00 0.00 1.00 780.00 28.00 0.00 700.00 739.00 94.00 60.00 127.00 130.00
T: 290.00 204.00 375.00 411.00 9.00 127.00 6.00 11.00 36.00 20.00 31.00 605.00 335.00 307.00 308.00
C: 138.00 82.00 129.00 81.00 774.00 1.00 3.00 1.00 0.00 6.00 17.00 49.00 193.00 122.00 148.00
0 1 2 3 4 5
G: 2.00 1.00 1.00 1.00 97.00 2.00
A: 9.00 1.00 1.00 97.00 1.00 94.00
T: 80.00 1.00 97.00 1.00 1.00 2.00
C: 9.00 97.00 1.00 1.00 1.00 2.00
0 1 2 3 4 5
G: 0.02 0.01 0.01 0.01 0.97 0.02
A: 0.09 0.01 0.01 0.97 0.01 0.94
T: 0.80 0.01 0.97 0.01 0.01 0.02
C: 0.09 0.97 0.01 0.01 0.01 0.02
abd-A
0 1 2 3 4 5 6 7
G: 0.46 0.07 0.01 0.02 0.03 0.03 0.02 0.16
A: 0.22 0.02 0.66 0.90 0.04 0.13 0.75 0.63
T: 0.24 0.59 0.01 0.03 0.88 0.80 0.21 0.18
C: 0.09 0.32 0.32 0.05 0.05 0.04 0.03 0.03
MA0001.1 AGL3
0 1 2 3 4 5 6 7 8 9
G: 1.00 0.00 3.00 4.00 1.00 0.00 5.00 3.00 28.00 88.00
A: 0.00 3.00 79.00 40.00 66.00 48.00 65.00 11.00 65.00 0.00
T: 2.00 19.00 11.00 50.00 29.00 47.00 22.00 81.00 1.00 6.00
C: 94.00 75.00 4.00 3.00 1.00 2.00 5.00 2.00 3.00 3.00
MA0001.1 AGL3
0 1 2 3 4 5 6 7 8 9
G: 1.00 0.00 3.00 4.00 1.00 0.00 5.00 3.00 28.00 88.00
A: 0.00 3.00 79.00 40.00 66.00 48.00 65.00 11.00 65.00 0.00
T: 2.00 19.00 11.00 50.00 29.00 47.00 22.00 81.00 1.00 6.00
C: 94.00 75.00 4.00 3.00 1.00 2.00 5.00 2.00 3.00 3.00
0 1 2 3 4 5 6 7 8
G: 0.02 0.02 0.03 0.02 0.02 0.03 0.05 0.04 0.22
A: 0.12 0.97 0.44 0.12 0.97 0.44 0.67 0.94 0.07
T: 0.15 0.00 0.31 0.15 0.00 0.31 0.14 0.01 0.61
C: 0.72 0.01 0.21 0.72 0.01 0.21 0.14 0.01 0.11
写作主题
说到出口，让我们总体来看一下出口功能。我们可以使用 format 内置函数在简单的JASVAR中编写主题 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 在JASVAR中写下主题 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]
位置权矩阵
的 .counts Motif对象的属性显示每个核苷酸在排列过程中的每个位置出现的频率。我们可以通过除以比对中的实例数量来对这个矩阵进行标准化，从而得出比对中每个位置上每个核苷酸的概率。我们将这些概率称为位置权重矩阵。然而，请注意，在文献中，该术语也可能用于指代特定于位置的评分矩阵，我们将在下面讨论该矩阵。
通常，在标准化之前将伪计数添加到每个位置。这避免了位置权重矩阵过度适合对齐中有限数量的基序实例，并且还可以防止概率变为零。要将固定伪计数添加到所有位置的所有核苷，请指定 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 论点例如，在GC含量为40%的背景下，使用
>>> 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
PSSM评分相对于特定背景的平均值和标准差通过 .mean 和 .std 方法.
>>> 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分数搜索匹配
寻找位置也很容易，从而产生与我们主题相反的高log赔率分数：
>>> 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
or a threshold (approximately) satisfying some relation between the false-positive rate and the false-negative rate (\(\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 ，分别。如果速度是一个问题并且您想要重复使用脉宽调制或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
这意味着主题之间的最佳PCC m 和 m_reb1 通过以下对齐获得：
m: bbTACGCbb
m_reb1: GTTACCCGG
哪里 b 代表背景分布。PCC本身大致 \(1-0.239=0.761\) .
De novo 主题发现
目前，Biopython仅对 de novo 模体发现也就是说，我们支持跑步 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解析器返回的图案可以完全像常规Motif对象（带有实例）一样对待，它们还通过添加有关实例的额外信息来提供一些额外功能。
>>> 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
有用的链接
Sequence motif 在维基百科上
PWM 在维基百科上
Consensus sequence 在维基百科上