質問:
FASTA / FASTQファイル内の未知のヌクレオチドの数を計算する最も速い方法は何ですか?
Kamil S Jaron
2017-06-01 21:47:03 UTC
view on stackexchange narkive permalink

以前は、基本的な統計情報が通常利用できる公開されているゲノム参照を使用していましたが、利用できない場合は、一度だけ計算する必要があるため、パフォーマンスについて心配する必要はありません。

最近、中規模のゲノム(〜Gbp)を持ついくつかの異なる種のシーケンスプロジェクトを開始し、異なるアセンブリパイプラインのテスト中に、生の読み取り(fastq)とアセンブリスキャフォールド(fasta)の両方で未知のヌクレオチドの数を何度も計算しました。したがって、計算を最適化したいと思いました。

  • 私にとっては、4行形式のfastqファイルを期待するのが妥当ですが、それでも一般的な解決策が好まれます
  • 解決策がgzipされたファイルでも機能するのであればいいでしょう

Q:fastaおよびfastqファイルの未知のヌクレオチド(N)の数を計算するための(パフォーマンス面での)最速の方法は何ですか? ?

おそらく、ヌクレオチドの総数とNの数(Nの比率を与えるため)が必要になる場合があります。これは、これに加えてfastaのいくつかの他のメトリックを提供するかなり高速な実装です:https://github.com/ENCODE-DCC/kentUtils/blob/master/src/utils/faSize/faSize.c
十 答え:
#1
+22
user172818
2017-06-02 02:46:23 UTC
view on stackexchange narkive permalink

FASTQの場合:

  seqtk fqchk in.fq | head -2  

ただし、正確な数ではなく、「N」塩基のパーセンテージが表示されます。

FASTAの場合:

  seqtk comp in.fa | awk '{x + = $ 9} END {print x}'  

このコマンドラインはFASTQでも機能しますが、awkが遅いため遅くなります。

編集:わかりました。@ BaCHのリマインダーに基づいて、ここに進みます(コンパイルするには kseq.hが必要です):

  / /コンパイルする:gcc -O2 -o count-N this-prog.c -lz#include <zlib.h>#include <stdio.h>#include <stdint.h>#include "kseq.h" KSEQ_INIT [256] = {0、1、2、3、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4 、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、5、4、4 、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、0、4、1、4、4、4、2、4 、4、4、4、4、4、4、4、4、4、4、4、4、3、4、4、4、4、4、4、4、4、4、4、4、4、4、0 、4、1、4、4、4、2、4、4、4、4、4、4、4、4、4、4、4、4、4、3、4、4、4、4、4、4、4 、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4 4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4 4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4 4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4 4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4、4 4、4、4、4、4、4、4、4、4、4、4、4}; int main(int argc、char * argv []){long i、n_n = 0、n_acgt = 0、n_gap = 0; gzFile fp; kseq_t * seq; if(argc == 1){fprintf(stderr、 "使用法:count-N <in.fa> \ n"); 1を返します。 } if((fp = gzopen(argv [1]、 "r"))== 0){fprintf(stderr、 "エラー:入力ファイルを開くことができません\ n"); 1を返します。 } seq = kseq_init(fp);
while(kseq_read(seq)> = 0){for(i = 0; i < seq->seq.l; ++ i){int c = dna5tbl [(unsigned char)seq->seq.s [i]]; if(c < 4)++ n_acgt; else if(c == 4)++ n_n;それ以外の場合++ n_gap; }} kseq_destroy(seq); gzclose(fp); printf( "%ld \ t%ld \ t%ld \ n"、n_acgt、n_n、n_gap); return 0;}  

FASTA / Qとgzip圧縮されたFASTA / Qの両方で機能します。以下はSeqAnを使用します:

  #include <seqan / seq_io.h>using namespace seqan; int main(int argc、char * argv []){if( argc == 1){std :: cerr << "使用法:count-N <in.fastq>" << std :: endl; 1を返します。 } std :: ios :: sync_with_stdio(false); CharString id; Dna5String seq; SeqFileIn seqFileIn(argv [1]); long i、n_n = 0、n_acgt = 0; while(!atEnd(seqFileIn)){readRecord(id、seq、seqFileIn); for(i = beginPosition(seq); i < endPosition(seq); ++ i)if(seq [i] < 4)++ n_acgt;それ以外の場合++ n_n; } std :: cout << n_acgt << '\ t'<< n_n << std :: endl; return 0;}  

