題:
隨機訪問FASTQ文件
winni2k
2017-06-19 18:15:23 UTC
view on stackexchange narkive permalink

我想從 log(n)時間複雜度(大O標記)或更少。一條記錄被定義為相當於FASTQ格式的四行。記錄不適合RAM,需要存儲在磁盤上。理想情況下,我想以壓縮格式存儲讀段。

我更喜歡不需要任何額外文件(例如參考基因組)的解決方案。

標題這個問題僅提及FASTQ,因為FASTQ是用於在磁盤上存儲未對齊讀取的通用格式。我對要求按時間複雜度順序 n 將數據進行一次有限的轉換為另一種文件格式的答案感到滿意。

Update

一個說明:希望以概率 1 / n 選擇隨機記錄。

相關問題:https://bioinformatics.stackexchange.com/q/453/292
十二 答案:
winni2k
2017-06-19 19:06:33 UTC
view on stackexchange narkive permalink

在恆定時間內訪問任意記錄

要在恆定時間內獲取隨機記錄,在恆定時間內獲取任意記錄就足夠了。

我在這裡有兩種解決方案:一種帶有 tabix ,另一種帶有 grabix 。我認為 grabix 解決方案更優雅,但我保留下面的 tabix 解決方案,因為 tabix 是比 grabix更成熟的工具

感謝 user172818建議 grabix

更新

此答案先前指出, tabix grabix log(n)時間內執行查找。在仔細研究了grapix的源代碼和tabix文件之後,我現在確信查找的複雜性與 n 無關。但是,這兩個工具都使用一個索引,該索引的大小與 n 成比例。因此,索引的加載順序為 n 。但是,如果我們認為索引的加載是“ ...將數據有限地轉換為另一種文件格式...”,那麼我認為這一答案仍然是有效的。如果要檢索多個記錄,則需要將索引存儲在內存中,也許使用諸如pysam或htslib之類的框架。

使用 grabix

  1. 使用 bgzip 壓縮。
  2. 為文件建立索引並使用 grabix
  3. 執行查找 ol>

    在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的文件格式感到滿意,則可以執行以下操作:

    1. 將每條記錄粘貼到一行中。
    2. 添加一個虛擬的染色體和行號作為第一列和第二列。
    3. bgzip 壓縮。
    4. 為文件建立索引並使用 tabix
    5. ol>

      請注意,我們需要刪除由 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) 

      Random在恆定時間內進行記錄

      現在,我們有了一種在 log(n)時間內檢索任意記錄的方法,檢索隨機記錄只是獲得一個良好的隨機數的問題發生器和採樣。以下是在python中執行此操作的一些示例代碼:

      使用grapix

       #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)) 

      使用tabix

       #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提出了此問題。

Sam Nicholls
2017-06-21 17:36:08 UTC
view on stackexchange narkive permalink

正如 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

Bio.SeqIO.index確實將索引作為Python字典保存在RAM中,但Bio.SeqIO.index_db將該索引放入可重用的SQLite3數據庫文件中。
@peterjc阿哈!很高興知道。
winni2k
2017-06-19 18:15:23 UTC
view on stackexchange narkive permalink

您可以洗一次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 命令來確保未對相同的記錄進行排序彼此相鄰。

kvg
2017-06-19 18:44:42 UTC
view on stackexchange narkive permalink

一種可能是:

  1. 重新格式化數據,以使每條記錄都是一行,其中包含讀取的描述,基礎和質量得分
  2. 填充每條記錄最大化每個字段中的長度,以使文件中的每個記錄都是相同的字節數
  3. 現在可以根據文件大小/記錄大小
  4. 選擇總記錄數來計算在0到記錄總數之間的隨機記錄數
  5. 對重新格式化的文件進行二進制搜索,直到獲得讀取的結果
  6. ol>

    這將為您獲取日誌(n)您想要的查找時間。

    當然,數據不會被壓縮。您可以對基數進行2位編碼並量化質量,以節省一些空間,但這將是有損的,也許不是您想要的。或者,您可以對重新格式化的數據進行塊gzip壓縮,並記錄文件中有多少塊以及每個塊中有多少次讀取(因為文件大小將不再反映文件中的記錄數)。然後,要獲取特定的讀數,您需要計算將出現在其中的塊編號,解壓縮該塊,然後返回適當的讀數。

Daniel Standage
2017-06-20 22:20:25 UTC
view on stackexchange narkive permalink

幾年前,賈里德·辛普森(Jared Simpson)在一篇博客文章中給出了對該問題(或類似問題:隨機讀取一部分)最徹底的處理方法之一。 http://simpsonlab.github.io/2015/05/19/io-performance/

如果您只想隨機讀取一次,則Jared的基準測試建議到文件中的隨機位置,然後檢索下一個完整的讀取應該是性能最高的選項。

