我想從 log(n)
時間複雜度(大O標記)或更少。一條記錄被定義為相當於FASTQ格式的四行。記錄不適合RAM,需要存儲在磁盤上。理想情況下,我想以壓縮格式存儲讀段。
我更喜歡不需要任何額外文件(例如參考基因組)的解決方案。
標題這個問題僅提及FASTQ,因為FASTQ是用於在磁盤上存儲未對齊讀取的通用格式。我對要求按時間複雜度順序 n
將數據進行一次有限的轉換為另一種文件格式的答案感到滿意。
Update
一個說明:希望以概率 1 / n
選擇隨機記錄。
我想從 log(n)
時間複雜度(大O標記)或更少。一條記錄被定義為相當於FASTQ格式的四行。記錄不適合RAM,需要存儲在磁盤上。理想情況下,我想以壓縮格式存儲讀段。
我更喜歡不需要任何額外文件(例如參考基因組)的解決方案。
標題這個問題僅提及FASTQ,因為FASTQ是用於在磁盤上存儲未對齊讀取的通用格式。我對要求按時間複雜度順序 n
將數據進行一次有限的轉換為另一種文件格式的答案感到滿意。
一個說明:希望以概率 1 / n
選擇隨機記錄。
要在恆定時間內獲取隨機記錄,在恆定時間內獲取任意記錄就足夠了。
我在這裡有兩種解決方案:一種帶有 tabix
,另一種帶有 grabix
。我認為 grabix
解決方案更優雅,但我保留下面的 tabix
解決方案,因為 tabix
是比 grabix更成熟的工具
。
感謝 user172818建議 grabix
。
此答案先前指出, tabix
和 grabix
在 log(n)
時間內執行查找。在仔細研究了grapix的源代碼和tabix文件之後,我現在確信查找的複雜性與 n
無關。但是,這兩個工具都使用一個索引,該索引的大小與 n
成比例。因此,索引的加載順序為 n
。但是,如果我們認為索引的加載是“ ...將數據有限地轉換為另一種文件格式...”,那麼我認為這一答案仍然是有效的。如果要檢索多個記錄,則需要將索引存儲在內存中,也許使用諸如pysam或htslib之類的框架。
grabix
bgzip
壓縮。 grabix
在bash中:
gzip -dc input.fastq.gz | bgzip -c > output.fastq.gzgrabix索引output.fastq.gz#在log(n)時間中檢索第5條記錄(基於1)#需要一些數學運算才能轉換索引(4 * 4 + 1、4 * 4 + 4)=(17,20)grabix抓取output.fastq.gz 17 20#計算該問題第二部分的記錄數export N_LINES = $(gzip -dc output.fastq.gz | wc -l)
tabix
tabix代碼更加複雜,並且基於這樣一個錯誤的假設: \ t
是FASTQ記錄中替換 \ n
的可接受字符。如果您對接近但不完全是FASTQ的文件格式感到滿意,則可以執行以下操作:
bgzip
壓縮。 tabix
請注意,我們需要刪除由 nl
引入的前導空格,並且需要引入虛擬染色體列以保持 tabix
happy:
gzip -dc input.fastq.gz |粘貼----| nl | sed's / ^ * //'| sed's / ^ / dummy \ t /'| bgzip -c > output.fastq.gztabix -s 1 -b 2 -e 2 output.fastq.gz#現在檢索log(n)timetabix output中的第5個記錄(基於1).fastq.gz虛擬:5-5 #此命令將檢索第5條記錄,並將其轉換回FASTQ格式。tabix output.fastq.gz dummy:5-5 | perl -pe's / ^ dummy \ t \ d + \ t //'| tr'\ t''\ n'#計算此問題第二部分的記錄數export N_RECORDS = $(gzip -dc output.fastq.gz | wc -l)
現在,我們有了一種在 log(n)
時間內檢索任意記錄的方法,檢索隨機記錄只是獲得一個良好的隨機數的問題發生器和採樣。以下是在python中執行此操作的一些示例代碼:
#random_read.pyimport osimport randomn_records = int(os .environ [“ N_LINES”])// 4rand_record_start = random.randrange(0,n_records)* 4 + 1rand_record_end = rand_record_start + 3os.system(“ grabix抓取output.fastq.gz {0} {1}”。format(rand_record_start ,rand_record_end))
#random_read.pyimport osimport randomn_records = int(os.environ [“ N_RECORDS”])rand_record_index = random.randrange(0,n_records)+ 1#非常難看,但是有效... os.system(“ tabix output.fastq.gz dummy:{0}-{0} | perl -pe's / ^ dummy \ t \ d + \ t //'| tr'\ t'' \ n'“。format(rand_record_index))
這對我有用:
python3.5 random_read.py
請注意, os.system
調用系統外殼程序,容易受到 shell注入漏洞。如果您正在編寫生產代碼,則可能需要採取額外的預防措施。
感謝 Chris_Rands提出了此問題。
正如 wkretzsch所建議的,這是一個實際的答案,我覺得這裡沒有明顯的解決方案。對 FASTQ
進行索引。
就像我通常不願意跳轉到需要腳本或框架的解決方案(相對於unix一樣)命令行工具),可悲的是沒有 samtools fqidx
(也許應該有),現有的答案顯示出很多麻煩。雖然它們可能起作用,但其中一些看起來很麻煩,並且有許多步驟可能會導致您犯錯。
因此,為了使事情簡單-一種快速而骯髒的替代方法可能是使用 biopython
,因為它已經已經實現了此功能來執行此操作,並且如果已安裝,則使用起來很簡單:
從生物導入SeqIOfq = SeqIO.index(“ myfastq.fq”,“ fastq”)
獲取索引後,您將獲得對以下任何內容的快速隨機訪問:
#通過讀取名稱的隨機訪問(我喜歡貓頭鷹)record = fq [“ HOOT”] record#> SeqRecord(seq = Seq('ACGTACGT ',SingleLetterAlphabet()),id ='HOOT',name ='HOOT',description ='HOOT',dbxrefs = [])#我們可以獲取sequencerecord.seq#> Seq('ACGTACGT',SingleLetterAlphabet()) #and qualityitiesrecord.letter_annotations#> {'phred_quality':[41,41,41,41,41,41,41,41,41]}
如果nt選擇任意隨機記錄,您可以使用 randrange
之類的東西在0和引用列表的長度之間進行選擇。
如果您想要多個記錄,則可能需要 sample
(以避免替換):
from random import sampleN = 100 random_indices = sample(xrange(len(key_list)),N)for random_indices中的key_i:random_readname = key_list [key_i] rand_record = fq.get(random_readname)#... 對於有價值的代碼,我認為 biopython
在RAM中保存此索引,因此,如果文件絕對龐大,則可能需要更聰明。如果是這種情況,您可以遍歷 FASTQ
一次,然後輸出讀取名稱,文件偏移量和長度-類似於 FASTA
fai
。 更新請參見 samtools fqidx
您可以洗一次FASTQ,然後根據需要從文件頂部讀取序列:
gzip -dc input.fastq .gz |粘貼----| shuf | tr'\ t''\ n'| gzip -c > output.fastq.gz
我建議在壓縮步驟中使用 pigz替代gzip。
這種方法的缺點是,在需要再次運行shuffle之前,您只會得到 n
個讀,並且顯然 shuf
擁有所有內容數據存儲在RAM中,因此如果FASTQ文件不符合問題中指定的大小,它將因內存不足錯誤而消失。
使用 sort -R
是複雜性 n log(n)
並使用臨時文件,因此它不應用完內存:
gzip -dc input.fastq.gz |粘貼----| nl |排序-R | perl -pe's / \ s * \ d + \ t //'| tr'\ t''\ n'| gzip -c > output.fastq.gz
必須使用 nl
和 perl
命令來確保未對相同的記錄進行排序彼此相鄰。
一種可能是:
這將為您獲取日誌(n)您想要的查找時間。
當然,數據不會被壓縮。您可以對基數進行2位編碼並量化質量,以節省一些空間,但這將是有損的,也許不是您想要的。或者,您可以對重新格式化的數據進行塊gzip壓縮,並記錄文件中有多少塊以及每個塊中有多少次讀取(因為文件大小將不再反映文件中的記錄數)。然後,要獲取特定的讀數,您需要計算將出現在其中的塊編號,解壓縮該塊,然後返回適當的讀數。
幾年前,賈里德·辛普森(Jared Simpson)在一篇博客文章中給出了對該問題(或類似問題:隨機讀取一部分)最徹底的處理方法之一。 http://simpsonlab.github.io/2015/05/19/io-performance/
如果您只想隨機讀取一次,則Jared的基準測試建議到文件中的隨機位置,然後檢索下一個完整的讀取應該是性能最高的選項。
但是,更一般而言,如果您希望隨機讀取一部分數據,則有許多因素會影響性能。
要強調此問題(如 1 / n
更新中所述),請考慮一個輸入文件,其中一個記錄的長度為500萬個鹼基,而一個記錄的長度為100個鹼基。您希望選擇兩個記錄中任何一個的概率相等。任何隨機查找方法都將壓倒性地挑選長記錄。
我希望索引記錄開始位置的確是這裡唯一可行的選擇,尤其是在可能進行多行記錄的情況下。創建一個包含每個記錄的起始位置的索引文件(大小相同的整數,例如64位),然後從索引文件中進行採樣(每個記錄的長度相同)以獲取起始位置。我設想此文件將僅包含開始位置。任何其他元數據(包括序列名稱)都需要在原始文件中進行搜索。
完成索引後,可以使用bgzip壓縮文件,並使用 -b
和 -s
選項。但是,我希望如果需要多個隨機記錄,壓縮將不會特別有效。
我編寫了一種名為 sample
的工具,您可以使用該工具進行隨機採樣,而無需將整個文件讀入內存。
它可以用於GNU shuf
因缺少足夠內存而失敗的情況。
它需要兩次通過文件才能進行隨機採樣,但第二次通過通常會更快,因為它使用 mmap
例程進行緩存的讀取。
如果您重複採樣,則重複採樣也會被 mmap
編輯(緩存),並且可以快速運行。
您可以在FASTQ上使用它像這樣的文件:
$ sample -k 1234 -l 4 in.fq > out.fq
它按每四個換行符(例如FASTQ文件的格式)將輸入文件解析為記錄,並將行偏移位置讀取到內存中。因此,內存開銷相對較低。
然後對這些行偏移量應用存儲庫採樣,以將隨機採樣(例如本示例中的 1234
記錄)寫到標準輸出中
這是另一種不需要任何索引的方法,使用 BEDOPS bedextract
對已排序的BED進行 log(n)
樣本文件。您的樣本包含概率為 1 / n
的隨機記錄。
此方法需要在文件中傳遞一個 O(n)
並將其轉換為一個BED文件:
$ cat records.fastq |粘貼----| awk'{print“ chrZ \ t” s“ \ t”(s + 1)“ $ 0}'> records.bed
將間隔存儲在單獨的文件中:
$ cut -f1-3 records.bed > interval.bed
要隨機抽取 k
元素,重新排列間隔文件並保留重新排列的元素的順序。
您可以使用我前面概述的 sample
工具 :
$ sample -k $ {K} -s interval.bed > interval-sample.bed
或者您可以 shuf
和 sort-bed
進行相同的操作:
$ shuf -n $ {K} interval.bed |排序床-> interval-sample.bed
這裡有一個 O(klog(k))
的成本,但如果 k <<< n
,即您正在使用全基因組規模輸入,則此成本在 lo上攤銷g(n)
搜索性能。
接下來,使用 bedextract
對記錄進行二進制搜索,然後進行線性化以返回FASTQ:
$ bedextract records.bed interval-sample.bed |切-f4 | tr'\ t''\ n'> sample.fq
對於Unix I / O流,可以通過一次完成:
$ sample -k $ {K} -s interval.bed | bedextract records.bed-|切-f4 | tr'\ t''\ n'> sample.fq
通過將排序順序烘焙到 records.bed
中,可以保證您可以執行二進制搜索,即 log(n)
。
注意:此外,通過線性化將FASTQ輸入到BED文件並查詢BED間隔,您有相等的概率選擇任何一個間隔(間隔== FQ記錄)。您可以繪製無偏樣本,而無需創建和存儲單獨的索引。
我編寫了一個名為 strandex的工具,該工具與整個FASTQ記錄相匹配,該記錄從解壓縮文件內的隨機偏移開始(因為非塊壓縮gzip文件不可搜索)。這是對高棉中 screed
FASTQ解析器的此評論的回應,起初只是概念的證明,但對我很有用和其他一些人。總體思路是:
CLI腳本和python庫可以通過 pip install strandex
,默認情況下,腳本使用從隨機數開始的文件偏移量(通過設置 -s
種子可複制)對腳本 -n
進行採樣,然後移動多個字節以統一採樣整個文件,具體取決於文件大小。這樣,過程並不是真正隨機的,而是足夠隨機的,並且避免了對文件的任何部分進行過多采樣。
如果指定 -n
讀取的文件少於文件總數,則 down 採樣,如果 -n
大於您獲得 up 採樣的讀取總數。
此方法相對於其他方法的優點是-這是一種單遍方法,並且由於它使用的是python C內部正則表達式引擎,因此速度非常快。
我相信其中一些答案會誤解原始問題。 OP(我相信)不希望找到特定的序列,因此不需要索引。我相信OP只是想要對數據集進行隨機抽樣。
我想我可以根據需要提供代碼,但是請記住,它可能必須用C語言完成。
可以使用bgzf壓縮的fastq文件來簡化此操作。是的,BAM已經使用了它,但是由於Illumina的軟件現在默認生產,因此應該已經有相當數量的可用。
那並不是嚴格意義上的O(log N),但與對整個文件進行解析並進行改組或以其他方式完全將數據整理成本來不是很有用的格式相比,它對內存和IO的要求較低。這還具有允許更快地選擇一個以上隨機條目的好處,因為您可以記憶事物。
一種快速而又骯髒的解決方案是將FASTQ文件轉換為兩個FASTA,其中一個存儲基礎,另一個存儲質量,然後使用標準方法隨機訪問FASTA。
唯一的麻煩是基本品質可以包含 >
。但是,只有當它是一行的第一個字符時,這才是問題,我們可以通過在每個質量序列前添加例如 I
來解決該問題。
如果還有其他人認為這可能是一個足夠的解決方案,我將擴大答案並添加所有命令。