Самый быстрый способ прочитать fastq с scikit-био
Я пытаюсь прочитать отформатированный текст fastq с помощью http://scikit-bio.org/.
Учитывая, что это довольно большой файл, выполнение операций выполняется довольно медленно.
В конечном итоге я пытаюсь удалить файл fastq в словарь:
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 для чтения последовательностей?
2 ответа
первое почему не проверить = ложное улучшение времени выполнения
verify=False является параметром, принятым API ввода / вывода scikit-bio в целом. Это не относится к конкретному формату файла. 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 в определенный тип биологической последовательности (ДНК, РНК или белок), вы можете передать constructor=skbio.DNA (например). Это выполнит проверку алфавита для соответствующего типа последовательности, который в настоящее время не может быть отключен при чтении. Поскольку у вас проблемы с производительностью, я не рекомендую проходить constructor поскольку это только замедлит вещи еще больше.
а во-вторых, есть ли более быстрый способ использования scikit-bio для чтения последовательностей?
Нет более быстрого способа чтения файлов FASTQ с помощью scikit-bio. Существует проблема [2], исследующая идеи, которые могут ускорить это, но эти идеи не были реализованы.
scikit-bio медленно читает файлы FASTQ, поскольку поддерживает чтение данных последовательности и показателей качества, которые могут занимать несколько строк. Это усложняет логику чтения и снижает производительность. Однако файлы FASTQ с многострочными данными больше не распространены; Illumina раньше создавала эти файлы, но теперь они предпочитают / рекомендуют писать записи FASTQ, которые состоят ровно из четырех строк (заголовок последовательности, данные последовательности, заголовок качества, показатели качества). Фактически, именно так scikit-bio записывает данные FASTQ. С этим более простым форматом записи намного быстрее и проще читать файл FASTQ. scikit-bio также медленно читает файлы FASTQ, поскольку декодирует и проверяет показатели качества. Он также хранит данные последовательности и показатели качества в skbio.Sequence объект, который имеет служебную нагрузку.
В вашем случае вам не нужно декодировать показатели качества, и, скорее всего, у вас есть файл FASTQ с простыми четырехстрочными записями. Вот Python 3-совместимый генератор, который читает файл FASTQ и выдает данные последовательности в виде строк Python:
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, содержащим записи, длина которых ровно четыре строки.
Это не относится к вашему вопросу, но я хотел заметить, что у вас может быть ошибка в вашей контр-логике. Когда вы видите последовательность в первый раз, ваш счетчик устанавливается на ноль вместо 1. Я думаю, что логика должна быть:
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
Если вы хотите посчитать количество вхождений каждой уникальной последовательности в файле fastq, я бы предложил попробовать анализатор банка pyGATB вместе с удобным Counter объект из collections модуль стандартной библиотеки.
#!/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 (согласно моим тестам в Bioinformatics SE: https://bioinformatics.stackexchange.com/a/380/292).
Этот счетчик должен вести себя как ваш seq_dic плюс несколько удобных методов.
Например, если я хочу напечатать 10 самых частых последовательностей вместе с их количеством:
print(*seq_counter.most_common(10), sep="\n")