400万回の150bp読み取りを伴うFASTQの場合:

  • Cバージョン:〜0.74秒
  • C ++バージョン:〜2.15秒
  • ルックアップテーブルのない古いCバージョン(前の編集を参照):〜2.65秒
Devon Ryanのベンチマークに基づくと、これまでで最速です...
最初のCの例に `#include `を追加できますか?
SeqAnに賛成することはできませんが、C ++コードがCコードよりもはるかに遅いことを深刻に懸念しています。それらは実質的に同一である必要があります。とはいえ、 `std :: ios :: sync_with_stdio(false)`を指定しなかったため、違いの一部が説明されている可能性があります。また、イテレータを使用することでループを簡略化できますが、これが「-O1」を超える設定の速度に影響を与えた場合は驚きます。
@KonradRudolph `sync_with_stdio`は、通常のファイルから読み取るときにほとんど役に立ちません。 seqanの内部がわからないので、なぜ遅いのか説明できません。そうは言っても、私はseqanのパフォーマンスについてあまり心配しません。 Gzipの解凍は、fastqの解析よりもはるかに低速です。 @DevonRyanが追加されました。
違いはreadfqにあります。その獣は*速い*。私は、HengのコードをC ++でそれほど醜いものに置き換えようとしましたが、単純な読み取りパフォーマンスではそれに近づくことができませんでした。クリーンなC ++で思いつくことができたのは、速度係数が2倍遅いことでした。
コンパイルフラグに `-lz`を追加する必要があります(少なくとも現在のUbuntu LTSの場合)` gcc -O2 -o count-N this-prog.c -lz`
あなたの数字はサンプル間のキャッシュのフラッシュに基づいていますか?そうでなければ、それらの数字はほぼ間違いなく正しくありません。
#2
+18
Devon Ryan
2017-06-02 03:20:21 UTC
view on stackexchange narkive permalink

5時間で、ベンチマークは投稿されていませんか?とてもがっかりしました。

fastqは同じになるので、比較をfastaファイルだけに制限します。これまでのところ、候補は次のとおりです。

  1. R ShortRead パッケージ(最速ではないにしても、確かに非常に便利な方法)。
  2. grep -v "^ >"のパイプライン| tr -cd A | wc -c
  3. grep -v "^ >"のパイプライン| grep -io A | wc -l
  4. seqtkcompのパイプライン| awk
  5. C / C ++のカスタム
  6. ol>

    Rコードを試してみると、何らかの理由で実際にはNレコードをカウントできないため、カウントしましたすべての A

    pythonが高速であると考えられる宇宙はないため、すべてのpythonソリューションを除外します。 C / C ++ソリューションでは、@ user172818が投稿したものを使用しました(作成しようとしたものはすべて同じ時間がかかりました)。

    100回繰り返して平均を取ります(はい、中央値を取る必要があります)。

    1. 1.7秒
    2. 0.65秒
    3. 15秒
    4. 1.2秒
    5. 0.48秒
    6. ol>

      当然のことながら、ストレートCまたはC ++のすべてが勝ちます。 tr を使用した grep は驚くほど優れています。これは、 grep が非常に最適化されているにもかかわらず、オーバーヘッドが増えると予想していたためです。 grep -io へのパイピングは、パフォーマンスを低下させます。通常のRオーバーヘッドを考えると、Rソリューションは驚くほど優れています。

      Update1:​​コメントで示唆されているように

    7. sum(x.count( 'A')for x in open( 'seqs.fa'、 'r')if x [0]!= '>') in python
    8. ol >

      これにより

    9. 1.6秒の時間が得られます(私は現在仕事をしているので、異なる時間を考慮してできる限り時間を調整しましたコンピューター)
    10. ol>

      Update2:コメントからのその他のベンチマーク

    11. awk -FA '! / ^ > "/ {cnt + = NF-1} END {print cnt} 'foo.fa
    12. perl -ne'if(!/ ^ > /){$ count + = tr / A //} END {print "$ count \ n"} 'foo.fa
    13. ol>

      これらの生成時間:

      1. 5.2秒
      2. 1.18秒
      3. ol>

        awkバージョンは私が思いつくことができる最高のハックです。おそらくもっと良い方法があります。

        アップデート3 :そうです、fastqファイルに関しては、これを実行しています。 3.3GBのgzipファイルを10回繰り返して(これは実行に少し時間がかかります)、最初は圧縮ファイルを処理するために簡単に変更できるコマンドにテストを制限するつもりはありません(結局のところ、非圧縮のfastqファイルは忌まわしいものです) 。

        1. seqtk fqchk
        2. sed -n '2〜4p' <(zcat foo.fastq.gz)| grep -io A | wc -l
        3. sed -n '1d; N; N; N; P; d'<(zcat Undetermined_S0_R1_001.fastq.gz)| tr -cd A | wc -c
        4. user1728181のCの例(C ++との比較はすでに提供されているため、含める必要はありません)。
        5. bioawk -c fastx '{n + = gsub(/ A /、 ""、$ seq)} END {print n}' foo.fastq.gz
        6. ol>

          平均時間は次のとおりです。

          1. 1分44秒
          2. 2分6秒
          3. 1分52秒
          4. 1分15秒
          5. 3分8秒
          6. ol>

            注意:これらは、4行を超えるエントリを持つfastqファイルを処理しないことがよくあります。ただし、そのようなファイルはバイオインフォマティクスの鳥であるため、実際には重要ではありません。また、 ShortRead を圧縮されたfastqファイルで動作させることができませんでした。それを修正する方法があると思います。