但是,更一般而言,如果您希望隨機讀取一部分數據,則有許多因素會影響性能。

緩存可以大大縮短查詢時間。如果他們在fseek之上測試基於mmap的方法,那將很有趣。特別是使用SSD硬件。
gringer
2017-06-19 18:57:45 UTC
view on stackexchange narkive permalink

要強調此問題(如 1 / n 更新中所述),請考慮一個輸入文件,其中一個記錄的長度為500萬個鹼基,而一個記錄的長度為100個鹼基。您希望選擇兩個記錄中任何一個的概率相等。任何隨機查找方法都將壓倒性地挑選長記錄。

我希望索引記錄開始位置的確是這裡唯一可行的選擇,尤其是在可能進行多行記錄的情況下。創建一個包含每個記錄的起始位置的索引文件(大小相同的整數,例如64位),然後從索引文件中進行採樣(每個記錄的長度相同)以獲取起始位置。我設想此文件將僅包含開始位置。任何其他元數據(包括序列名稱)都需要在原始文件中進行搜索。

完成索引後,可以使用bgzip壓縮文件,並使用 -b -s 選項。但是,我希望如果需要多個隨機記錄,壓縮將不會特別有效。

關於您的第三段:我可能丟失了一些內容,但是索引不是`n`指令嗎?我認為這可以作為問題的預處理步驟嗎?
n被定義為第一句中“未比對的測序讀段”的數目。我不確定如何更清楚。但是,既然您知道這是什麼意思,請立即提出建議。
我將問題中的“記錄”與“讀取”混為一談。我將進行更具體的編輯。我採用的是biopython的措辭,它把序列記錄稱為序列,質量列表和描述的組合。
“任何需要通讀所有行的內容都將至少需要n,因此,除此以外的任何事情都將具有大於n的複雜度。”-這不是大O表示法的複雜度,這就是我在這裡使用的。我也將澄清這一點。
我真的很喜歡帶有塊壓縮的隨機尋道方法。我的問題是,如果記錄的大小不完全相同,則不能保證對記錄進行統一採樣。
Alex Reynolds
2017-06-20 01:19:14 UTC
view on stackexchange narkive permalink

我編寫了一種名為 sample 的工具,您可以使用該工具進行隨機採樣,而無需將整個文件讀入內存。

它可以用於GNU shuf 因缺少足夠內存而失敗的情況。

它需要兩次通過文件才能進行隨機採樣,但第二次通過通常會更快,因為它使用 mmap 例程進行緩存的讀取。

如果您重複採樣,則重複採樣也會被 mmap 編輯(緩存),並且可以快速運行。

您可以在FASTQ上使用它像這樣的文件:

  $ sample -k 1234 -l 4 in.fq > out.fq  

它按每四個換行符(例如FASTQ文件的格式)將輸入文件解析為記錄,並將行偏移位置讀取到內存中。因此,內存開銷相對較低。

然後對這些行偏移量應用存儲庫採樣,以將隨機採樣(例如本示例中的 1234 記錄)寫到標準輸出中

