Осциллирующая скорость обработки в скрипте python с использованием pysam.TabixFile для комментирования показаний

Первоначальный вопрос

Я пишу сценарий биоинформатики в python (3.5), который анализирует большой (отсортированный и проиндексированный) байт- файл, представляющий последовательность чтения, выровненную по геному, связывает геномную информацию («аннотации») с этими чтениями и подсчитывает типы встречающихся аннотаций , Я измеряю скорость, с которой мои скрипты обрабатывают выровненные показания (по партиям 1000 просмотров), и я получаю следующие вариации скорости:

скорость чтения с использованием tabix

Что может объяснить этот шаблон?

Моя интуиция заставила бы меня сделать ставку на некоторую структуру данных, которая постепенно становится медленнее, когда она становится плотнее, но время от времени она расширяется.

Похоже, что использование памяти значимо, хотя (почти через 2 часа мой сценарий по-прежнему использует только 0,1% памяти моего компьютера, в соответствии с htop ).

Как работает мой код (см. В конце для фактического кода)

Я использую модуль pysam для обработки синтаксического анализа. Метод AlignmentFile.fetch дает мне итератор, предоставляющий информацию о последовательных выровненных чтениях в виде объектов AlignedSegment .

Я связываю аннотации с чтениями на основе их координат выравнивания и файла аннотации в формате gtf (сжимается с помощью bgzip и индексируется с помощью tabix ). Я использую метод TabixFile.fetch (все еще из pysam ), чтобы получить эти аннотации, я их фильтрую и frozenset сводку из них в виде frozenset строк ( process_annotations , не показаны ниже, возвращает такой frozenset ), в генераторе функция, которая внутренне пересекает итератор AlignedSegment.

Я питаю генерируемые фенизоны к объекту Counter . Может ли счетчик отвечать за наблюдаемое поведение скорости?

Как я могу узнать, как избежать этого регулярного замедления?


Дополнительные тесты

Следуя рекомендациям в комментариях, я профилировал весь свой анализ с помощью cProfile и обнаружил, что наибольшее время работы было потрачено при доступе к данным аннотации ( pysam/ctabix.pyx:579(__cnext__) , см. График вызовов позже), что, если я понимаю Правильно ли это код Cython, связанный с библиотеками samtools C. Казалось, что причиной наблюдаемого замедления было бы трудно понять.

