轉(zhuǎn)載

引言

參考基因組是現(xiàn)代遺傳學(xué)的基石，更是生物信息學(xué)研究的基礎(chǔ)。不論你研究的是DNA、RNA，還是蛋白，你都無法避開參考基因組而進(jìn)行你的研究。在生物信息分析過程中，參考基因組最常用的一個(gè)場景就是序列比對---將測序的read比對回參考基因組上。那么，如果我們想知道參考基因組某個(gè)位置或者某個(gè)區(qū)域的堿基序列，我們應(yīng)該如何快速的解決這個(gè)問題呢？本文將介紹三種解決方案。

參考基因組格式

參考基因組通常以FASTA格式存在。在生物信息學(xué)中，F(xiàn)ASTA格式是一種用于記錄核酸序列或肽序列的文本格式，其中的核酸或氨基酸均以單個(gè)字母編碼呈現(xiàn)。該格式同時(shí)還允許在序列之前定義名稱和編寫注釋。這一格式最初由FASTA軟件包定義，但現(xiàn)今已是生物信息學(xué)領(lǐng)域的一項(xiàng)標(biāo)準(zhǔn)。（維基百科）

FASTA文件中存儲(chǔ)的一條完整的序列，由兩部分組成：

1. 開頭的單行描述行

描述行以“>”開頭，用以和數(shù)據(jù)行分開，后緊接的內(nèi)容為該序列的標(biāo)識(shí)符，該行剩余部分則為序列的描述（標(biāo)識(shí)符與描述均非必須）

2. 緊隨其后的序列行

序列行可以是多行，序列的結(jié)束以下一條序列的“>”出現(xiàn)為標(biāo)識(shí)。

FASTA格式示例

提取位置堿基

接下來我們以chr5:8397384-8397384為例，分別介紹三種方法提取某個(gè)位置的堿基。

01?利用samtools faidx

samtools faidx常常用來對參考基因組建立索引，但它還有個(gè)鮮為人知的功能，就是序列提取，如下：

02 利用bedtools getfasta

bedtools說明文檔中對getfasta的描述是“Extract DNA sequences into a fasta file based on feature coordinates.”顯而易見，bedtools getfasta的功能就是根據(jù)坐標(biāo)信息提取序列信息。操作如下：

bedtools getfasta有三個(gè)必選參數(shù)：-fi即參考基因組fasta文件；-bed即需要提取的位置坐標(biāo)信息，格式：chr\tstart\tend；-fo：輸出文件。

有一點(diǎn)需要說明，bedtools接收的是bed文件，而bed文件時(shí)0-based，因此我們要提取chr5:8397384-8397384位置的堿基，輸入文件pos.bed中需要寫入：chr5\t8397383\t8397384

03 利用pysam模塊

直接上代碼吧！

性能比較

三種方法耗時(shí)如下：

可以看出，samtools和bedtools的性能很好，python的性能就比較尷尬了。其實(shí)，個(gè)人比較推薦bedtools，比較容易進(jìn)行批量處理，把想處理的位置信息寫到輸出文件，然后就可以輕松的進(jìn)行序列提取。

版權(quán)聲明：本文為CSDN博主「whenfree」的原創(chuàng)文章，遵循CC 4.0 BY-SA版權(quán)協(xié)議，轉(zhuǎn)載請附上原文出處鏈接及本聲明。

原文鏈接：https://blog.csdn.net/whenfree/article/details/85305616