素晴らしいベンチマークです。テストfastaファイルのサイズはどれくらいですか?また、時間があればPythonを追加してみませんか。そうです、最速ではありませんが、比較のために見てみると面白いでしょう。 `sum(x.count( 'A')for x in open( 'seqs.fa' 、 'r')if x [0]!= '>') `、おそらく` BioPython`はこのために遅くなります
@Chris_Rands dm6(非圧縮)を使用したので、それほど大きくはありません。 Pythonは明らかに遅いので、気にしたくなかったと思いましたが、追加しても問題はないと思います。
@Chris_RandsPythonベンチマークで更新しました。それは私が推測したよりも優れているRと同様に機能します。
ありがとう! Pythonの `str.count`はCで実装されているので、それほど驚くことではありません。また、高速の` awk`のみの解決策があるはずだと感じていますが、それを理解することはできません。
Perl: `` `perl -ne'if(!/ ^> /){$ count + = tr / A //} END {print" $ count \ n "} 'dm6.fa```
さて、質問が更新されたので、これらを変更して、一般的なFASTQの場合に機能するようにすることができます(複数行のシーケンス+品質文字列、 '@'が品質行に表示される可能性があります)。
私のawk-fuは、実用的な例を思い付くには十分ですが、痛々しいほど遅くないものを思い付くには十分ではありません。
@gringer今日は生産的になります...
私の `bioawk -c fastx '{n + = gsub(/ A /、" "、$ seq)} END {print n}'`を試してみませんか?それがあなたのデータの他のソリューションとどのように比較されるかを知りたいです。
@bli:うん、今は他に5人走っている...
Rコードは、正しく記述すればNをカウントできます。 ;-)(私の悪い)
awkバージョンとperlバージョンでは、テストfastaファイルで同じ結果が得られません。awkバージョンでは、C。elegansゲノムにperlバージョンよりも1つ多い「A」が見つかります。シーケンスの1つが「A」で始まるか終わるかによると思います。
いくつかのメソッド(bioawkを含む)とのfastq比較を追加しました。
@bliawkコードには `-1`が含まれている必要があります。代わりに `bioawk`を使用する多くの理由の中で!
bioawkではzcatは必要ないと思います。速度は向上しますか?
@bliそれは必要ではなく、それは私の側の間違いです(私は実際にベンチマークでzcatを使用しませんでした)。それをキャッチしてくれてありがとう!
なぜ `N`ではなく` A`を数えているのですか?それとも私は何かが足りないのですか?
元のShortReadバリアントはNを処理できず、パフォーマンスの比較は関係なく同じであるためです。 「Rコードを試してみると、なんらかの理由で実際にはNレコードをカウントできないため、すべてAをカウントしました。」
はい、申し訳ありませんが、私が最初に答えを読んだとき、私はそれを逃しました。
あなたはこう言います:「100回繰り返して平均を取る(はい、中央値を取るべきです):」しかし、平均は中央値ではなく平均です-それであなたはどちらを使いましたか?中央値はスキューの影響を受けにくいと思いますが、この場合、競合を引き起こしている他のユーザーからの非常に高い瞬間CPUまたはIO負荷のあるシステムで実行している場合を除き、実行時間は正規分布しています(この場合は平均で問題ありません)。 。
中央値ではなく平均値であるため、注意が必要です。はい、実行時間は正規分布である必要がありますが、時折奇妙なIOスパイクなどが発生します。
これらすべてのテストの間にシステムキャッシュをフラッシュしましたか?
いや、実際には、ファイルはすべての実行でキャッシュにありました。
テスト間でキャッシュをフラッシュする必要があります。テスト結果から平均または中央値を引き出すかどうかに関係なく、そうでなければ(まあ、間違って、本当に)あなたの数は誤解を招くでしょう。
#3
+9
Karel Brinda
2017-06-01 22:22:20 UTC
view on stackexchange narkive permalink

