Pyfastx:一個快速隨機讀取基因組數據的Python模塊

2021-02-20 簡說基因

今天介紹一個同門師兄開發的 Python 模塊：pyfastx，用於快速隨機訪問基因組序列文件。作品發表在生信頂刊上，必須強行安利一波。

師兄任職於成都大學，專注於生物信息學研究，是真正的Too maker。

Pyfastx: a robust python module for fast random access to sequences from plain and gzipped FASTA/Q file

文章地址：https://doi.org/10.1093/bib/bbaa368

代碼地址：https://github.com/lmdu/pyfastx

模塊的使用文檔非常詳細，感興趣的朋友可以直接到 pyfastx 官網查看使用說明。

本文僅作拋磚引玉，首先我們來看一下 pyfastx 的特點。

功能很多，覆蓋了平時序列文件操作的常見需求。Pyfastx 內部含有多個功能模塊，比如:

FASTX 接口，為迭代 Fasta/q 文件提供統一的接口FASTA 接口，迭代或隨機訪問 Fasta 文件FASTQ 接口 ，迭代或隨機訪問 Fastq 文件Sequence 類，提供 Fasta 記錄的常用操作安裝

目前，pyfastx 支持 Python 3.5 以上的版本，通過pip即可安裝。

pip install pyfastx
FASTX 模塊FASTA 文件迭代
迭代 Fasta 文件時，返回一個元組（name, seq, comment），其中 comment 是標題欄第一個空格後面的內容。
>>> fa = pyfastx.Fastx('tests/data/test.fa.gz')
>>> for name,seq,comment in fa:
>>>     print(name)
>>>     print(seq)
>>>     print(comment)

>>> #always output uppercase sequence
>>> for item in pyfastx.Fastx('tests/data/test.fa', uppercase=True):
>>>     print(item)

>>> #Manually specify sequence format
>>> for item in pyfastx.Fastx('tests/data/test.fa', format="fasta"):
>>>     print(item)
FASTQ 文件迭代
迭代 Fastq 文件時，返回一個元組（name, seq, qual, comment），其中 comment 是標題欄第一個空格後面的內容。
>>> fq = pyfastx.Fastx('tests/data/test.fq.gz')
>>> for name,seq,qual,comment in fq:
>>>     print(name)
>>>     print(seq)
>>>     print(qual)
>>>     print(comment)
FASTA 模塊
讀取 Fasta 文件，並且支持隨機訪問其中的任意序列。
這裡要說明一下順序迭代和隨機讀取的區別。順序迭代顧名思義就是從一個文件的開始逐條記錄往後讀，直至最後一條記錄。
隨機讀取就是能夠直接訪問指定的序列，不需要從頭讀到尾。怎麼實現呢？一般是先要為序列建立索引，pyfastx 也不例外。
FASTA 文件讀取
>>> import pyfastx
>>> fa = pyfastx.Fasta('test/data/test.fa.gz')
>>> fa
<Fasta> test/data/test.fa.gz contains 211 seqs
FASTA 文件迭代
Fasta 文件中每條序列最重要的就是名稱和序列信息了，這兩個信息可以方便地通過迭代返回。
>>> import pyfastx
>>> for name, seq in pyfastx.Fasta('test/data/test.fa.gz', build_index=False):
>>>     print(name, seq)
此外也可以直接返回 FASTA 對象。
>>> import pyfastx
>>> for seq in pyfastx.Fasta('test/data/test.fa.gz'):
>>>     print(seq.name)
>>>     print(seq.seq)
>>>     print(seq.description)
FASTA 類
FASTA 對象有許多屬性和方法可供使用，如計算 GC 含量、計算 N50/L50、提取任意序列等。
以提取指定序列為例，FASTA 不僅可以提取指定序列，還可以指定序列的某一區間。
>>> # get subsequence with start and end position
>>> interval = (1, 10)
>>> fa.fetch('JZ822577.1', interval)
'CTCTAGAGAT'

