我認為很難快速完成此工作-與這個問題一樣，處理大型壓縮的FASTQ文件大多受IO限制。相反，我們可以專注於確保得到正確的答案。

人們經常嘲笑他們，但這是一個編寫良好的解析器值得它在黃金中發揮作用的地方。亨利（Heng Li）向我們提供了此 FASTQ Parser in C。

我下載了示例tarball並修改了示例代碼（對不起，我的C ...）：

  #include <zlib.h>＃include <stdio.h>＃include“ kseq.h” KSEQ_INIT（gzFile，gzread）int main（int argc，char * argv []）{gzFile fp; kseq_t * seq;國際如果（argc == 1）{fprintf（stderr，“用法：％s <in.seq> \ n”，argv [0]）；返回1; } fp = gzopen（argv [1]，“ r”）; seq = kseq_init（fp）; int seqcount = 0;長序列= 0;而（（l = kseq_read（seq））> = 0）{seqcount = seqcount +1; seqlen = seqlen +（長）strlen（seq->seq.s）; } kseq_destroy（seq）; gzclose（fp）; printf（“序列數：％d \ n”，seqcount）； printf（“序列中的鹼基數：％ld \ n”，seqlen）；返回0;}  

然後 make 和 kseq_test foo.fastq.gz 。

對於我的示例文件（〜35m讀到〜75bp）花費了：

 真實0m49.670suser 0m49.364ssys 0m0.304s  

與您的示例相比：

 真實0m43.616suser 1m35.060ssys 0m5.240s  

Konrad的解決方案（在我手中）：

 真實0m39.682suser 1m11.900ssys 0m5.112s  

（順便說一下，只需將數據文件zcat到/ dev / null）：

 真實0m38.736suser 0m38.356ssys 0m0.308s  

因此，我的速度非常接近，但可能會更符合標準。同樣，該解決方案為您提供了更多數據處理靈活性。

幾乎可以肯定地優化了我的可怕C。

使用Github的 kseq.h 進行相同的測試，如註釋中所建議：

今天早上我的機器承受的負載不同，所以我已經重新測試。掛鐘時間：

OP：0m44.813s

Konrad：0m40.061s

zcat > / dev / null ：0m34 .508s

kseq.h（Github）：0m32.909s

因此，最新版本的 kseq.h 比簡單地對文件進行zcat處理要快（始終在我的測試中...）。

以下內容的速度快兩倍以上；但是， wc 也會計算換行符。因此，我們需要從基本計數中減去行計數：

  fix_base_count（）{本地計數=（$（cat））echo“ $ {counts [0]} $（（$ { counts [1]}-$ {counts [0]}））“} gungun -c” $ file“ awk'NR％4 == 2' wc -cl \ | fix_base_count  

但是，該字符也會計算換行符。因此，我們需要從中減去行數：

Simon註釋中的所有警告都適用：這假設“簡單” FASTQ格式，其中每條記錄正好由四行組成。我認為這對於Illumina定序器和下游工俱生成的所有文件都是正確的。

pigz | awk | wc是最快的方法

首先要使用FASTQ進行基準測試，最好使用具有已知答案的特定實際示例。我選擇了以下文件：

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz

作為我的測試文件，正確的答案是：

 讀取次數：67051220讀取的鹼基數：6034609800  

接下來，我們想找到最快的計數方法，所有時間都是在其他未加載的系統上，用bash time 收集的10次運行的平均掛鐘時間（實際）。 / p>

zgrep

  zgrep。 ERR047740_1.filt.fastq.gz | awk'NR％4 == 2 {c ++; l + = length（$ 0）} END {print“讀取次數：” c; print“讀取的鹼基數：” l}' 

這是最慢的方法，平均運行時間為125.35秒

gzip awk

使用 gzip ，我們又獲得了大約10秒鐘的時間：

  gzip -dc ERR047740_1.filt.fastq.gz | awk'NR％4 == 2 {c ++; l + = length（$ 0）} END {print“讀取次數：” c; print“讀取的鹼基數：” l}' 

平均運行時間為116.69秒