これはかなり速いはずだと思います:

FASTA:

  grep -v "^ >" seqs.fa | tr -cd N | wc -c  

FASTQ:

  sed -n '1d; N; N; N; P; d'seqs.fq | tr -cd N | wc -c  

さまざまなアプローチを使用してBASHの文字をカウントする方法については、この SOに関する回答を参照してください。

残念ながら、FASTQの回答には、基本的にIakovと同じ問題があります。つまり、FASTQレコードの行数は可変です。
元の標準ではFASTQの行分割が認められていることを私は知っています。しかし、実際にそのようなファイルを見たことがありますか?
私はそうしていませんが、サンガーシーケンシング/ PacBioデータを使ったことがありません。そこで出会うかもしれないと思います。短い読み取りデータの場合、「レコード= 4行」の仮定は安全だと思います。
長時間読み取られたデータはfastqを使用しません。現在、PacBioはbam(はい、本当にbam)とONT fast5(h5アーカイブ)を使用しています。 4行のfastqを想定するのは非常に合理的だと思います。
#4
+8
Iakov Davydov
2017-06-01 22:00:26 UTC
view on stackexchange narkive permalink

FASTQ

指摘されたように、fastqは複雑になる可能性があります。ただし、レコードごとに4行ある単純なケースでは、bashで考えられる解決策の1つは次のとおりです。

  sed -n '2〜4p' seqs.fastq | grep -io N | wc -l  
  • sed -n '2〜4p' は4行ごとに出力します
  • grep -o N は、一致するシンボルごとに N の行を出力します
  • wc -l は行をカウントします

このpythonアプローチはより速く機能すると思います:

  cat seqs.fastq | python3 -c "import sys; print(sum(line.upper()。count( 'N')for line in sys.stdin))"  

FASTA

Coreutils:

  grep -v "^ >" seqs.fasta | grep -io N | wc -l  

Python:

  cat seqs.fasta | python3 -c "import sys; print(sum(line.upper()。count( 'N')for line in sys.stdin if not line.startswith( '>')))"  
これで、NではなくNを含む行を数えています。また、FASTQレコードは4行を超える場合があります。 :-(
-ogrepを指定した@KonradRudolphは個々のシンボルを出力します。fastqコードを修正する方法を考えます
`@`は高品質の文字列行を開始することもできるため、FASTQは注意が必要です。 http://chat.stackexchange.com/transcript/message/37809200#37809200以降の行を参照してください。参考までに、FASTQを解析し、識別子からメイト情報を削除するPerlスクリプトを次に示します。それはおそらく出発点として役立つ可能性があります。これは最小限のコードに非常に近いと思います:https://gist.github.com/klmr/81a2b86cd93c706d28611f40bdfe2702
@KonradRudolph,これが単純な場合にのみ機能することを示すために、回答を更新しました
#5
+6
bli
2017-06-02 14:28:46 UTC
view on stackexchange narkive permalink

bioawkの使用

bioawkを使用( C. elegans ゲノムの「A」を数える、「N」がないように見えるため)このファイル)、私のコンピューター:

   $ time bioawk -c fastx '{n + = gsub(/ A /、 ""、$  span> seq)} END {print n} 'genome.fa 32371810real 0m1.645suser 0m1.548ssys 0m0.088s  