>>> # get subsequences with a list of start and end position
>>> intervals = [(1, 10), (50, 60)]
>>> fa.fetch('JZ822577.1', intervals)
'CTCTAGAGATTTTAGTTTGAC'

>>> # get subsequences with reverse strand
>>> fa.fetch('JZ822577.1', (1, 10), strand='-')
'ATCTCTAGAG'
當然，FASTA 對象還有更加花式的序列提取方式。
>>> # get sequence like a dictionary by identifier
>>> s1 = fa['JZ822577.1']
>>> s1
<Sequence> JZ822577.1 with length of 333

>>> # get sequence like a list by index
>>> s2 = fa[2]
>>> s2
<Sequence> JZ822579.1 with length of 176

>>> # get last sequence
>>> s3 = fa[-1]
>>> s3
<Sequence> JZ840318.1 with length of 134

>>> # check a sequence name weather in FASTA file
>>> 'JZ822577.1' in fa
True
Sequence 類
FASTA 文件的記錄，被封裝成 Sequence 類，為操作 FASTA 文件的記錄提供接口。例如，你可以對它進行切片。
>>> # get a sub seq from sequence
>>> s = fa[-1]
>>> ss = s[10:30]
>>> ss
<Sequence> JZ840318.1 from 11 to 30

>>> ss.name
'JZ840318.1:11-30'

>>> ss.seq
'CTTCTTCCTGTGGAAAGTAA'

>>> ss = s[-10:]
>>> ss
<Sequence> JZ840318.1 from 125 to 134

>>> ss.name
'JZ840318.1:125-134'

>>> ss.seq
'CCATGTTGGT'
或者對 Sequence 進行反向、互補或反向互補。
>>> # get sliced sequence
>>> fa[0][10:20].seq
'GTCAATTTCC'

>>> # get reverse of sliced sequence
>>> fa[0][10:20].reverse
'CCTTTAACTG'

>>> # get complement of sliced sequence
>>> fa[0][10:20].complement
'CAGTTAAAGG'

>>> # get reversed complement sequence, corresponding to sequence in antisense strand
>>> fa[0][10:20].antisense
'GGAAATTGAC'
FASTQ 模塊FASTQ 文件讀取
讀取 Fastq 文件，並支持隨機訪問，前提是先要構建索引。
>>> import pyfastx
>>> fq = pyfastx.Fastq('tests/data/test.fq.gz')
>>> fq
<Fastq> tests/data/test.fq.gz contains 100 reads
FASTQ 文件迭代
Fastq 每一條記錄有 4 行，其中 comment 通常總為+號，因此有價值的是name, seq, qual三項信息。
>>> import pyfastx
>>> for name,seq,qual in pyfastx.Fastq('tests/data/test.fq.gz', build_index=False):
>>>     print(name)
>>>     print(seq)
>>>     print(qual)
也可以直接返回 FASTQ 對象。
>>> import pyfastx
>>> for read in pyfastx.Fastq('test/data/test.fq.gz'):
>>>     print(read.name)
>>>     print(read.seq)
>>>     print(read.qual)
>>>     print(read.quali)
FASTQ 類
FASTQ 類封裝了 Fastq 文件，為 Fastq 文件常用操作提供方便的接口。
>>> # get read counts in FASTQ
>>> len(fq)
800

>>> # get total bases
>>> fq.size
120000

>>> # get GC content of FASTQ file
>>> fq.gc_content
66.17471313476562

>>> # get composition of bases in FASTQ
>>> fq.composition
{'A': 20501, 'C': 39705, 'G': 39704, 'T': 20089, 'N': 1}

>>> # New in pyfastx 0.6.10
>>> # get average length of reads
>>> fq.avglen
150.0

>>> # get maximum lenth of reads
>>> fq.maxlen
150

>>> # get minimum length of reas
>>> fq.minlen
150

>>> # get maximum quality score
>>> fq.maxqual
70

>>> # get minimum quality score
>>> fq.minqual
35

>>> # get phred which affects the quality score conversion
>>> fq.phred
33

