序列註解物件
第 序列物件 章介紹了序列類別。在 Seq
類別「之上」的是序列記錄或 SeqRecord
類別,定義在 Bio.SeqRecord
模組中。此類別允許將更高階的功能(例如識別碼和特徵(作為 SeqFeature
物件))與序列關聯,並在整個序列輸入/輸出介面 Bio.SeqIO
中使用,完整說明請見第 序列輸入/輸出 章。
如果你只會處理像是 FASTA 檔案的簡單資料,你現在可能可以跳過本章。另一方面,如果你將要使用來自 GenBank 或 EMBL 檔案等具有豐富註解的序列資料,此資訊非常重要。
雖然本章應涵蓋有關 SeqRecord
和 SeqFeature
物件的大部分內容,你可能還想閱讀 SeqRecord
wiki 頁面 (https://biopython.dev.org.tw/wiki/SeqRecord) 以及內建文件 (Bio.SeqRecord
和 Bio.SeqFeature
)
>>> from Bio.SeqRecord import SeqRecord
>>> help(SeqRecord)
SeqRecord 物件
SeqRecord
(序列記錄) 類別定義在 Bio.SeqRecord
模組中。此類別允許將更高階的功能(例如識別碼和特徵)與序列關聯(請參閱第 序列物件 章),並且是 Bio.SeqIO
序列輸入/輸出介面的基本資料類型(請參閱第 序列輸入/輸出 章)。
SeqRecord
類別本身非常簡單,並將以下資訊作為屬性提供
- .seq
序列本身,通常是
Seq
物件。- .id
用於識別序列的主要 ID - 一個字串。在大多數情況下,這類似於登錄編號。
- .name
序列的「通用」名稱/ID - 一個字串。在某些情況下,這會與登錄編號相同,但也可能是一個複製名稱。我認為這類似於 GenBank 記錄中的 LOCUS ID。
- .description
序列的人類可讀描述或表達性名稱 - 一個字串。
- .letter_annotations
使用(受限制的)字典保存每個字母的註解,其中包含關於序列中字母的其他資訊。鍵是資訊的名稱,資訊包含在值中,作為與序列本身長度相同的 Python 序列(即列表、元組或字串)。這通常用於品質分數(例如,第 FASTQ 檔案的簡單品質篩選 節)或二級結構資訊(例如,來自斯德哥爾摩/PFAM 比對檔案)。
- .annotations
關於序列的其他資訊字典。鍵是資訊的名稱,資訊包含在值中。這允許將更多「非結構化」資訊添加到序列中。
- .features
具有關於序列特徵的更多結構化資訊的
SeqFeature
物件列表(例如,基因組上基因的位置或蛋白質序列上的網域)。序列特徵的結構將在下方第 特徵、位置與方位物件 節中說明。- .dbxrefs
資料庫交叉引用列表,以字串形式呈現。
建立 SeqRecord
使用 SeqRecord
物件並不複雜,因為所有資訊都以類別的屬性呈現。通常你不會「手動」建立 SeqRecord
,而是使用 Bio.SeqIO
為你讀入序列檔案(請參閱第 序列輸入/輸出 章以及下面的範例)。但是,建立 SeqRecord
可以非常簡單。
從頭開始建立 SeqRecord 物件
要建立 SeqRecord
,你至少需要一個 Seq
物件
>>> from Bio.Seq import Seq
>>> simple_seq = Seq("GATC")
>>> from Bio.SeqRecord import SeqRecord
>>> simple_seq_r = SeqRecord(simple_seq)
此外,你也可以將 id、name 和 description 傳遞到初始化函式,但如果沒有,它們將會設定為指示它們為未知的字串,並且可以隨後修改
>>> simple_seq_r.id
'<unknown id>'
>>> simple_seq_r.id = "AC12345"
>>> simple_seq_r.description = "Made up sequence I wish I could write a paper about"
>>> print(simple_seq_r.description)
Made up sequence I wish I could write a paper about
>>> simple_seq_r.seq
Seq('GATC')
如果你想要將 SeqRecord
輸出到檔案,包含識別碼非常重要。你通常會在建立物件時包含此識別碼
>>> from Bio.Seq import Seq
>>> simple_seq = Seq("GATC")
>>> from Bio.SeqRecord import SeqRecord
>>> simple_seq_r = SeqRecord(simple_seq, id="AC12345")
如上所述,SeqRecord
有一個字典屬性 annotations
。這用於任何不屬於其他更特定屬性的雜項註解。添加註解很容易,只需直接處理註解字典即可
>>> simple_seq_r.annotations["evidence"] = "None. I just made it up."
>>> print(simple_seq_r.annotations)
{'evidence': 'None. I just made it up.'}
>>> print(simple_seq_r.annotations["evidence"])
None. I just made it up.
使用每個字母的註解類似,letter_annotations
是一個類似字典的屬性,可讓你指定任何長度與序列相同的 Python 序列(即字串、列表或元組)
>>> simple_seq_r.letter_annotations["phred_quality"] = [40, 40, 38, 30]
>>> print(simple_seq_r.letter_annotations)
{'phred_quality': [40, 40, 38, 30]}
>>> print(simple_seq_r.letter_annotations["phred_quality"])
[40, 40, 38, 30]
dbxrefs
和 features
屬性只是 Python 列表,應分別用於儲存字串和 SeqFeature
物件(本章稍後討論)。
從 FASTA 檔案建立 SeqRecord 物件
此範例使用一個相當大的 FASTA 檔案,其中包含來自 NCBI 最初下載的鼠疫耶爾辛菌生物變種小鼠 str. 91001 質體 pPCP1 的完整序列。此檔案包含在 Biopython 單元測試的 GenBank 資料夾下,或在我們的網站上線上 NC_005816.fna 。
檔案開頭如下所示 - 你可以檢查是否只存在一筆記錄(即只有一行以大於符號開頭)
>gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus ... pPCP1, complete sequence
TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGGGGGTAATCTGCTCTCC
...
回到第 快速開始 – 你可以用 Biopython 做什麼? 章,你會看到函式 Bio.SeqIO.parse(...)
用於將檔案中的所有記錄作為 SeqRecord
物件迴圈處理。 Bio.SeqIO
模組有一個姊妹函式用於僅包含一筆記錄的檔案,我們將在此處使用(詳細資訊請參閱第 序列輸入/輸出 章)
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.fna", "fasta")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG'), id='gi|45478711|ref|NC_005816.1|', name='gi|45478711|ref|NC_005816.1|', description='gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])
現在,讓我們單獨看看這個 SeqRecord
的主要屬性 - 從 seq
屬性開始,它會給你一個 Seq
物件
>>> record.seq
Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG')
接下來,是識別碼和描述
>>> record.id
'gi|45478711|ref|NC_005816.1|'
>>> record.name
'gi|45478711|ref|NC_005816.1|'
>>> record.description
'gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence'
如您在上方所見,FASTA 記錄的標題行中第一個字(移除大於符號後)會被用作 id
和 name
屬性。整個標題行(移除大於符號後)則會用作記錄描述。這是經過考量的設計,一部分是為了向後相容性,但也適用於像這樣的 FASTA 檔案:
>Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1
TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGGGGGTAATCTGCTCTCC
...
請注意,讀取 FASTA 檔案時,其他註解屬性都不會被填入。
>>> record.dbxrefs
[]
>>> record.annotations
{}
>>> record.letter_annotations
{}
>>> record.features
[]
在這個例子中,我們的範例 FASTA 檔案來自 NCBI,他們對於 FASTA 行的格式有相當明確的慣例。這表示可以解析這些資訊,並提取例如 GI 號碼和 accession。然而,來自其他來源的 FASTA 檔案各有不同,因此這在一般情況下是不可行的。
來自 GenBank 檔案的 SeqRecord 物件
如同之前的範例,我們將檢視整個鼠疫耶爾森菌 (Yersinia pestis biovar Microtus) str. 91001 質體 pPCP1 的序列,該序列最初從 NCBI 下載,但這次是以 GenBank 檔案的形式呈現。同樣地,這個檔案包含在 Biopython 單元測試的 GenBank 資料夾中,或可於我們的網站線上取得 NC_005816.gb。
此檔案包含單一記錄(即只有一個 LOCUS 行),並以以下內容開頭:
LOCUS NC_005816 9609 bp DNA circular BCT 21-JUL-2008
DEFINITION Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete
sequence.
ACCESSION NC_005816
VERSION NC_005816.1 GI:45478711
PROJECT GenomeProject:10638
...
我們將再次使用 Bio.SeqIO
來讀取此檔案,且程式碼幾乎與上述用於 FASTA 檔案的程式碼相同(詳情請參閱序列輸入/輸出 章節)。
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=['Project:58037'])
>>> record.seq
Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG')
name
來自 LOCUS 行,而 id
則包含版本後綴。描述來自 DEFINITION 行。
>>> record.id
'NC_005816.1'
>>> record.name
'NC_005816'
>>> record.description
'Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence'
GenBank 檔案沒有任何逐字母的註解。
>>> record.letter_annotations
{}
大部分的註解資訊會記錄在 annotations
字典中,例如:
>>> len(record.annotations)
13
>>> record.annotations["source"]
'Yersinia pestis biovar Microtus str. 91001'
dbxrefs
清單會從任何 PROJECT 或 DBLINK 行填入。
>>> record.dbxrefs
['Project:58037']
最後,也許最有趣的是,功能表中的所有條目(例如基因或 CDS 功能)都會以 SeqFeature
物件的形式記錄在 features
清單中。
>>> len(record.features)
41
我們將在下一節 功能、位置和定位物件 中討論 SeqFeature
物件。
功能、位置和定位物件
SeqFeature 物件
序列功能是描述序列的重要部分。一旦您超越序列本身,就需要某種方法來組織和輕鬆取得關於序列的更多「抽象」資訊。雖然可能無法開發一個涵蓋所有內容的通用序列功能類別,但 Biopython 的 SeqFeature
類別嘗試盡可能地封裝關於序列的資訊。此設計主要基於 GenBank/EMBL 功能表,因此如果您了解它們的外觀,您可能會更容易掌握 Biopython 類別的結構。
每個 SeqFeature
物件的關鍵概念是描述父序列(通常是 SeqRecord
物件)上的區域。該區域使用位置物件描述,通常是兩個位置之間的範圍(請參閱下方位置和定位 章節)。
SeqFeature
類別有許多屬性,因此我們先列出這些屬性及其一般功能,然後在本章稍後逐步說明範例,展示這些屬性如何應用於真實範例。SeqFeature
的屬性如下:
- .type
這是功能類型的文字描述(例如,會是類似「CDS」或「gene」之類的文字)。
- .location
您正在處理的序列上
SeqFeature
的位置,請參閱下方位置和定位 章節。SeqFeature
將其大部分功能委派給位置物件,並包含許多用於位置屬性的快捷屬性。- .ref
是
.location.ref
的縮寫 – 位置所參考的任何(不同)參考序列。通常只是 None。- .ref_db
是
.location.ref_db
的縮寫 – 指定.ref
中任何識別碼所參考的資料庫。通常只是 None。- .strand
是
.location.strand
的縮寫 – 功能所在的序列鏈。對於雙鏈核苷酸序列,這可能是 \(1\) 表示頂鏈,\(-1\) 表示底鏈,\(0\) 表示鏈很重要但未知,或None
表示無關緊要。對於蛋白質或單鏈序列,則為 None。
- .qualifiers
這是一個 Python 字典,包含關於功能的額外資訊。鍵是關於值中包含的資訊的一些簡潔單字描述,而值是實際的資訊。例如,限定符的常見鍵可能是「evidence」,而值可能是「計算 (非實驗)」。這只是一種讓查看功能的人知道它尚未經實驗(即在實驗室中)確認的方法。請注意,其他的值將是一個字串清單(即使只有一個字串)。這是 GenBank/EMBL 檔案中功能表的反映。
- .sub_features
這曾經用來表示具有複雜位置的功能,例如 GenBank/EMBL 檔案中的「joins」。隨著
CompoundLocation
物件的引入,此功能已被棄用,現在應該忽略。
位置和定位
每個 SeqFeature
物件的關鍵概念是描述父序列上的區域,我們使用位置物件來描述,通常描述兩個位置之間的範圍。為了明確我們使用的術語,請看以下說明:
- position
指的是序列上的單一定位,可能是模糊的或明確的。例如,5、20、
<100
和>200
都是定位。- location
定位是由某些定位所界定的序列區域。例如,
5..20
(即 5 到 20)是一個定位。
我只是提一下,因為有時我會在這兩者之間感到困惑。
SimpleLocation 物件
除非您處理真核基因,否則大多數 SeqFeature
定位都非常簡單 - 您只需要起始和結束座標以及鏈。這基本上就是基本的 SimpleLocation
物件所做的一切。
當然,在實踐中,事情可能會更複雜。首先,我們必須處理由幾個區域組成的複合定位。其次,定位本身可能是模糊的(不精確)。
CompoundLocation 物件
Biopython 1.62 引入了 CompoundLocation
,作為重組由多個區域組成的複雜定位之表示方式的一部分。主要用途是處理 EMBL/GenBank 檔案中的「join」定位。
模糊定位
到目前為止,我們只使用簡單的定位。處理功能定位時的一個複雜問題來自於定位本身。在生物學中,很多時候事情並非完全確定(儘管我們實驗室生物學家盡力使它們確定!)。例如,您可能會進行二核苷酸引物實驗,並發現 mRNA 轉錄物的起始點位於兩個位點之一。這是非常有用的資訊,但複雜之處在於如何將其表示為定位。為了幫助我們處理這個問題,我們有了模糊定位的概念。基本上,有幾種類型的模糊定位,因此我們有五個類別來處理它們:
- ExactPosition
顧名思義,這個類別表示沿序列指定為精確的定位。這表示為一個數字,您可以查看物件的
position
屬性來取得定位。- BeforePosition
這個類別表示發生在指定位點之前的模糊定位。在 GenBank/EMBL 表示法中,這表示為類似
<13
的內容,表示實際位置位於小於 13 的某處。若要取得指定的上限,請查看物件的position
屬性。- AfterPosition
與
BeforePosition
相反,這個類別表示發生在指定位點之後的定位。在 GenBank 中表示為>13
,並且與BeforePosition
類似,您可以查看物件的position
屬性來取得邊界編號。- WithinPosition
此類別偶爾用於 GenBank/EMBL 定位,它模擬發生在兩個指定核苷酸之間的定位。在 GenBank/EMBL 表示法中,這會表示為
(1.5)
,表示定位位於 1 到 5 的範圍內。- OneOfPosition
此類別偶爾用於 GenBank/EMBL 定位,它處理存在多個可能值的定位,例如,如果起始密碼子不明確,並且有兩個候選基因起始位置,您可以使用此類別。或者,這可能明確地作為兩個相關的基因功能處理。
- UnknownPosition
這個類別處理位置不明的定位。這不會在 GenBank/EMBL 中使用,但與 UniProt 中使用的 '?' 功能座標相對應。
以下範例說明如何建立具有模糊端點的定位:
>>> from Bio import SeqFeature
>>> start_pos = SeqFeature.AfterPosition(5)
>>> end_pos = SeqFeature.BetweenPosition(9, left=8, right=9)
>>> my_location = SeqFeature.SimpleLocation(start_pos, end_pos)
請注意,Biopython 1.59 中某些模糊定位的詳細資訊有所變更,特別是對於 BetweenPosition 和 WithinPosition,您現在必須明確指出哪個整數位置應使用於切片等。對於起始位置,這通常是較低(左側)的值,而對於結束位置,這通常是較高(右側)的值。
如果您列印出 SimpleLocation
物件,您可以取得資訊的良好表示:
>>> print(my_location)
[>5:(8^9)]
我們可以使用定位的起始和結束屬性來存取模糊的起始和結束位置:
>>> my_location.start
AfterPosition(5)
>>> print(my_location.start)
>5
>>> my_location.end
BetweenPosition(9, left=8, right=9)
>>> print(my_location.end)
(8^9)
如果您不想處理模糊定位,而只想取得數字,它們實際上是整數的子類別,因此應該像整數一樣運作:
>>> int(my_location.start)
5
>>> int(my_location.end)
9
同樣地,為了方便建立位置而不用擔心模糊位置的問題,您可以直接將數字傳遞給 FeaturePosition
建構函式,然後會得到 ExactPosition
物件。
>>> exact_location = SeqFeature.SimpleLocation(5, 9)
>>> print(exact_location)
[5:9]
>>> exact_location.start
ExactPosition(5)
>>> int(exact_location.start)
5
以上是關於在 Biopython 中處理模糊位置的大部分細節。它的設計宗旨是讓處理模糊位置不會比處理精確位置複雜太多,希望您覺得確實如此!
位置測試
您可以將 Python 關鍵字 in
與 SeqFeature
或位置物件一起使用,以查看父坐標的鹼基/殘基是否在特徵/位置內。
例如,假設您有一個感興趣的 SNP,並且您想知道這個 SNP 在哪些特徵內,並且假設這個 SNP 的索引為 4350(Python 計數!)。這是一個簡單的暴力破解解決方案,我們只是在迴圈中逐一檢查所有特徵
>>> from Bio import SeqIO
>>> my_snp = 4350
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> for feature in record.features:
... if my_snp in feature:
... print("%s %s" % (feature.type, feature.qualifiers.get("db_xref")))
...
source ['taxon:229193']
gene ['GeneID:2767712']
CDS ['GI:45478716', 'GeneID:2767712']
請注意,來自 GenBank 或 EMBL 檔案並使用 joins 定義的基因和 CDS 特徵是外顯子的聯集,它們不涵蓋任何內含子。
特徵或位置描述的序列
SeqFeature
或位置物件不直接包含序列,而是位置(請參閱位置和位址章節)描述如何從父序列取得此序列。例如,考慮一個(短)基因序列,其位置在反向鏈上為 5:18
,在 GenBank/EMBL 符號中,使用從 1 開始的計數,則為 complement(6..18)
,如下所示
>>> from Bio.Seq import Seq
>>> from Bio.SeqFeature import SeqFeature, SimpleLocation
>>> seq = Seq("ACCGAGACGGCAAAGGCTAGCATAGGTATGAGACTTCCTTCCTGCCAGTGCTGAGGAACTGGGAGCCTAC")
>>> feature = SeqFeature(SimpleLocation(5, 18, strand=-1), type="gene")
您可以取得父序列,將其切片以提取 5:18
,然後取反向互補。特徵位置的起點和終點類似於整數,因此可以這樣做
>>> feature_seq = seq[feature.location.start : feature.location.end].reverse_complement()
>>> print(feature_seq)
AGCCTTTGCCGTC
這是一個簡單的範例,所以還不算太糟 – 但是一旦您必須處理複合特徵(joins),就會變得相當混亂。因此,SeqFeature
物件有一個 extract
方法來處理所有這些(並且自 Biopython 1.78 起,可以透過提供參照序列的字典來處理轉錄剪接)
>>> feature_seq = feature.extract(seq)
>>> print(feature_seq)
AGCCTTTGCCGTC
SeqFeature
或位置的長度與其描述的序列區域的長度相符。
>>> print(len(feature_seq))
13
>>> print(len(feature))
13
>>> print(len(feature.location))
13
對於 SimpleLocation
物件,長度只是起點和終點位置之間的差異。但是,對於 CompoundLocation
,長度是組成區域的總和。
比較
SeqRecord
物件可能非常複雜,但這是一個簡單的範例
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> record1 = SeqRecord(Seq("ACGT"), id="test")
>>> record2 = SeqRecord(Seq("ACGT"), id="test")
當您嘗試比較這些「相同」的記錄時會發生什麼事?
>>> record1 == record2
或許令人驚訝的是,舊版本的 Biopython 會將 Python 的預設物件比較用於 SeqRecord
,這表示只有在這些變數指向記憶體中的同一個物件時,record1 == record2
才會傳回 True
。在這個範例中,record1 == record2
在這裡會傳回 False
!
>>> record1 == record2 # on old versions of Biopython!
False
從 Biopython 1.67 開始,像 record1 == record2
這樣的 SeqRecord
比較將會明確引發錯誤,以避免人們被此問題所困擾
>>> record1 == record2
Traceback (most recent call last):
...
NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest.
相反地,您應該檢查您感興趣的屬性,例如識別符和序列
>>> record1.id == record2.id
True
>>> record1.seq == record2.seq
True
請注意,比較複雜的物件會很快變得複雜(另請參閱比較 Seq 物件章節)。
參考文獻
另一個與序列相關的常見註解是對處理該序列的期刊或其他已出版著作的參考。我們在 Biopython 中有一種相當簡單的方法來表示參考文獻 – 我們有一個 Bio.SeqFeature.Reference
類別,它將有關參考文獻的相關資訊儲存為物件的屬性。
這些屬性包括您在參考文獻中會看到的內容,例如 journal
、title
和 authors
。此外,它還可以保存 medline_id
和 pubmed_id
以及關於參考文獻的 comment
。這些都可以簡單地作為物件的屬性來存取。
參考文獻也有一個 location
物件,因此它可以指定參考文獻所參考的序列上的特定位置。例如,您可能有一本期刊正在處理位於 BAC 上的特定基因,並且想要指定它只準確地參考這個位置。location
是一個可能模糊的位置,如位置和位址章節中所述。
任何參考物件都作為列表儲存在 SeqRecord
物件的 annotations
字典中,其鍵為「references」。就這樣而已。參考文獻的目的是易於處理,並希望足夠通用以涵蓋許多使用案例。
format 方法
SeqRecord
類別的 format()
方法會給出一個字串,其中包含使用 Bio.SeqIO
支援的其中一種輸出檔案格式(例如 FASTA)格式化的記錄
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> record = SeqRecord(
... Seq(
... "MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD"
... "GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK"
... "NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM"
... "SSAC"
... ),
... id="gi|14150838|gb|AAK54648.1|AF376133_1",
... description="chalcone synthase [Cucumis sativus]",
... )
>>> print(record.format("fasta"))
應給出
>gi|14150838|gb|AAK54648.1|AF376133_1 chalcone synthase [Cucumis sativus]
MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD
GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK
NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM
SSAC
此 format
方法採用單一強制性引數,即 Bio.SeqIO
支援的小寫字串作為輸出格式(請參閱序列輸入/輸出章節)。但是,Bio.SeqIO
可以寫入的某些檔案格式需要多個記錄(通常適用於多序列比對格式),因此無法透過此 format()
方法運作。另請參閱將您的 SeqRecord 物件作為格式化的字串取得章節。
切片 SeqRecord
您可以切片 SeqRecord
,以取得涵蓋序列一部分的新 SeqRecord
。這裡重要的是,任何每個字母的註解也會被切片,並且會保留完全落在新序列中的任何特徵(並調整它們的位置)。
例如,採用先前使用的同一個 GenBank 檔案
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=['Project:58037'])
>>> len(record)
9609
>>> len(record.features)
41
在此範例中,我們將重點放在 pim
基因,YP_pPCP05
。如果您直接查看 GenBank 檔案,您會發現此基因/CDS 的位置字串為 4343..4780
,或在 Python 計數中為 4342:4780
。從查看檔案中,您可以計算出它們是檔案中的第十二個和第十三個條目,因此在 Python 從零開始的計數中,它們是 features
列表中的條目 \(11\) 和 \(12\)
>>> print(record.features[20])
type: gene
location: [4342:4780](+)
qualifiers:
Key: db_xref, Value: ['GeneID:2767712']
Key: gene, Value: ['pim']
Key: locus_tag, Value: ['YP_pPCP05']
>>> print(record.features[21])
type: CDS
location: [4342:4780](+)
qualifiers:
Key: codon_start, Value: ['1']
Key: db_xref, Value: ['GI:45478716', 'GeneID:2767712']
Key: gene, Value: ['pim']
Key: locus_tag, Value: ['YP_pPCP05']
Key: note, Value: ['similar to many previously sequenced pesticin immunity protein entries of Yersinia pestis plasmid pPCP, e.g. gi| 16082683|,ref|NP_395230.1| (NC_003132) , gi|1200166|emb|CAA90861.1| (Z54145 ) , gi|1488655| emb|CAA63439.1| (X92856) , gi|2996219|gb|AAC62543.1| (AF053945) , and gi|5763814|emb|CAB531 67.1| (AL109969)']
Key: product, Value: ['pesticin immunity protein']
Key: protein_id, Value: ['NP_995571.1']
Key: transl_table, Value: ['11']
Key: translation, Value: ['MGGGMISKLFCLALIFLSSSGLAEKNTYTAKDILQNLELNTFGNSLSHGIYGKQTTFKQTEFTNIKSNTKKHIALINKDNSWMISLKILGIKRDEYTVCFEDFSLIRPPTYVAIHPLLIKKVKSGNFIVVKEIKKSIPGCTVYYH']
讓我們將此父記錄從 4300 切片到 4800(足以包含 pim
基因/CDS),並查看我們獲得多少特徵
>>> sub_record = record[4300:4800]
>>> sub_record
SeqRecord(seq=Seq('ATAAATAGATTATTCCAAATAATTTATTTATGTAAGAACAGGATGGGAGGGGGA...TTA'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])
>>> len(sub_record)
500
>>> len(sub_record.features)
2
我們的子記錄只有兩個特徵,即 YP_pPCP05
的基因和 CDS 條目
>>> print(sub_record.features[0])
type: gene
location: [42:480](+)
qualifiers:
Key: db_xref, Value: ['GeneID:2767712']
Key: gene, Value: ['pim']
Key: locus_tag, Value: ['YP_pPCP05']
>>> print(sub_record.features[1])
type: CDS
location: [42:480](+)
qualifiers:
Key: codon_start, Value: ['1']
Key: db_xref, Value: ['GI:45478716', 'GeneID:2767712']
Key: gene, Value: ['pim']
Key: locus_tag, Value: ['YP_pPCP05']
Key: note, Value: ['similar to many previously sequenced pesticin immunity protein entries of Yersinia pestis plasmid pPCP, e.g. gi| 16082683|,ref|NP_395230.1| (NC_003132) , gi|1200166|emb|CAA90861.1| (Z54145 ) , gi|1488655| emb|CAA63439.1| (X92856) , gi|2996219|gb|AAC62543.1| (AF053945) , and gi|5763814|emb|CAB531 67.1| (AL109969)']
Key: product, Value: ['pesticin immunity protein']
Key: protein_id, Value: ['NP_995571.1']
Key: transl_table, Value: ['11']
Key: translation, Value: ['MGGGMISKLFCLALIFLSSSGLAEKNTYTAKDILQNLELNTFGNSLSHGIYGKQTTFKQTEFTNIKSNTKKHIALINKDNSWMISLKILGIKRDEYTVCFEDFSLIRPPTYVAIHPLLIKKVKSGNFIVVKEIKKSIPGCTVYYH']
請注意,它們的位置已調整以反映新的父序列!
雖然 Biopython 對於特徵(以及任何每個字母的註解)已做了一些合理且希望直觀的事情,但對於其他註解,不可能知道這是否仍然適用於子序列。為了避免猜測,除了分子類型之外,.annotations
和 .dbxrefs
會從子記錄中省略,並且由您自行決定是否酌情傳輸任何相關資訊。
>>> sub_record.annotations
{'molecule_type': 'DNA'}
>>> sub_record.dbxrefs
[]
您可能希望保留其他條目,例如生物體?請注意,複製整個註解字典時,您的部分序列不再是環狀 DNA – 它現在是線性的
>>> sub_record.annotations["topology"] = "linear"
關於記錄的 id
、name
和 description
也可以提出相同的觀點,但為了實用起見,這些都保留了
>>> sub_record.id
'NC_005816.1'
>>> sub_record.name
'NC_005816'
>>> sub_record.description
'Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence'
不過,這很好地說明了這個問題,我們新的子記錄不是質體的完整序列,所以描述是錯誤的!讓我們修正這個問題,然後使用format 方法章節中上述的 format
方法,將子記錄視為簡化的 GenBank 檔案
>>> sub_record.description = (
... "Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, partial"
... )
>>> print(sub_record.format("genbank")[:200] + "...")
LOCUS NC_005816 500 bp DNA linear UNK 01-JAN-1980
DEFINITION Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, partial.
ACCESSION NC_005816
VERSION NC_0058...
新增 SeqRecord 物件
您可以將 SeqRecord
物件加在一起,以產生新的 SeqRecord
。這裡重要的是,也會新增任何常見的每個字母的註解,所有特徵都會被保留(並調整它們的位置),並且也會保留任何其他常見的註解(例如 id、名稱和描述)。
為了舉例說明每個字母的註解,我們將使用 FASTQ 檔案中的第一個記錄。 序列輸入/輸出章節將說明 SeqIO
函數
>>> from Bio import SeqIO
>>> record = next(SeqIO.parse("example.fastq", "fastq"))
>>> len(record)
25
>>> print(record.seq)
CCCTTCTTGTCTTCAGCGTTTCTCC
>>> print(record.letter_annotations["phred_quality"])
[26, 26, 18, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 22, 26, 26, 26, 26, 26, 26, 26, 23, 23]
假設這是 Roche 454 的數據,並且根據其他資訊,您認為 TTT
應該只有 TT
。我們可以先在「額外」的第三個 T
前後切割 SeqRecord
,來建立一個新的編輯過的紀錄。
>>> left = record[:20]
>>> print(left.seq)
CCCTTCTTGTCTTCAGCGTT
>>> print(left.letter_annotations["phred_quality"])
[26, 26, 18, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 22, 26, 26, 26, 26]
>>> right = record[21:]
>>> print(right.seq)
CTCC
>>> print(right.letter_annotations["phred_quality"])
[26, 26, 23, 23]
現在將兩個部分加在一起。
>>> edited = left + right
>>> len(edited)
24
>>> print(edited.seq)
CCCTTCTTGTCTTCAGCGTTCTCC
>>> print(edited.letter_annotations["phred_quality"])
[26, 26, 18, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 22, 26, 26, 26, 26, 26, 26, 23, 23]
簡單又直觀吧?我們希望如此!您可以用更簡短的方式完成:
>>> edited = record[:20] + record[21:]
現在,舉一個包含 features 的例子,我們將使用 GenBank 檔案。假設您有一個環狀基因組:
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=['Project:58037'])
>>> len(record)
9609
>>> len(record.features)
41
>>> record.dbxrefs
['Project:58037']
>>> record.annotations.keys()
dict_keys(['molecule_type', 'topology', 'data_file_division', 'date', 'accessions', 'sequence_version', 'gi', 'keywords', 'source', 'organism', 'taxonomy', 'references', 'comment'])
您可以像這樣移動原點:
>>> shifted = record[2000:] + record[:2000]
>>> shifted
SeqRecord(seq=Seq('GATACGCAGTCATATTTTTTACACAATTCTCTAATCCCGACAAGGTCGTAGGTC...GGA'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])
>>> len(shifted)
9609
請注意,這並非完美,因為某些註釋(如資料庫交叉引用)、除了分子類型之外的所有註釋,以及其中一個 features(來源 feature)都已遺失。
>>> len(shifted.features)
40
>>> shifted.dbxrefs
[]
>>> shifted.annotations.keys()
dict_keys(['molecule_type'])
這是因為 SeqRecord
切割步驟對於保留哪些註釋非常謹慎(錯誤地傳播註釋可能會導致重大問題)。如果您想要保留資料庫交叉引用或註釋字典,則必須明確地執行此操作。
>>> shifted.dbxrefs = record.dbxrefs[:]
>>> shifted.annotations = record.annotations.copy()
>>> shifted.dbxrefs
['Project:58037']
>>> shifted.annotations.keys()
dict_keys(['molecule_type', 'topology', 'data_file_division', 'date', 'accessions', 'sequence_version', 'gi', 'keywords', 'source', 'organism', 'taxonomy', 'references', 'comment'])
另請注意,在像這樣的範例中,您可能應該更改紀錄識別符,因為 NCBI 參考的是原始未修改的序列。
反向互補 SeqRecord 物件
Biopython 1.57 的新功能之一是 SeqRecord
物件的 reverse_complement
方法。此方法試圖在易用性與擔心如何處理反向互補紀錄中的註釋之間取得平衡。
對於序列,此方法使用 Seq 物件的反向互補方法。任何 features 都會被轉移,並且會重新計算位置和鏈。同樣地,任何每個字母的註釋也會被複製但反轉(對於品質分數等典型範例來說,這是合理的)。但是,大部分註釋的轉移是有問題的。
例如,如果紀錄 ID 是一個登錄號,則該登錄號實際上不應適用於反向互補序列,並且預設轉移識別符很容易在下游分析中導致細微的資料損毀。因此,預設情況下,SeqRecord
的 id、name、description、annotations 和資料庫交叉引用都不會被轉移。
SeqRecord
物件的 reverse_complement
方法採用多個與紀錄屬性對應的可選參數。將這些參數設定為 True
表示複製舊值,而 False
表示捨棄舊值並使用預設值。您也可以提供新的期望值來代替。
考慮以下範例紀錄:
>>> from Bio import SeqIO
>>> rec = SeqIO.read("NC_005816.gb", "genbank")
>>> print(rec.id, len(rec), len(rec.features), len(rec.dbxrefs), len(rec.annotations))
NC_005816.1 9609 41 1 13
在這裡,我們取得反向互補,並指定一個新的識別符 – 但請注意,大多數註釋都已捨棄(但 features 沒有)。
>>> rc = rec.reverse_complement(id="TESTING")
>>> print(rc.id, len(rc), len(rc.features), len(rc.dbxrefs), len(rc.annotations))
TESTING 9609 41 0 0