bioawkは、便利な解析オプションを備えたawkの拡張機能です。たとえば、 -c fastx は、シーケンスをfastaおよびfastq形式( gzip圧縮バージョンを含む)で $ seq として使用できるようにするため、上記のコマンドは fastq形式も堅牢に処理する必要があります

gsub コマンドは通常のawk関数です。置換文字の数を返します( https://unix.stackexchange.com/a/169575/55127から取得)。 awk内をハックの少ない方法でカウントする方法についてのコメントを歓迎します。

このfastaの例では、 grep tr codeのほぼ3倍の速度です。 >、 wc パイプライン:

  $ time grep -v "^ >" Genome.fa | tr -cd "A" | wc -c32371810real 0m0.609suser 0m0.744ssys 0m0.184s  

(タイミングを繰り返しましたが、上記の値は代表的なもののようです)

pythonの使用

readfq

この回答に記載されているリポジトリから取得したPythonreadfqを試しました: https://bioinformatics.stackexchange.com/a/365/292

  $ time python -c "from sys import stdin; from readfq import readfq; print sum((seq.count( 'A')for _、seq、_ in readfq(stdin)))" <ゲノム.fa32371810real 0m0.976suser 0m0.876ssys 0m0.100s  

"pure python"

この回答から採用されたものとの比較: https:// bioinformatics。 stackexchange.com/a/363/292

  $ time python -c "from sys import stdin; print sum((line.count( 'A')for line in stdin line.startswith( '>'))) "<ゲノムではない場合.fa32371810real0m1.143suser 0m1.100s
sys 0m0.040s  

(コンピューター上のpython2.7よりもpython3.6の方が遅い)

pyGATB

学習したばかりです(08 / 06/2017) GATBにはfasta / fastqパーサーが含まれており、最近PythonAPIがリリースされました。私は昨日それを使って現在の質問に対する別の答えをテストしようとしましたが、バグを見つけました。このバグは修正されたので、 pyGATBベースの回答を次に示します。

  $ time python3 -c "from gatb import Bank; print(sum((seq。 sequence.count(b \ "A \")for seq in Bank(\ "genome.fa \")))) "32371810real 0m0.663suser 0m0.568ssys 0m0.092s  

sequence.decode( "utf-8")。count( "A")を実行することもできますが、これは少し遅いようです。)

ここではpython3.6を使用しましたが(pyGATBはpython3のみのようです)、これは他の2つのpythonアプローチ(報告されたタイミングはpython 2.7で取得されます)よりも高速です。これは、 grep tr wc パイプラインとほぼ同じ速度です。

Biopython

さらに多くの比較を行うために、Biopythonの SeqIO.parse を使用したソリューションを次に示します(python2.7を使用):

  $ time python -c "from Bio import SeqIO; print(sum((rec.seq.count(\ "A \")for rec in SeqIO.parse(\ "genome.fa \"、\ "fasta \")))) "32371810real 0m1.632suser 0m1.532ssys 0m0.096s  

これは「純粋なPython」ソリューションよりも少し遅いですが、おそらくより堅牢です。

でわずかな改善があるようです下位レベルの SimpleFastaParser を使用するという@peterjcの提案:

  $ time python -c "from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count(タイトルは「A」)、SimpleFastaParser(open( 'genome.fa')))) "32371810real 0m1.618suser 0m1.500ssys 0m0.116s  

(一連のタイミングと代表的なものをとろうとしましたが、重複したwiがたくさんあります上位レベルのパーサーのタイミング。)

scikit-bio

今日(2017年6月23日)、シーケンスリーダーを備えた scikit-bioについて読みました。

  $ time python3 -c "import skbio; print(sum(seq.count( 'A')for seq in skbio.io.read( 'genome.fa'、format = 'fasta'、verify = False))) "32371810real 0m3.643suser 0m3.440ssys 0m1.228s  

pyfastx

