使用 Bio.motifs 進行序列基序分析
本章概述 Biopython 中包含的 Bio.motifs 套件的功能。它適用於參與序列基序分析的人，因此我假設您熟悉基序分析的基本概念。如果有些地方不清楚，請查看第 有用的連結 節以獲取一些相關連結。
本章大部分內容描述了 Biopython 1.61 及更高版本中包含的新的 Bio.motifs 套件，它正在取代 Biopython 1.50 中引入的舊的 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"),
... ]
然後我們可以建立一個基序物件，如下所示
>>> m = motifs.create(instances)
建立此基序的實例儲存在 .alignment 屬性中
>>> print(m.alignment.sequences)
[Seq('TACAA'), Seq('TACGC'), Seq('TACAC'), Seq('TACCC'), Seq('AACCC'), Seq('AATGC'), Seq('AATGC')]
列印基序物件會顯示建立它的實例
>>> print(m)
TACAA
TACGC
TACAC
TACCC
AACCC
AATGC
AATGC
基序的長度定義為序列長度，這對於所有實例都應該相同
>>> len(m)
5
基序物件有一個屬性 .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 的慣例，並採用以下引數
- substitution_matrix
比較序列時使用的評分矩陣。預設情況下，它是
None，在這種情況下，我們僅計算每個字母的頻率。您可以使用Bio.Align.substitution\_matrices中提供的替換矩陣，而不是預設值。常見的選擇是蛋白質的 BLOSUM62（也稱為 EBLOSUM62）和核苷酸的 NUC.4.4（也稱為 EDNAFULL）。注意：目前，此方法尚未針對預設值None以外的值實作。- plurality
達到一致所需的正向匹配數閾值，除以列中的總計數。如果
substitution_matrix是None，則此引數也必須是None，否則會忽略；否則會引發ValueError。如果substitution_matrix不是None，則 plurality 的預設值為 0.5。- identity
定義一致值所需的相同值數，除以列中的總計數。如果相同值的數目小於 identity 乘以列中的總計數，則在一致序列中使用未定義的字元（核苷酸為
N，氨基酸序列為X）。如果identity為 1.0，則只有相同的字母列才會產生一致。預設值為零。- setcase
正向匹配的閾值，除以列中的總計數，高於該閾值，則一致為大寫，低於該閾值，則一致為小寫。預設情況下，這等於 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}\) – 在第 \(j\) 列中字母 \(i\) 的觀察到的頻率（正規化）；（參見下文）；
\(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 屬性取得 motif m 的每一欄的相對熵
>>> m.relative_entropy
array([1.01477186, 2. , 1.13687943, 0.44334329, 1.40832722])
這些值使用以 2 為底的對數計算，因此單位為位元 (bits)。第二欄（僅由 A 核苷酸組成）具有最高的相對熵；第四欄（由 A、C 或 G 核苷酸組成）具有最低的相對熵。motif 的相對熵可以通過對各欄求和來計算
>>> print(f"Relative entropy is {sum(m.relative_entropy):0.5f}")
Relative entropy is 6.00332
建立序列標誌圖
如果我們可以連上網路，就可以建立一個 weblogo
>>> m.weblogo("mymotif.png")
我們應該將標誌圖以 PNG 格式儲存到指定檔案中。
讀取 motifs
手動從 instances 建立 motifs 有點枯燥，因此有一些用於讀寫 motifs 的 I/O 函式會很有用。目前沒有任何真正完善的 motifs 儲存標準，但有幾種格式比其他格式更常用。
JASPAR
JASPAR 是最流行的 motif 資料庫之一。除了 motif 序列資訊外，JASPAR 資料庫還為每個 motif 儲存了許多元資訊。Bio.motifs 模組包含一個特殊的類別 jaspar.Motif，其中此元資訊表示為屬性
matrix_id- 唯一的 JASPAR motif ID，例如 'MA0004.1'name- TF 的名稱，例如 'Arnt'collection- motif 所屬的 JASPAR 集合，例如 'CORE'tf_class- 此 TF 的結構類別，例如 'Zipper-Type'tf_family- 此 TF 所屬的家族，例如 'Helix-Loop-Helix'species- 此 TF 所屬的物種，可能有多個值，這些值指定為分類 ID，例如 10090tax_group- 此 motif 所屬的分類超群，例如 'vertebrates'acc- TF 蛋白的登錄號，例如 'P53762'data_type- 用於建構此 motif 的資料類型，例如 'SELEX'medline- 支持此 motif 的文獻的 Pubmed ID，可能有多個值，例如 7592839pazar_id- PAZAR 資料庫中 TF 的外部參考，例如 'TF0000003'comment- 包含關於 motif 建構的註釋的自由格式文字
jaspar.Motif 類別繼承自通用的 Motif 類別，因此提供了任何 motif 格式的所有功能 - 讀取 motifs、寫入 motifs、掃描序列以尋找 motif instances 等。
JASPAR 以多種不同的方式儲存 motifs，包括三種不同的純文字檔案格式和一個 SQL 資料庫。所有這些格式都有助於建構計數矩陣。但是，上述可用元資訊的數量會隨格式而異。
JASPAR sites 格式
三種純文字檔案格式中的第一種包含 instances 的列表。例如，以下是 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 instances。
我們可以從這些 instances 建立一個 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
此 motif 的計數矩陣會自動從 instances 計算
>>> 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 也直接以計數矩陣的形式提供 motifs，而無需提供建立它的 instances。此 pfm 格式僅儲存單個 motif 的計數矩陣。例如，以下是包含人類 SRF 轉錄因子計數矩陣的 JASPAR 檔案 SRF.pfm
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
我們可以為此計數矩陣建立一個 motif，如下所示
>>> 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
由於此 motif 是直接從計數矩陣建立的，因此它沒有與之相關聯的 instances
>>> print(srf.alignment)
None
我們現在可以要求這兩個 motifs 的一致序列
>>> print(arnt.counts.consensus)
CACGTG
>>> print(srf.counts.consensus)
GCCCATATATGG
與 instances 檔案一樣，此格式中不儲存任何元資訊。
JASPAR 格式 jaspar
jaspar 檔案格式允許在單個檔案中指定多個 motifs。在此格式中，每個 motif 記錄都包含一個標頭行，後跟四行定義計數矩陣的行。標頭行以 > 字元開頭（類似於 Fasta 檔案格式），後跟唯一的 JASPAR 矩陣 ID 和 TF 名稱。以下範例顯示了包含三個 motifs Arnt、RUNX1 和 MEF2A 的 jaspar 格式檔案
>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 ]
motifs 讀取方式如下
>>> 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 motif 會產生計數資料和可用的元資訊。
存取 JASPAR 資料庫
除了剖析這些純文字檔案格式外，我們還可以從 JASPAR SQL 資料庫檢索 motifs。與純文字檔案格式不同，JASPAR 資料庫允許儲存 JASPAR Motif 類別中定義的所有可能的元資訊。如何設定 JASPAR 資料庫已超出本文檔的範圍（請參閱主要 JASPAR 網站）。使用 Bio.motifs.jaspar.db 模組從 JASPAR 資料庫讀取 motifs。首先使用 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 獲取單個 motif。請注意，JASPAR ID 由一個基本 ID 和一個版本號組成，並以小數點分隔，例如 'MA0004.1'。fetch_motif_by_id 方法允許您使用完全指定的 ID 或僅使用基本 ID。如果僅提供基本 ID，則會傳回 motif 的最新版本。
>>> arnt = jdb.fetch_motif_by_id("MA0004")
列印 motif 會顯示 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
我們也可以通過名稱獲取 motifs。名稱必須完全匹配（目前不支援部分匹配或資料庫萬用字元）。請注意，由於無法保證名稱是唯一的，因此 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 方法允許您獲取符合一組指定條件的 motifs。這些條件包括任何上述的元資訊，以及某些矩陣屬性，例如最小資訊內容 (min_ic 在以下範例中)、矩陣的最小長度或用於建構矩陣的最小站點數。僅傳回通過所有指定條件的 motifs。請注意，與允許多個值的元資訊對應的選擇條件可以指定為單個值或值列表，例如以下範例中的 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 模組相容。因此，關於背景和偽計數的預設選擇以及資訊內容的計算方式和序列中搜尋 instances 的方式的某些細節都基於此相容性標準。這些選擇在下面的特定小節中說明。
- 背景選擇Perl
TFBS模組似乎允許選擇自訂背景機率（儘管文件聲明假設均勻背景）。但是，預設值是使用均勻背景。因此，建議您使用均勻背景來計算位置特異性評分矩陣 (PSSM)。這是使用 Biopythonmotifs模組時的預設值。 - 偽計數選擇依預設，Perl
TFBS模組使用等於 \(\sqrt{N} * \textrm{bg}[\textrm{nucleotide}]\) 的偽計數，其中 \(N\) 表示用於建構矩陣的序列總數。若要應用相同的偽計數公式，請使用jaspar.calculate\_pseudcounts()函式設定 motifpseudocounts屬性>>> motif.pseudocounts = motifs.jaspar.calculate_pseudocounts(motif)
請注意，計數矩陣的列中可能具有不等的序列數。偽計數計算使用組成矩陣的平均序列數。但是，當在計數矩陣上呼叫
normalize時，列中的每個計數值都會除以組成該特定列的序列總數，而不是平均序列數。這與 PerlTFBS模組不同，因為正規化不是作為單獨的步驟完成的，因此在 pssm 的整個計算過程中都使用了平均序列數。因此，對於具有不等列計數的矩陣，motifs模組計算的 PSSM 會與 PerlTFBS模組計算的 pssm 有所不同。 - 矩陣資訊內容的計算矩陣的訊息內容 (IC) 或特異性是使用
PositionSpecificScoringMatrix類別的mean方法計算的。然而，值得注意的是，在 Perl 的TFBS模組中，預設行為是在不先應用偽計數的情況下計算 IC，即使預設情況下，PSSM 是如上所述使用偽計數計算的。 - 搜尋實例使用 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 格式的檔案。例如，列之間的空格數可能不同，或者可能使用 Tab 而不是空格。使用 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 引數。例如，針對 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 物件。
搜尋實例
motif 最常見的用途是在某些序列中尋找其出現的實例。為了本節的目的，我們將使用像這樣的虛擬序列
>>> test_seq = Seq("TACACTGCATTACAACCCAAGCATTA")
>>> len(test_seq)
26
搜尋完全匹配
尋找實例最簡單的方法是尋找 motif 真實實例的完全匹配項
>>> 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 分數搜尋匹配項
尋找導致針對我們 motif 的高對數勝算分數的位置同樣容易
>>> 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
負位置指的是在測試序列的反向鏈上找到的 motif 實例，並遵循 Python 關於負索引的慣例。因此，在 pos 的 motif 實例位於 test_seq[pos:pos+len(m)]，無論 pos 的值是正數還是負數。
您可能會注意到閾值參數，這裡任意設定為 \(3.0\)。這是以 \(log_2\) 為單位，因此我們現在只尋找在 motif 模型下發生的機率比在背景中高出八倍的詞語。預設閾值為 \(0.0\)，它會選擇所有看起來比背景更像 motif 的詞語。
您也可以計算序列中所有位置的分數
>>> 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 分數的分佈。由於分數分佈的空間隨著 motif 長度呈指數增長，我們使用給定精度的近似值來保持計算成本可控
>>> distribution = pssm.distribution(background=background, precision=10**4)
distribution 物件可用於確定許多不同的閾值。我們可以指定要求的偽陽性率（在背景生成的序列中「找到」motif 實例的機率）
>>> threshold = distribution.threshold_fpr(0.01)
>>> print("%5.3f" % threshold)
4.009
或偽陰性率（「未找到」從 motif 生成的實例的機率）
>>> 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
例如，就我們的 motif 而言，您可以取得閾值，讓您獲得與搜尋具有 \(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
每個 motif 物件都有一個關聯的位置特異性評分矩陣
為了方便使用 PSSM 搜尋潛在的 TFBS，位置權重矩陣和位置特異性評分矩陣都與每個 motif 關聯。以 Arnt motif 為例
>>> 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
比較 motif
一旦我們有多個 motif，我們可能想比較它們。
在我們開始比較 motif 之前，我應該指出 motif 邊界通常相當隨意。這表示我們經常需要比較不同長度的 motif，因此比較需要涉及某種對齊方式。這表示我們必須考慮兩件事
motif 的對齊
比較對齊的 motif 的某些函數
為了對齊 motif，我們使用 PSSM 的無間隙對齊，並將零取代矩陣開頭和結尾的任何遺失列。這表示我們實際上將背景分佈用於 PSSM 中遺失的列。距離函數接著會傳回 motif 之間的最小距離，以及對齊中的對應偏移。
為了舉例說明，讓我們先載入另一個與我們的測試 motif m 相似的 motif
>>> 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
為了使 motif 具有可比性，我們為偽計數和背景分佈選擇與我們的 motif 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
我們將使用皮爾森相關性來比較這些 motif。由於我們希望它類似於距離度量，我們實際上會取 \(1-r\)，其中 \(r\) 是皮爾森相關係數 (PCC)
>>> distance, offset = pssm.dist_pearson(pssm_reb1)
>>> print("distance = %5.3g" % distance)
distance = 0.239
>>> print(offset)
-2
這表示 motif m 和 m_reb1 之間最佳的 PCC 是透過下列對齊方式取得
m: bbTACGCbb
m_reb1: GTTACCCGG
其中 b 代表背景分佈。PCC 本身大約是 \(1-0.239=0.761\)。
從頭 motif 尋找
目前，Biopython 對於從頭 motif 尋找的支援有限。也就是說，我們支援執行 xxmotif 以及剖析 MEME。由於 motif 尋找工具的數量快速增長，因此歡迎提供新的剖析器。
MEME
假設您已在您選擇的序列上使用您喜歡的參數執行 MEME，並將輸出儲存在檔案 meme.out 中。您可以執行下列程式碼來檢索 MEME 報告的 motif
>>> 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>]
除了最想要的 motif 清單之外，結果物件還包含更多實用資訊，可透過具有不言自明的名稱的屬性存取
.alphabet.datafile.sequences.version.command
MEME 剖析器傳回的 motif 可以像一般的 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