Konrad的gzip awk wc變體

  fix_base_count（）{本地計數=（$（cat））echo“ $ {counts [0]} $（（$$ counts [1]}-$ {counts [0]}））”} gzip -dc ERR047740_1 .filt.fastq.gz \ | awk'NR％4 == 2' wc -cl \ | fix_base_count  

在此測試文件上運行的速度比解決方案的 gzip awk 變體慢，平均運行時間為122.28秒。

kseq_test使用來自 klib

的最新 kseq.h 進行以下編譯的代碼： gcc -O2 -o kseq_test kseq_test.c -lz 其中 kseq_test.c 是Simon對Heng Li的FASTQ解析器的改編。

kseq_test ERR047740_1.filt.fastq.gz

平均運行時間為99.14秒，比到目前為止基於 gzip 核心實用程序的解決方案要好，但是我們可以做得更好！

piz awk

使用Mark Adler的 pigz代替 gzip ，請注意， pigz 在 gzip 之上為我們提供了速度上的提高，因為它除了主要的deflate線程外，還使用另外3個線程進行讀取，寫入和校驗和計算，請參見手冊頁有關詳細信息。

  pigz -dc ERR047740_1.filt.fastq.gz | awk'NR％4 == 2 {c ++; l + = length（$ 0）} END {print“讀取次數：” c; print“讀取的鹼基數：” l}' 

現在平均運行時間為93.86秒，比基於kseq的C代碼快約5秒，但我們可以

pigz awk wc

接下來，我們將 pigz 替換為Konrad的 wc awk 基於解決方案的code>變體。

  fix_base_count（）{本地計數=（$（cat））echo“ $ {counts [0]} $（ （$ {counts [1]}-$ {counts [0]}））“} gzip -dc ERR047740_1.filt.fastq.gz \ | awk'NR％4 == 2' wc -cl \ | fix_base_count  

現在平均運行時間降至83.03秒，這比基於kseq的解決方案快約16秒，比OPs zgrep 為基礎的解決方案。

接下來作為基準，讓我們看看運行時間中有多少是由於對輸入的 fastq.gz 文件進行解壓縮而導致的。

gzip單獨

gzip -dc ERR047740_1.filt.fastq.gz > / dev / null

平均運行時間： 105.95秒，因此基於 gzip 的解決方案（其中還包括 zcat 和 zgrep ，因為這些由 gzip 提供）永遠不會比 kseq_test 快。

pigz

pigz -dc ERR047740_1.filt.fastq.gz > / dev / null

平均運行時間：77.66秒，因此很明顯，另外三個線程用於讀取，寫入和校驗和計算提供了有用的優勢。更重要的是，利用 awk |時，此加速效果更大。基於wc 的解決方案，原因尚不清楚，但是我希望這是由於額外的寫線程所致。

有趣的是，所有線程的平均CPU使用率對於各種答案都非常有用，我已經使用GNU time / usr / bin / time --verbose p整理了這些統計信息基於>

zgrep 的解決方案133％-必須以某種方式超過一個線程

gzip |基於awk 的解決方案99％-所有基於 gzip 的解決方案都以99％的CPU使用率運行單線程

pigz | awk 147％

gzip | awk | wc 99％與 gzip

pgiz | awk | wc 155％

kseq_test 99％

gzip > dev / null 99％

pigz > dev / null 155％

雖然 pigz 中的主deflate線程將在100％CPU負載下運行，但額外的3個不會相當多的內核佔用了100％的內核（通過平均150％的CPU使用率證明），但是它們確實可以減少運行時間。

我正在使用Ubuntu 16.04.2 LTS ** ，我的 gzip ， zcat ， zgrep 版本均為gzip 1.6，而 pigz 版本為2.3.1。 gcc 是5.4.0版

**我認為我的補丁程序級別實際上是16.04.4，但是我沒有重啟170天：p

我的 fastx-length.pl腳本獲得了相當快的結果，並且具有能夠處理多行FASTQ文件並顯示其他讀取長度QC統計信息的額外好處：

 時間zcat albacored_all.fastq.gz | /bioinf/scripts/fastx-length.pl > / dev / null總序列：301135總長度：283.902419 Mb最長序列：5.601 kb最短序列：6 b平均長度：942 b中位長度：999 bN50：111835序列； L50：1.103 kb N90：245243序列; L90：608 breal 0m8,802suser 0m16,584ssys 0m0,260s  