今日(2020年8月27日)、シーケンスリーダーを備えた pyfastxについて読みました。

  $ time python3 -c "from pyfastx import Fasta; print(sum(seq.count( 'A')for(_、seq)in Fasta( 'genome.fa'、build_index = False)))" 32371810real 0m0.781suser 0m0.740ssys 0m0.040s  

(タイミングは、3年前の以前のタイミングと同じマシンでpython 3.6を使用して行われるため、比較可能である必要があります。)

fastx_read from mappy

私は(27/08/2020)から fastx_read のベンチマークを追加していなかったことに気づきました。 mappyは、かなり前からよく使用しています。

  $ time python3 -c "from mappy import fastx_read; print (sum(seq.count( 'A')for(_、seq、_)in fastx_read( 'genome.fa'))) "32371810real 0m0.731suser 0m0.648ssys 0m0.080s  

fastq.gzファイルのベンチマーク

ここで使用したのと同じファイルの「N」を数えて、上記のソリューション(またはその適応)のいくつかをテストしました: https:// bioinformatics .stackexchange.com / a / 400/292

bioawk

   $ time bioawk -c fastx '{ n + = gsub(/ N /、 ""、$  span> seq)} END {print n} 'SRR077487_2.filt.fastq.gz306072real 1m9.686suser 1m9.376ssys 0m0.304s  

pigz + readfqpythonモジュール

readfqは文句を言わず、圧縮されたfastqを直接渡すと非常に高速ですが、何か問題が返されるので、手動で解凍を行うことを忘れないでください。

ここで pigz を試してみました:

  $ time pigz -dc SRR077487_2.filt.fastq.gz | python -c "from sys import stdin; from readfq import readfq; print sum((seq.count( 'N')for _、seq、_ in readfq(stdin)))" 306072real 0m52.347suser 1m40.716ssys 0m8.604s  

コンピューターには16コアがありますが、 pigz の制限要因はディスクからの読み取りであると思われます。プロセッサーはフルスピードで実行するにはほど遠いです。

そして gzip の場合:

  $ time gzip -dc SRR077487_2.filt.fastq.gz | python -c "from sys import stdin; from readfq import readfq; print sum((seq.count( 'N')for _、seq、_ in readfq(stdin)))" 306072real 0m49.448suser 1m31.984ssys 0m2.312s  

ここでは、gzipとpythonの両方がフルプロセッサを使用していました。リソース使用量のバランスはわずかに改善されました。私はデスクトップコンピューターで作業しています。コンピューティングサーバーは pigz をより有効に活用すると思います。

pyGATB

  $ time python3 -c "from gatb import Bank; print( sum((seq.sequence.count(b \ "N \")for seq in Bank(\ "SRR077487_2.filt.fastq.gz \")))) "306072real 0m46.784suser 0m46.404ssys 0m0.296s  

biopython

  $ time python -c "from gzip import open as gzopen; from Bio.SeqIO.QualityIO import FastqGeneralIterator; print(sum(seq.count( ' N ')for(_、seq、_)in FastqGeneralIterator(gzopen(' SRR077487_2.filt.fastq.gz ')))) "306072real 3m18.103suser 3m17.676ssys 0m0.428s  

pyfastx(2020年8月27日追加)

  $ time python3 -c "from pyfastx import Fastq; print(sum(seq.count( 'N')for(_、seq、_ )in Fastq( 'SRR077487_2.filt.fastq.gz'、build_index = False))) "306072real 0m51.140suser 0m50.784ssys 0m0.352s  

fastx_read mappyから(2020年8月27日追加)

  $ time python3 -c "from mappy import fastx_read; pri nt(sum(seq.count( 'N')for(_、seq、_)in fastx_read( 'SRR077487_2.filt.fastq.gz'))) "306072real 0m52.336suser 0m52.024ssys 0m0.300s  

結論

pyGATBとbioawkはどちらも、透過的に圧縮(gzip圧縮されているかどうか)と形式の違い(fastaまたはfastq)を処理します。 pyGATBはまったく新しいツールですが、私がテストした他のpythonモジュールと比較して効率的です。

(更新27/08/2020)pyfastxとmappyも非常に便利で、pyGATBよりわずかに遅いだけです。 。