>>> # Guess fastq quality encoding system
>>> # New in pyfastx 0.4.1
>>> fq.encoding_type
['Sanger Phred+33', 'Illumina 1.8+ Phred+33']
Read 類
FASTQ 文件的每一條記錄是一個 read，Read 類為操作 Fastq 記錄提供了接口。
獲取 Read 對象
>>> #get read like a dict by read name
>>> r1 = fq['A00129:183:H77K2DMXX:1:1101:4752:1047']
>>> r1
<Read> A00129:183:H77K2DMXX:1:1101:4752:1047 with length of 150

>>> # get read like a list by index
>>> r2 = fq[10]
>>> r2
<Read> A00129:183:H77K2DMXX:1:1101:18041:1078 with length of 150

>>> # get the last read
>>> r3 = fq[-1]
>>> r3
<Read> A00129:183:H77K2DMXX:1:1101:31575:4726 with length of 150

>>> # check a read weather in FASTQ file
>>> 'A00129:183:H77K2DMXX:1:1101:4752:1047' in fq
True
獲取 Read 對象的信息。
>>> r = fq[-10]
>>> r
<Read> A00129:183:H77K2DMXX:1:1101:1750:4711 with length of 150

>>> # get read order number in FASTQ file
>>> r.id
791

>>> # get read name
>>> r.name
'A00129:183:H77K2DMXX:1:1101:1750:4711'

>>> # get read full header line, New in pyfastx 0.6.3
>>> r.description
'@A00129:183:H77K2DMXX:1:1101:1750:4711 1:N:0:CAATGGAA+CGAGGCTG'

>>> # get read length
>>> len(r)
150

>>> # get read sequence
>>> r.seq
'CGAGGAAATCGACGTCACCGATCTGGAAGCCCTGCGCGCCCATCTCAACCAGAAATGGGGTGGCCAGCGCGGCAAGCTGACCCTGCTGCCGTTCCTGGTCCGCGCCATGGTCGTGGCGCTGCGCGACTTCCCGCAGTTGAACGCGCGCTA'

>>> # get raw string of read, New in pyfastx 0.6.3
>>> print(r.raw)
@A00129:183:H77K2DMXX:1:1101:1750:4711 1:N:0:CAATGGAA+CGAGGCTG
CGAGGAAATCGACGTCACCGATCTGGAAGCCCTGCGCGCCCATCTCAACCAGAAATGGGGTGGCCAGCGCGGCAAGCTGACCCTGCTGCCGTTCCTGGTCCGCGCCATGGTCGTGGCGCTGCGCGACTTCCCGCAGTTGAACGCGCGCTA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FF,FFFFFFFFFFFFFFFFFFFFFFFFFF,F:FFFFFFFFF:

>>> # get read quality ascii string
>>> r.qual
'FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FF,FFFFFFFFFFFFFFFFFFFFFFFFFF,F:FFFFFFFFF:'

>>> # get read quality integer value, ascii - 33 or 64
>>> r.quali
[37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 25, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 25, 37, 37, 11, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 11, 37, 25, 37, 37, 37, 37, 37, 37, 37, 37, 37, 25]

>>> # get read length
>>> len(r)
150
命令行接口
$ pyfastx -h

usage: pyfastx COMMAND [OPTIONS]

A command line tool for FASTA/Q file manipulation

optional arguments:
  -h, --help     show this help message and exit
  -v, --version  show program's version number and exit

Commands:

    index        build index for fasta/q file
    stat         show detailed statistics information of fasta/q file
    split        split fasta/q file into multiple files
    fq2fa        convert fastq file to fasta file
    subseq       get subsequences from fasta file by region
    sample       randomly sample sequences from fasta or fastq file
    extract      extract full sequences or reads from fasta/q file