您的工具“ sample”使用儲層採樣,如果我沒有記錯的話,則為n階。但是,該問題要求使用順序`log(n)`算法。這可能是一個XY問題,因為該工具只需要一次通過就可以檢索“ k”條記錄(在您的情況下為1234)。但是,在我的特殊情況下,我不提前知道k,這就是要求log(n)複雜度的原因。
您可以省略-k來重新整理整個文件,或者在運行時派生它時將其指定為shell變量(即-k $ {K}`)。
Alex Reynolds
2017-06-20 02:30:21 UTC
view on stackexchange narkive permalink

這是另一種不需要任何索引的方法,使用 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記錄)。您可以繪製無偏樣本,而無需創建和存儲單獨的索引。

Matt Shirley
2017-06-20 05:43:57 UTC
view on stackexchange narkive permalink

我編寫了一個名為 strandex的工具,該工具與整個FASTQ記錄相匹配,該記錄從解壓縮文件內的隨機偏移開始(因為非塊壓縮gzip文件不可搜索)。這是對高棉 screed FASTQ解析器的此評論的回應,起初只是概念的證明,但對我很有用和其他一些人。總體思路是:

  1. 尋找文件中的隨機偏移量
  2. 讀取文件的一部分
  3. 測試正則表達式模式 @。+ [\ n \ r] +。+ [\ n \ r] + \ +。*?[\ n \ r]。+ [\ n \ r] 匹配
  4. 如果是,請使用正則表達式捕獲組提取匹配的FASTQ記錄
  5. 如果不匹配,請轉到步驟2
  6. ol>

    CLI腳本和python庫可以通過 pip install strandex ,默認情況下,腳本使用從隨機數開始的文件偏移量(通過設置 -s 種子可複制)對腳本 -n 進行採樣,然後移動多個字節以統一採樣整個文件,具體取決於文件​​大小。這樣,過程並不是真正隨機的,而是足夠隨機的,並且避免了對文件的任何部分進行過多采樣。

    如果指定 -n 讀取的文件少於文件總數,則 down 採樣,如果 -n 大於您獲得 up 採樣的讀取總數。

    此方法相對於其他方法的優點是-這是一種單遍方法,並且由於它使用的是python C內部正則表達式引擎,因此速度非常快。

涼!感謝您的回答。此工具是否可以保證對大小不均的記錄進行統一採樣?
我想您可以以保留記錄的概率'1 / L'對記錄執行拒絕採樣,其中'L'是記錄中的字符數...
不幸的是,@wkretzsch沒有。如果已將已修剪為不同長度的FASTQ文件配對,則該工具甚至會失敗。
Edward Kirton
2018-03-01 11:41:23 UTC
view on stackexchange narkive permalink

我相信其中一些答案會誤解原始問題。 OP(我相信)不希望找到特定的序列,因此不需要索引。我相信OP只是想要對數據集進行隨機抽樣。

  1. 如果所需的讀取次數很小,請確保可以在文件上查找並返回下一個完整的記錄。記住要讀取的讀取ID,無需替換。但是,在磁盤上跳轉的速度很慢,並且每次讀取的都是4k或任何大小的磁盤緩衝區。
  2. 如果所需的讀取次數很大,請讀取整個文件並決定
  3. 我的操作方式(針對Illumina):使用文件大小,讀取長度和標頭長度來估算讀取次數,確定我想要的讀取次數,跳過第一批讀取(由於邊緣效應,例如20-100k;您可以使用查找或僅從開始讀取並丟棄),並順序讀取下一批。我想這取決於您將讀取內容用作什麼;也就是說,如果您要從樣品中隨機讀取以進行生物學分析,或者您希望從運行中隨機讀取以進行QC分析。此方法更適合前者,但對後者也足夠好。
  4. ol>
Devon Ryan
2017-06-19 18:51:49 UTC
view on stackexchange narkive permalink

我想我可以根據需要提供代碼,但是請記住,它可能必須用C語言完成。

可以使用bgzf壓縮的fastq文件來簡化此操作。是的,BAM已經使用了它,但是由於Illumina的軟件現在默認生產,因此應該已經有相當數量的可用。

  1. 計算文件中的塊數。
  2. 從這些塊中隨機選擇一個。
  3. 對該塊中的條目執行儲層選擇。
  4. ol>

    那並不是嚴格意義上的O(log N),但與對整個文件進行解析並進行改組或以其他方式完全將數據整理成本來不是很有用的格式相比,它對內存和IO的要求較低。這還具有允許更快地選擇一個以上隨機條目的好處,因為您可以記憶事物。

Karel Brinda
2017-06-19 20:02:16 UTC
view on stackexchange narkive permalink

一種快速而又骯髒的解決方案是將FASTQ文件轉換為兩個FASTA,其中一個存儲基礎,另一個存儲質量,然後使用標準方法隨機訪問FASTA。

唯一的麻煩是基本品質可以包含 > 。但是,只有當它是一行的第一個字符時,這才是問題,我們可以通過在每個質量序列前添加例如 I 來解決該問題。

如果還有其他人認為這可能是一個足夠的解決方案,我將擴大答案並添加所有命令。

我認為,如果可以從此方案中重建整個隨機選擇的FASTQ記錄(包括描述),這將是一個可接受的解決方案?
是。每個記錄的長度都是相同的(假設“ I” hack包括一個對應的“-” hack序列),因此,如果知道文件A在序列記錄的開始位置在哪裡,那麼相同的位置將是開始位置文件B中的序列記錄。
@gringer,閱讀的說明如何?
@wkretzsch您有什麼特別的念頭嗎?如果是,您能為此創造一個主旨嗎?基本上我見過的所有FASTQ的+和@行都相同。
在這種成對的情況下,質量和順序的標頭將是相同的。
[此讀物](https://gist.github.com/wkretzsch/97fd6f182ea40a4dde937316fa026905)來自瓶裝Ashkenazi猶太三人組的基因組。它的@行包含多個字符,但+行僅包含一個字符。此外,對我來說,所有@行都具有相同長度的假設似乎[evil](http://www.cs.technion.ac.il/users/yechiel/c++-faq/defn-evil.html)。
標題行的長度不是相同的。每次讀取的標題行與其質量對應的行相同。無論如何,對於每個序列來說,這也不會以相同的概率處理更新的請求。


該問答將自動從英語翻譯而來。原始內容可在stackexchange上找到,我們感謝它分發的cc by-sa 3.0許可。
Loading...