Biopythonの低レベル文字列FASTAパーサー番号を追加できますか? `` time python -c "from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count( 'A')for title、seq in SimpleFastaParser(open( 'genome.fa'))))" ``
質問のとおり、すべての例でN(およびおそらくn)ではなく文字Aがカウントされるのはなぜですか?
テスト用に考えた最初の大きなファイルにたまたま「N」がなかったので、文字「A」を数えます。
Pythonのロードとライブラリのインポートのオーバーヘッドのため、この小さなサンプルファイルではすべてのPythonの例のスコアが低くなると思います。しかし、ユーザーが1つのゲノムだけを数えたい場合、それはすべて同じように公正なテストです。
#6
+6
Matt Bashton
2017-06-03 01:28:27 UTC
view on stackexchange narkive permalink

gzip圧縮されたFASTQの解凍が主な問題です

些細なFASTAではなく、実際のgzip圧縮されたFASTQ(OPが提案したように有益です)を開始点として使用する場合、実際の問題は実際に解凍されます。 Nをカウントしないファイル、この場合はCプログラム count-N は、最速のソリューションではなくなりました。

さらに、ベンチマークに特定のファイルを使用するとよいでしょう。実際にはNがあります。これは、いくつかのメソッドでは、Nではなくより頻繁に発生するAsをカウントすることで、非常に興味深い実行時間の違いが得られるためです。

そのようなファイルの1つは次のとおりです。

http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/SRR077487_2.filt.fastq.gz

さまざまなものも確認する価値がありますソリューションは正解を返します。上記のファイルには306072Nが含まれているはずです。

次に、 / dev / null にリダイレクトされたこのファイルの解凍は zcat および gzip (どちらも私のシステムではgzip 1.6です) MarkAdlerの pigz のようなgzipの並列実装は、解凍に4つのスレッドを使用しているように見えます。すべてのタイミングは、実際の(実時間)を報告する平均10回の実行を表します。

  time pigz -dc SRR077487_2.filt.fastq.gz > / dev / nullreal 0m29.0132stime gzip -dcSRR077487_2。 filt.fastq.gz > / dev / nullreal 0m40.6996s  

2つの間に約11.7秒の違いがあります。次に、FASTQで実行され、正しい答えが得られる1ライナーのベンチマークを試みた場合(4行ではないFASTQにまだ遭遇しておらず、誰がこれらのファイルを生成するのかを真剣に考えています!)

  time pigz -dc SRR077487_2.filt.fastq.gz | awk'NR%4 == 2 {print $ 1} '| tr -cd N | wc -c306072real 0m34.793s  

ご覧のとおり、 pigz ベースの解凍と比較して、カウントにより合計実行時間が最大5.8秒長くなります。さらに、この時間の差は、 gzip をgzip解凍のみより約6.7秒上で使用すると、より高くなります。

  time gzip -dc SRR077487_2.filt.fastq.gz | awk'NR%4 == 2 {print $ 1} '| tr -cd N | wc -c 306072real 0m44.399s  

pigz awk tr wc ベースのソリューションは、 count-N よりも約4.5秒高速です。ベースのCコードソリューション:

  time count-N SRR077487_2.filt.fastq.gz 2385855128 306072 0real 0m39.266s  

この違いは、何度でも再実行します。 Cベースのソリューションでpthreadを使用するか、pthreadを変更して、pigzから標準を削除できるとしたら、パフォーマンスも向上すると思います。

別の代替 pigz grep のベンチマークバリアントは、 tr ベースのバリアントとほぼ同じ時間がかかるようです。

  time pigz -dc SRR077487_2.filt.fastq.gz | awk'NR%4 == 2 {print $ 1} '| grep -o N | wc -l306072real 0m34.869s  

上記の seqtk ベースのソリューションは著しく遅いことに注意してください。

  time seqtk comp SRR077487_2 .filt.fastq.gz | awk '{x + = $ 9} END {print x}' 306072real 1m42.062s  

ただし、 seqtk comp は他のソリューションよりも少し多くのことを行っていることに注意してください。 。

#7
+5
Konrad Rudolph
2017-06-01 22:10:03 UTC
view on stackexchange narkive permalink

正直なところ、(特にFASTQの場合)最も簡単な方法は、おそらくR / Bioconductorなどの専用の解析ライブラリを使用することです。

  suppressMessages(library(ShortRead))seq = readFasta(commandArgs (TRUE)[1])#またはreadFastqcat(colSums(alphabetFrequency(sread(seq))[、 'N'、drop = FALSE])、 '\ n') 