В попытке ускорить мой скрипт я попробовал другое решение на pybedtools интерфейса python pybedtools с bedtools, которое также может извлекать аннотации из файлов gtf ( https://daler.github.io/pybedtools/index.html ).

скорость

Улучшение скорости очень важно. Вот фактические результаты команды и времени (два фактически выполнялись параллельно):

 $ time python3 -m cProfile -o tests/total_pybedtools.prof ~/src/bioinfo_utils/small_RNA_seq_annotate.py -b results/bowtie2/mapped_C_elegans/WT_1_21-26_on_C_elegans_sorted.bam -g annotations/all_annotations.gtf -a "pybedtools" -l total_pybedtools.log > total_pybedtools.out real 5m48.474s user 5m48.204s sys 0m1.336s $ time python3 -m cProfile -o tests/total_tabix.prof ~/src/bioinfo_utils/small_RNA_seq_annotate.py -b results/bowtie2/mapped_C_elegans/WT_1_21-26_on_C_elegans_sorted.bam -g annotations/all_annotations.gtf.gz -a "tabix" -l total_tabix.log > total_tabix.out real 195m40.990s user 194m54.356s sys 0m47.696s 

(Следует отметить: результаты аннотации немного отличаются друг от друга. Может быть, я должен проверить, как я обрабатываю координаты.)

Профиль скорости не имеет ранее наблюдаемых длиннопериодных капель:

скорость чтения с использованием pybedtools

Моя проблема с скоростью решена, но меня все еще интересует ретроспективный взгляд на то, почему табликский подход имеет эти скорости . По этой причине я добавил тег «bioinformatics» и «samtools».

Графы вызова

Для записи я сгенерировал графики вызовов, используя gprof2dot для результатов профилирования:

 $ gprof2dot -f pstats tests/total_pybedtools.prof \ | dot -Tpng -o tests/total_pybedtools_callgraph.png $ gprof2dot -f pstats tests/total_tabix.prof \ | dot -Tpng -o tests/total_tabix_callgraph.png 

Вот график вызовов для подхода на основе tabix:

График вызова при использовании tabix

Для подхода на основе pybedtools:

График вызова при использовании pybedtools

Код

Вот основная часть моего текущего кода:

 @contextmanager def annotation_context(annot_file, getter_type): """Yields a function to get annotations for an AlignedSegment.""" if getter_type == "tabix": gtf_parser = pysam.ctabix.asGTF() gtf_file = pysam.TabixFile(annot_file, mode="r") fetch_annotations = gtf_file.fetch def get_annotations(ali): """Generates an annotation getter for *ali*.""" return fetch_annotations(*ALI2POS_INFO(ali), parser=gtf_parser) elif getter_type == "pybedtools": gtf_file = open(annot_file, "r") # Does not work because somehow gets "consumed" after first usage #fetch_annotations = BedTool(gtf_file).all_hits # Much too slow #fetch_annotations = BedTool(gtf_file.readlines()).all_hits # https://daler.github.io/pybedtools/topical-low-level-ops.html fetch_annotations = BedTool(gtf_file).as_intervalfile().all_hits def get_annotations(ali): """Generates an annotation list for *ali*.""" return fetch_annotations(Interval(*ALI2POS_INFO(ali))) else: raise NotImplementedError("%s not available" % getter_type) yield get_annotations gtf_file.close() def main(): """Main function of the program.""" parser = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument( "-b", "--bamfile", required=True, help="Sorted and indexed bam file containing the mapped reads." "A given read is expected to be aligned at only one location.") parser.add_argument( "-g", "--gtf", required=True, help="A sorted, bgzip-compressed gtf file." "A corresponding .tbi tabix index should exist.") parser.add_argument( "-a", "--annotation_getter", choices=["tabix", "pybedtools"], default="tabix", help="Method to use to access annotations from the gtf file.") parser.add_argument( "-l", "--logfile", help="File in which to write logs.") args = parser.parse_args() if not args.logfile: logfilename = "%s.log" % args.annotation_getter else: logfilename = args.logfile logging.basicConfig( filename=logfilename, level=logging.DEBUG) INFO = logging.info DEBUG = logging.debug WARNING = logging.warning process_annotations = make_annotation_processor(args.annotation_getter) with annotation_context(args.gtf, args.annotation_getter) as get_annotations: def generate_annotations(bamfile): """Generates annotations for the alignments in *bamfile*.""" last_t = perf_counter() for i, ali in enumerate(bamfile.fetch(), start=1): if not i % 1000: now = perf_counter() INFO("%d alignments processed (%.0f alignments / s)" % ( i, 1000.0 / (now - last_t))) #if not i % 50000: # gc.collect() last_t = perf_counter() yield process_annotations(get_annotations(ali), ali) with pysam.AlignmentFile(args.bamfile, "rb") as bamfile: annot_stats = Counter(generate_annotations(bamfile)) print(*reversed(annot_stats.most_common()), sep="\n") return 0 по @contextmanager def annotation_context(annot_file, getter_type): """Yields a function to get annotations for an AlignedSegment.""" if getter_type == "tabix": gtf_parser = pysam.ctabix.asGTF() gtf_file = pysam.TabixFile(annot_file, mode="r") fetch_annotations = gtf_file.fetch def get_annotations(ali): """Generates an annotation getter for *ali*.""" return fetch_annotations(*ALI2POS_INFO(ali), parser=gtf_parser) elif getter_type == "pybedtools": gtf_file = open(annot_file, "r") # Does not work because somehow gets "consumed" after first usage #fetch_annotations = BedTool(gtf_file).all_hits # Much too slow #fetch_annotations = BedTool(gtf_file.readlines()).all_hits # https://daler.github.io/pybedtools/topical-low-level-ops.html fetch_annotations = BedTool(gtf_file).as_intervalfile().all_hits def get_annotations(ali): """Generates an annotation list for *ali*.""" return fetch_annotations(Interval(*ALI2POS_INFO(ali))) else: raise NotImplementedError("%s not available" % getter_type) yield get_annotations gtf_file.close() def main(): """Main function of the program.""" parser = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument( "-b", "--bamfile", required=True, help="Sorted and indexed bam file containing the mapped reads." "A given read is expected to be aligned at only one location.") parser.add_argument( "-g", "--gtf", required=True, help="A sorted, bgzip-compressed gtf file." "A corresponding .tbi tabix index should exist.") parser.add_argument( "-a", "--annotation_getter", choices=["tabix", "pybedtools"], default="tabix", help="Method to use to access annotations from the gtf file.") parser.add_argument( "-l", "--logfile", help="File in which to write logs.") args = parser.parse_args() if not args.logfile: logfilename = "%s.log" % args.annotation_getter else: logfilename = args.logfile logging.basicConfig( filename=logfilename, level=logging.DEBUG) INFO = logging.info DEBUG = logging.debug WARNING = logging.warning process_annotations = make_annotation_processor(args.annotation_getter) with annotation_context(args.gtf, args.annotation_getter) as get_annotations: def generate_annotations(bamfile): """Generates annotations for the alignments in *bamfile*.""" last_t = perf_counter() for i, ali in enumerate(bamfile.fetch(), start=1): if not i % 1000: now = perf_counter() INFO("%d alignments processed (%.0f alignments / s)" % ( i, 1000.0 / (now - last_t))) #if not i % 50000: # gc.collect() last_t = perf_counter() yield process_annotations(get_annotations(ali), ali) with pysam.AlignmentFile(args.bamfile, "rb") as bamfile: annot_stats = Counter(generate_annotations(bamfile)) print(*reversed(annot_stats.most_common()), sep="\n") return 0 

(Я использовал контекст-менеджер и другие функции более высокого порядка ( make_annotation_processor и функции, которые он вызывает), чтобы упростить использование различных подходов к аннотации в том же скрипте.)

     
    Interesting Posts for Van-Lav

    создание структуры световых данных из многострочной записи

    Связывание значения объектов внутри функции (закрытие)

    Невозможно установить / обновить пакеты с помощью python-pip (« Not a directory»)

    Errno 10060] Не удалось выполнить попытку подключения, потому что связанная сторона не ответила должным образом через некоторое время

    Как перетасовать список с распределением Гаусса

    socket.send работает только один раз в коде python для эхо-клиента

    CherryPy и RESTful web api

    Можем ли мы использовать регулярные выражения для проверки наличия нечетного числа каждого типа символов?

    Опция автоинкремента для индекса Pandas DataFrame

    Pandas Python: объединение двух строк в один фрейм данных

    Python – как получить список самопеременных в классе состоит из N-self

    Использование self. * В качестве значения по умолчанию для метода

    Установленный Virtualenv и активация virtualenv не работают

    python pandas timeseries plot, как установить xlim и xticks за пределами ts.plot ()?

    Обнаружение бесконечной рекурсии на языках Python или динамических языках

    Python - лучший язык программирования в мире.