相對於您提供的腳本：

  zcat albacored_all.fastq.gz | awk'NR％4 == 2 {c ++; l + = length（$ 0）} END {print“讀取次數：” c; print“讀取的鹼基數：” l}'讀取的鹼基數：301135讀取的鹼基數：283902419real 0m8,382suser 0m10,216ssys 0m0,332s  

Cat to / dev / null 進行比較：

 時間zcat albacored_all.fastq.gz > / dev / nullreal 0m7,877suser 0m7,856ssys 0m0,020s  

我懷疑使用 bioawk的程序可能會更快（並且類似FASTQ兼容）。

我已經使用 klib

中的kseq.h實現了 seqtk_counts只需幾行代碼：

  #include <stdio.h>＃include <stlib.h>＃include“ kseq.h” KSEQ_INIT（gzFile，gzread）int main（int argc，char * argv []）{gzFile fp; kseq_t * seq; int l = 0總數= 0； int64_t行= 0；如果（argc == 1）{fprintf（stderr，“用法：％s <fastq> <sample> \ n”，argv [0]）； return 1;} fp = strcmp（argv [1]，“-”）？ gzopen（argv [1]，“ r”）：gzdopen（fileno（stdin），“ r”）; seq = kseq_init（fp）;而（（l = kseq_read（seq））> = 0）{總計+ = seq ->seq.l;行+ = 1;} printf（“％s \ t％lld \ t％lld \ n”，argv [2]，（long long）行，（long long）總數）; ;返回0;  

}

編譯它：

  gcc -O2 seqtk_counts.c -o seqtk_counts -Iklib -lz  

用法：

  seqtk_counts foo.fasq.gz foo orcat foo.fasq.gz | seqtk_counts-foo  

使用pyGATB

（我使用與 https://bioinformatics.stackexchange.com/a/400/292中相同的文件，與中相同的工作站https://bioinformatics.stackexchange.com/a/380/292）

  $ time python3 -c“ from gatb import Bank; seq_lens = [len（seq）for seq in Bank（'SRR077487_2.filt.fastq.gz'）]; print（'讀取次數：％d'％len（seq_lens），'讀取的鹼基數量：％d'％sum（seq_lens），sep =' \ n'）“讀取次數：23861612讀取的鹼基數：2386161200real 0m41.122suser 0m40.788ssys 0m0.312s  

它比bioawk快得多：

  $ time bioawk -c fastx'{nb_seq + = 1; nb_char + = length（$ seq）} END {打印“讀取次數：” nb_seq“ \ n讀取的鹼基數：” nb_char}'SRR077487_2.filt.fastq.gz讀取的數目：23861612讀取的鹼基數：2386161200real 1m3.182suser 1m2.916ssys 0m0.268s  

但不超過OP示例：

  $ time zgrep。 SRR077487_2.filt.fastq.gz | awk'NR％4 == 2 {c ++; l + = length（$ 0）} END {print“讀取次數：” c; print“讀取的鹼基數：” l}'讀取的鹼基數：23861612讀取的鹼基數：2386161200real 0m47.127suser 1m36.292ssys 0m6.796s  

或 the wc 的解決方案：

  $ fix_base_count（）{>本地計數=（$（cat））>回顯“ $ {counts [0]} $ （（（$ {counts [1]}-$ {counts [0]}）））“ >} $ time gunzip -c SRR077487_2.filt.fastq.gz | awk'NR％4 == 2'| wc -cl | fix_base_count23861612 2386161200real 0m44.915suser 1m12.000ssys 0m6.972s  

我沒有與基於C的解決方案進行比較。

zcat / dev / null 引用的代碼>如下：

  $ time zcat SRR077487_2.filt.fastq.gz > / dev / nullreal 0m39。 745suser 0m39.464ssys 0m0.252s  

我仍然對pyGATB的速度印象深刻

如果數據在SRA中，則存在sra-stat實用程序，該實用程序返回讀取，鹼基和質量分佈。這些都存儲在SRA文件中。使用--quick獲取存儲的統計信息或使用--statistics計算每個readgroup / barcode細分的其他值。sra-stat --quick SRR077487