これは可能性があります最速ではありませんが、関連するメソッドが高度に最適化されているため、かなり高速です。最大のオーバーヘッドは、実際にはおそらく ShortRead のロードによるものです。これは不当なスローグです。

#8
+5
BaCh
2017-06-01 22:15:23 UTC
view on stackexchange narkive permalink

それがあなたが求めている生の速度であるなら、あなた自身の小さなC / C ++プログラムを書くことはおそらくあなたがする必要があることです。

幸いなことに、最悪の部分(高速で信頼できるパーサー)はすでに取り組んでいます:HengLiの readfqはおそらく最速のFASTA / FASTQパーサーです。

そして使いやすく、GitHubの例は簡単に拡張して何をすることができます。あなたが必要です。 (分析するファイル名の)パラメーターパーサーを追加し、単純な「N」カウンターループをプログラムして、結果をstdoutに出力するだけです。完了しました。

それ以外の場合、必要なのが単純な場合は、Rベースのソリューションに関するKonradの回答を参照してください。または、Biopythonを使用した非常に単純なスクリプト。

絶対に凶悪なコードである神。 ☹️生の速度については、[SeqAn](http://docs.seqan.de/seqan/master/?p=SeqFileIn)(C ++)はおそらく同等であり、はるかに高品質だと思います。
それは周りで最もきれいではありません。繰り返しになりますが、必要なことを実行し、ライブラリ全体の依存関係はありません。 SeqAnとreadfqのベンチマークはありますか?
readfq ansSeqAnがGATBのfasta / fastqパーサーとどのように比較されるのか知りたいです:http://gatb-core.gforge.inria.fr/doc/api/classgatb_1_1core_1_1bank_1_1impl_1_1BankFasta.html#details実際にコーディングすることはできませんただし、C / C ++では。
#9
+3
Matt Shirley
2017-06-15 20:47:20 UTC
view on stackexchange narkive permalink

FASTAファイルの場合、 pyfaidx v0.4.9.1で比較的効率的な方法を実装しました。この投稿により、以前のコードは非常に遅く、簡単に置き換えることができました。

  $ pip install pyfaidx $ time faidx -inucleotide〜 / Downloads / hg38.faname start end ATCGN otherschr1 1 248956422 67070277 67244164 48055043 48111528 18475410 chr10 1 133797422 38875926 39027555 27639505 27719976 534460 chr11 1 135086622 39286730 39361954 27903257 27981801 552880 chr11_KI270721v1_random 1 100316 18375 19887 31042 31012 0 ... chrX 1 156040895 37240 8274 12978 7232 87560実数2m28.251suser2m15.548ssys 0m9.951s  
#10
  0
Edward Kirton
2018-03-01 11:55:22 UTC
view on stackexchange narkive permalink

以下は、合計が必要な場合に十分な速度のfastaファイルとfastqファイルのperl1ライナーです。

$ perl -ne'BEGIN {my $ n = 0}; / ^ > /または$ n + = s / N // gi; END {print "N = $ n \ n"} 'contigs.fasta

FASTAファイルの各行について、ヘッダー(/ ^> /)でない場合は、シーケンス行です。置換Ns。これは、行われた置換の数を返し、合計に追加されます。

$ zcat read.fastq | perl -ne'BEGIN {my $ n = 0}; / ^ [ATCGN] + $ / iおよび$ n + = s / N // gi; END {print "N = $ n \ n"} '

/ ^ [ATCGN] + $ / i(つまり、すべてのNT)と一致するFASTQファイルのシーケンス行ごとに、次に、同じ方法で行われた置換をカウントします。

pigzがインストールされている場合は、zcat(またはgunzip -c)の代わりにpigzを使用できます。

14M 250bpIlluminaの非圧縮fastqファイル読み取りは40秒で完了しました。

各perl行が何をするのかを説明できれば、答えが改善され、特に両方のアプローチの違いがコメントされます。


このQ&Aは英語から自動的に翻訳されました。オリジナルのコンテンツはstackexchangeで入手できます。これは、配布されているcc by-sa 3.0ライセンスに感謝します。
Loading...