scikit-bioでfastqを読むための最速の方法

Question

scikit-bioを使用してfastq形式のテキストファイルを読み込もうとしています。

かなり大きなファイルであるため、操作の実行は非常に遅いです。

最終的に、私は辞書にFASTQファイルをdereplicateしようとしています：時間の

f = 'Undetermined_S0_L001_I1_001.fastq' 
seqs = skbio.io.read(f, format='fastq') 

seq_dic = {} 
for seq in seqs: 
    seq = str(seq) 
    if seq in seq_dic.keys(): 
     seq_dic[seq] +=1 
    else: 
     seq_dic[seq] = 1 

ほとんどをここでは、ファイルの読み込み時に使用されます。

%%time 
f = 'Undetermined_S0_L001_I1_001.fastq' 
seqs = skbio.io.read(f, format='fastq') 

for seq in itertools.islice(seqs, 100000): 
    seq 

CPU times: user 46.2 s, sys: 334 ms, total: 46.5 s 
Wall time: 47.8 s 

私の理解では、ということですシーケンスを検証しても実行時間は改善されませんが、そのようには見えません。

%%time 
f = 'Undetermined_S0_L001_I1_001.fastq' 
seqs = skbio.io.read(f, format='fastq', verify=False, variant='illumina1.8') 

for seq in itertools.islice(seqs, 100000): 
    seq 

CPU times: user 47 s, sys: 369 ms, total: 47.4 s 
Wall time: 48.9 s 

私の質問は、なぜverify=Falseが実行時間を改善していないのか、第二にscikit-bioを使ってシーケンスを読み込む方が速いのですか？

Answer 1

最初の理由= Falseの改善の実行時間を確認されていない

verify=False一般scikitバイオのI/O APIによって受け入れられたパラメータです。特定のファイル形式に固有のものではありません。 verify=Falseは、scikit-bioにファイル形式のスニッファを呼び出して、ファイルがユーザーが指定した形式であることを再確認しないように指示します。ドキュメントから、[1]：

verify : bool, optional 
    When True, will double check the format if provided. 

のでverify=Falseはシーケンスデータの検証をオフにしません。ファイル形式のスニファ検証を無効にします。 verify=Falseで最小のパフォーマンス向上が得られます。

seqs = skbio.io.read(f, format='fastq', verify=False, variant='illumina1.8')は、skbio.Sequenceのオブジェクトを生成します。 skbio.Sequenceにアルファベットが含まれていないため、シーケンスのアルファベットの検証は実行されません。そのため、パフォーマンスのボトルネックはどこにありません。 FASTQファイルを特定の種類の生物学的配列（DNA、RNA、またはProtein）に読みたい場合は、たとえばconstructor=skbio.DNAを渡すことができます。 は、関連する配列タイプのアルファベット検証を実行しますが、現在は読み込み時にはトグルオフできません。パフォーマンスの問題があるので、constructorを渡すことはお勧めしません。

2番目の方法は、scikit-bioを使ってシーケンスを読み込む方が速い方法ですか？

scikit-bioでFASTQファイルを読み込む方法はありません。これをスピードアップできるアイデアを探求する問題がありますが、それらのアイデアは実装されていません。

scikit-bioは、複数の行にまたがるシーケンスデータと品質スコアの読み取りをサポートしているため、FASTQファイルの読み込み速度が遅いです。これは読み取りロジックを複雑にし、パフォーマンスが低下します。しかし、複数行のデータを持つFASTQファイルはもう一般的ではありません。イルミナはこれらのファイルの作成に使用されましたが、正確に4行（シーケンスヘッダー、シーケンスデータ、品質ヘッダー、品質スコア）のFASTQレコードの作成を推奨/推奨しています。実際、これはscikit-bioがFASTQデータを書き込む方法です。この簡単なレコードフォーマットでは、FASTQファイルを読むのがはるかに高速で簡単です。scikit-bioは、品質スコアをデコードして検証するため、FASTQファイルの読み込み速度も遅くなります。また、シーケンスデータと品質スコアは、パフォーマンスオーバーヘッドを伴うskbio.Sequenceオブジェクトに格納されます。

あなたの場合、デコードされた品質スコアは必要ありません。単純な4行レコードのFASTQファイルがある可能性があります。ここでFASTQファイルを読み込んで、Pythonの文字列として配列データを生成するPythonの3互換ジェネレータです：

import itertools 

def read_fastq_seqs(filepath): 
    with open(filepath, 'r') as fh: 
     for seq_header, seq, qual_header, qual in itertools.zip_longest(*[fh] * 4): 
      if any(line is None for line in (seq_header, seq, qual_header, qual)): 
       raise Exception(
        "Number of lines in FASTQ file must be multiple of four " 
        "(i.e., each record must be exactly four lines long).") 
      if not seq_header.startswith('@'): 
       raise Exception("Invalid FASTQ sequence header: %r" % seq_header) 
      if qual_header != '+\n': 
       raise Exception("Invalid FASTQ quality header: %r" % qual_header) 
      if qual == '\n': 
       raise Exception("FASTQ record is missing quality scores.") 

      yield seq.rstrip('\n') 

あなたのファイルが正確にされているレコードを含む有効なFASTQファイルで特定している場合は、ここでの検証チェックを削除することができます4行の長さ。

これはあなたの質問とは関係ありませんが、私はあなたのカウンターロジックにバグがあるかもしれないことに気付きました。

if seq in seq_dic: # calling .keys() is necessary 
    seq_dic[seq] +=1 
else: 
    seq_dic[seq] = 1 

[1] http://scikit-bio.org/docs/latest/generated/skbio.io.registry.read.html

[2] https://github.com/biocore/scikit-bio/issues/907

Answer 2

：あなたが初めてのためのシーケンスが表示されたら、あなたのカウンターは、私はロジックがあるべきだと思うの代わりに1のゼロに設定されていますfastqファイル内の各ユニークシーケンスの出現回数をカウントする場合は、標準ライブラリのcollectionsオブジェクトの便利なCounterオブジェクトとともに、BankパーサpyGATBを試してみることをお勧めします。

#!/usr/bin/env python3 

from collections import Counter 
from gatb import Bank 

# (gunzipping is done transparently) 
seq_counter = Counter(seq.sequence for seq in Bank("SRR077487_2.filt.fastq.gz")) 

これは、Pythonモジュール（バイオインフォマティクスSEの私のベンチマークによると：https://bioinformatics.stackexchange.com/a/380/292）のためには非常に効率的です。

このカウンタは、seq_dic,plus some convenient methodsのように動作する必要があります。例えば

、私は一緒にカウントして、10の最も頻繁シーケンスを印刷する場合：この@jairideout

print(*seq_counter.most_common(10), sep="\n")