Pyfastx 還提供了額外接口，方便在命令行下直接操作 Fasta/q 文件。好了，鑑於時間和篇幅有限，這部分大家自行探索吧。
最後，Pyfastx 雖然已經發表文章，但還在持續維護更新。希望大家多多使用，有什麼建議可以跟作者反饋。
好的工具和用戶是共同成長的，祝大家科研順利。
作者簡介：杜聯明，博士，成都大學特聘副研究員。畢業於四川大學生物信息學專業，主要從事基因組學與轉錄組學研究，生物信息學資料庫和數據分析軟體開發。相關成果發表在 Bioinformatics，Molecular Ecology Resources，Aging，Journal of Heredity，Genome Research，Molecular Biology and Evolution 等國際知名期刊上。

相關焦點

  • python海量數據快速查詢的技巧
    對於小文件而言，這樣的操作編碼簡單，運行速度也比較滿意，但是對於大型資料庫而言，將資料庫存為字典這個動作是非常耗費時間的，而且每次運行代碼都要執行這樣的操作，導致效率大大降低。想要改善這一狀況，有以下兩種解決辦法1.
  • Python數據分析:pandas讀取和寫入數據
    pandas讀取這些數據文件的方法如表格所示：01讀取寫入文本文件read_csv()方法用來讀取 csv格式的數據文件，read_table()方法則是讀取通用分隔符分隔的數據文件，它們的參數相同。使用cat命令顯示文件內容：import pandas as pddf = pd.read_csv('01.csv')df當使用read_table()時，運行代碼後出現一個Warning，使用的是最新的版本python3.7。
  • python隨機模塊22個函數詳解(上)
    作者：小伍哥來源： AI入門學習今天給大家介紹下python中的隨機模塊，隨機數可以用於數學，遊戲，安全等領域中，還經常被嵌入到算法中，用以提高算法效率，並提高程序的安全性。平時數據分析各種分布的數據構造也會用到。
  • Python讀取ini配置文件
    python看過我之前文章的同學可能知道，最近一直在做百度語音合成的功能，進一步的延伸功能，此次是批量生成文章的語音文件。格式如下：；注釋說明此文件應用場景[DATABASE]host = 127port = 3306[TYPE]cat = 0我們簡單的寫兩個配置參數信息，下面來看一下如何讀取信息。
  • biopython簡介
    biopython和bioperl, biojava項目類似，都是Open Bioinformatics Foundation組織的項目之一，旨在提供一個編程接口，方便生物信息數據的處理。
  • Python基礎教程(一) - 快速入門
    /usr/bin/python為Linux系統下Python解釋器的路徑，通常python解釋器的路徑安裝在/usr/local/bin或/usr/bin目錄下。程序輸入和raw_input()內建函數從用戶得到數據輸入的最好方式使用raw_input()函數，它讀取標準輸入，並將讀取到的數據賦值給指定的變量。
  • 如何序列化處理數據?真的會用pickle嗎?python中常用模塊詳解
    在編程中我們經常會遇到這樣的需求：將一串數據或者有效信息暫存在本地(一般情況下這種數據量並不大，當數據量較大時，可以選擇使用temp文件存儲方式，參見一個項目引發的思考，如何操作臨時文件？Python中的內置模塊實現)。
  • 使用biopython可視化染色體和基因元件
    基因組結構元件的可視化有多種方式，比如IGV等基因組瀏覽器中以track為單位的展示形式，亦或以circos為代表的圈圖形式，比如在細胞器基因組組裝中，基因元件常用圈圖形式展示，示例如下
  • 基於python+OpenCV模塊的人臉識別定位技術
    本文將基於OpenCV模塊，在windows作業系統上，利用python語言，進行人臉識別技術的研究。當然OpenCV的應用領域很廣，除了人臉識別之外，它還支持圖像分割、動作識別、視頻處理等技術。代碼分析下面我們對代碼進行分析，代碼如圖所示：一共不超過15行，當然這是建立在別人已有的數據上做的，如果自己寫的話，不會這麼簡單，我們這只是調用了別人的接口，而這個接口是開源的，共享的。代碼第1行導入opencv模塊。
  • python sys模塊的常見用法匯總
    python的內置模塊sys，提供了系統相關的一些變量和函數，在實際開發中，常見的有以下幾種用法1.6. sys.path該變量存儲了python尋找模塊的路徑7. sys.module該變量是一個字典，存儲了已經導入的模塊
  • Python數據讀取之生成器(generator)
    簡單來說，生成器是一個函數，它返回一個我們可以迭代的對象（迭代器），迭代器一次返回一個值較使用列表將所有數據都加載到內存中，生成器節省了大量內存空間。深度學習的數據讀取部分一般都需要使用迭代器。，而g是一個generator。
  • Python學習120課 pandas簡介kaggle下載數據及pandas讀取外部數據
    【每天幾分鐘，從零入門python編程的世界！】numpy的基本的東西我們學習差不多了，後面具體應用中遇到問題具體分析，然後去深入了解遇到的新的知識點就行。現在我們開始學習pandas，pandas一般用的更多，pandas是基於numpy去寫的。pandas是一個專門做數據結構和數據分析的庫。
  • Python視頻教程網課編程零基礎入門數據分析網絡爬蟲全套Python...
    因篇幅有限,以下展示的只是課程裡部分內容如對python課程有更多疑問 請諮詢客服 1零基礎入門全能班 01 –python簡介 02 第一個程序 03-python執行方式和pycharm設置 04-程序的注釋和算術運算符 05 程序執行原理 06變量的使用以及類型 07
  • 數據分析-numpy庫快速了解
    1.numpy是什麼庫NumPy是一個開源的Python科學計算基礎庫，包含： 一個強大的N維數組對象 ndarray 廣播功能函數 整合C/C++/Fortran代碼的工具 線性代數、傅立葉變換、隨機數生成等功能
  • 玩轉Python 中的隨機數
    隨機生成 0 到 1 之間的浮點數random.random() 方法會返回 [0.0, 1.0) 之間的浮點數，注意，這是一個左閉右開的區間，隨機數可能會是 0 但不可能為 1 。隨機生成 a 與 b 之間的整數使用 random.randint(a , b) 方法，你可以生成一個 a 與 b 之間的隨機整數，也就是 [a, b] 。
  • python隨機函數random分配應用,隨機分配8名老師到3個教室中
    羽憶教程最近遇到一個問題，要隨機分配8名老師到3個辦公室中，這時小編想要了python中的隨機函數random來進行分配工作，感覺小編像個月老一樣。python隨機函數python隨機函數在python中，想要生成隨機數
  • Python的武器庫05:numpy模塊(下)
    說到程式語言python，有一個著名的格言"餘生太短，只用python"。如果要分析為什麼會存在這麼一句格言？python的語法並不簡單，有複雜難懂的部分，之所以有這樣一句格言，是因為python中有很多強大的模塊，就像一個武器庫。
  • 數據科學的Python軟體包
    平臺無關Python可用於包括window，mac，Linux和Unix在內的各種平臺，因此一次編寫的代碼可以在另一個平臺上運行而無需進行任何更改。巨大的社區支持Python具有廣泛的社區支持之一，在dev op社區上存在各種活躍的論壇，python開發人員在該論壇上發布他們的錯誤，而社區則試圖幫助他們。
  • python安全開發軍規之四:使用安全的隨機數生成器
    背景日常開發中，必然會碰到需要生成隨機數的需求，比如生成圖片驗證碼，簡訊驗證碼……隨機數生成既然是這麼簡單的一個功能，開發必然也很簡單，我們看看怎麼生成一個隨機數，這裡以隨機生成1-100的整數為例。普通程式設計師的寫法import randomrandom.randint(1,100)只用了兩行代碼，程式設計師小Z就寫出了一個隨機數。
  • 使用樹莓派(Raspberry Pi)的真正隨機數生成器
    使用電視上的靜態信號將Raspberry Pi變成了真正的隨機數生成器。在國外，我們不再接收模擬地面廣播，因此在電視上查找靜態信號就像將其放在模擬頻道上一樣簡單。 我使用的設置是插入Raspberry Pi的eSecure USB 8MP網絡攝像頭，我將其指向電視。