Giter Site home page Giter Site logo

platonb / high-perf-bio Goto Github PK

View Code? Open in Web Editor NEW
9.0 2.0 0.0 2.69 MB

Open-source toolkit that simplifies and speeds up work with bioinformatics data. high-perf-bio allows you easily upload VCF, BED or arbitrary tables to DB and execute popular bioinformatic queries without MongoDB Query Language knowledge.

License: GNU General Public License v3.0

Python 98.34% Shell 1.66%
bioinformatics mongodb vigg python pymongo vcf bed tsv

high-perf-bio's Introduction

Синопсис.

Биоинформатика — это не только BWA/GATK/Samtools/VEP/продолжите-список-сами, но и беспросветная табличная рутина. Поэтому любой биоинформатик хранит жирную папку со скриптами, а также богатую историю команд, задействующих coreutils. Форумы завалены bash/awk-однострочниками для решения мелких и не очень задач парсинга. Мало того, что парсеров целый зоопарк, так они ещё и медленные все. Есть, правда, tabix и bedtools, которые частично закрывают проблему производительности благодаря самописным B-tree-алгоритмам. Но они обеспечивают быстрый рандомный доступ только по геномным координатам. А если, предположим, надо найти набор rsIDs и отсечь результаты по p-value?

HIGH-PERF-BIO ПОЗВОЛЯЕТ СОЗДАВАТЬ ПОЛНОЦЕННЫЕ БАЗЫ ДАННЫХ НА ОСНОВЕ ТАБЛИЦ И ЗА ДОЛИ СЕКУНДЫ ВЫПОЛНЯТЬ НАИБОЛЕЕ ПОПУЛЯРНЫЕ В БИОИНФОРМАТИКЕ ДЕЙСТВИЯ ☚

Для VCF создаются многоуровневые документы, что открывает возможность глубоко парсить INFO и генотипные поля.

Проиндексировать можно какие угодно поля и группы полей. Индекс редко используется и занимает много места? Можно его удалить, и, если что, потом создать заново отдельным менеджером — reindex.

Флагманские компоненты проекта — три аннотатора (annotate, dock, ljoin), уникальность которых заключается в свободе выбора характеризуемого столбца и отсутствии привязки к определённому кэшу данных. Грубо говоря, можно аннотировать что угодно по чему угодно.

Что касается фильтрации, то в инструментарии имеется удивительная по своей универсальности программа — query. Можете написать любые MongoDB-запросы на целую научную работу вперёд и хлопнуть их разом. Сочинять запросы - одно удовольствие, т.к. синтаксис MongoDB страшно прост и отлично задокументирован.

По состоянию на начало 2022 года я ещё не видел высокоуровневые программы-группировщики данных. Видимо, считается, что суровый дата-сайнтист всякий раз группирует руками в pandas или какой-нибудь СУБД. Но если вдруг хочется одной командой — такое умеют, скорее всего, лишь count и split.

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

Компоненты.

Ядро.

Наиболее универсальные и фичастые программы для взаимодействия с БД.

Программа Основная функциональность
annotate получает характеристики табличного столбца по всем коллекциям БД; поддерживается пересечение координат
concatenate объединяет коллекции; есть опциональная фича удаления дублей документов
count считает количество и частоту элементов каждого набора взаимосвязанных значений разных полей
create создаёт новую или дописывает старую БД по VCF, BED или нестандартным таблицам; производит индексацию коллекций
dock получает характеристики табличного столбца по всем коллекциям и самим таблицам; поддерживается пересечение координат
ljoin фильтрует одни коллекции по наличию/отсутствию значений в других; возможна работа с координатами
info выводит имена или основные свойства имеющихся БД
query выполняет наборы произвольных запросов по всем коллекциям
reindex удаляет и строит индексы коллекций
split дробит коллекции по значениям какого-либо поля

Плагины.

Узкоспециализированные парсеры БД с минимумом настроек.

Программа Основная функциональность
revitalize_id_column восполняет столбец ID VCF-файлов идентификаторами вариантов из БД

Скрипты.

Полезные инструменты для работы с таблицами (не базами).

Программа Основная функциональность
count_lines выдаёт различные метрики, касающиеся количества табличных строк
gen_test_files создаёт из одной таблицы N меньших таблиц

Архитектура.

Что происходит с обрабатываемыми данными с точки зрения программиста: high-perf-bio_architecture

Преимущества.

  • Высокая скорость работы:
    • параллельная обработка нескольких файлов или коллекций;
    • в основе — известная своей отменной производительностью СУБД MongoDB.
  • Простота запуска:
    • самому тулкиту не требуется установка, можно запускать и с флешки;
    • обязательные аргументы - только пути к папкам и имя БД.
  • Наглядность:
    • есть инструмент вывода примерной структуры интересующей БД.
    • для просмотра баз также существует бесплатный MongoDB Compass.
  • Автоматизация:
    • одна команда — обработка всех коллекций базы.
    • полностью автоматическое чтение и создание VCF/BED.
  • Гибкость:
    • индексируемые, а значит быстро парсимые, вложенные поля.
  • Конвейерность:
    • результаты можно перенаправлять в другую БД, минуя создание конечных файлов (но есть баг).
  • Техподдержка:
    • -h: вывод тщательно продуманной справки по синтаксису команд.
    • Тьюториал прямо в этом ридми.
    • Всегда бесплатная консультация.

Недостатки.

  • Высокая производительность возможна лишь при вытаскивании небольших частей коллекций. Преимущество в скорости теряется, если запросить, к примеру, половину dbSNP.
  • Проблемы с производительностью при запросах по нескольким полям, несмотря на задействование индексов.
  • Индексы целиком размещаются в RAM. Готовьтесь к тому, что придётся запихивать в оперативку весящий 15.6 ГБ индекс поля ID dbSNP.
  • Наличие виртуальных ядер процессора не обеспечивает прироста скорости при распараллеливании. hpb_create_benchmark
  • Инструмент ljoin не использует индексы при пересечении/вычитании по геномным координатам. Создатели MongoDB знают о лежащей в основе проблеме, и, надеюсь, когда-нибудь исправят.
  • При скидывании результатов в другую БД сбивается последовательность полей. Голосуйте за исправление этого бага тут (может потребоваться VPN).

Перед началом работы.

Чудо-команда установки.

wget https://github.com/PlatonB/high-perf-bio/archive/refs/heads/master.zip && unzip -o master.zip && rm master.zip && cd high-perf-bio-master && bash __install_3rd_party.sh

Тулкит скачан и распакован, зависимости разрешены. В принципе, теперь можно приступать к эксплуатации. Но при желании ознакомиться с нюансами установки читайте идущие следом разделы.

Установка сторонних компонентов.

Программа Предназначение
MongoDB собственно, размещает данные в базы, индексирует поля и выполняет запросы
PyMongo позволяет сабжевому Python-тулкиту бесшовно работать с MongoDB-базами
Streamlit основа для веб-интерфейса этого инструментария

Установка MongoDB сложновата из-за невероятно абсурдного решения признать новую лицензию MongoDB несвободной, повлекшего не менее идиотское удаление продукта из штатных репозиториев линуксовых дистрибутивов. Но в декабре-2022 в составе тулкита появился скрипт, ставящий MongoDB и другие необходимые high-perf-bio сторонние решения одной командой. Выполните 5 элементарных действий:

  1. Скачайте архив с инструментарием.
  2. Распакуйте его в любую папку.
  3. В терминале перейдите в папку high-perf-bio-master:
cd /path/to/high-perf-bio-master
  1. Разрешите зависимости:
bash __install_3rd_party.sh
  1. Перезагрузитесь, иначе не заработает GUI тулкита.

Удаление сторонних компонентов.

Это делается так же просто, как и установка:

  1. cd /path/to/high-perf-bio-master
  2. bash __uninstall_3rd_party.sh

Примечания.

  • Если хотите установить MongoDB вручную, то на официальном сайте СУБД есть подробные инструкции.
  • Скрипт-инсталлятор, найдя папку miniconda3 в домашнем каталоге, установит Streamlit с помощью Conda. В противном случае будет задействован Pip.
  • Скрипт-инсталлятор ставит PyMongo исключительно через Pip, ибо Conda-версия этого драйвера неработоспособна.
  • Скрипт-деинсталлятор не удаляет Pip-зависимости Streamlit.
  • Скрипт-деинсталлятор беспощадно сносит /etc/mongod.conf и другие конфиги. Если вы их правили, забэкапьте.
  • Windows. Теоретически, после установки MongoDB и whl-пакетов PyMongo + Streamlit программа должна работать. Но у меня сейчас этой ОС нет, и я пока не проверял. Надеюсь, кто-нибудь поделится опытом в Issues.

Произвольное расположение MongoDB-данных.

(опционально)

К примеру, у вас есть до неприличия маленький SSD и вполне просторный HDD (далее - V1). И так случилось, что MongoDB заталкивает базы и логи в первое из перечисленных мест. Как сделать, чтобы эти данные хостились на V1?

  1. Создайте папку для БД и лог-файл.
mkdir /mnt/V1/mongodb
touch /mnt/V1/mongod.log

P.S. Можно было и через файл-менеджер;).

  1. Предельно аккуратно укажите пути к свежеупомянутым папке и файлу в главном конфиге MongoDB:
sudo nano /etc/mongod.conf
<...>
storage:
  dbPath: /mnt/V1/mongodb
<...>
systemLog:
  <...>
  path: /mnt/V1/mongod.log
<...>

Сохраните изменения (CTRL+S) и выйдите из редактора (CTRL+X).

  1. Предоставьте необходимые права доступа нашим переселенцам:
sudo chown mongodb -R /mnt/V1/mongodb
sudo chown mongodb -R /mnt/V1/mongod.log

MongoDB Compass.

(опционально)

Рекомендую устанавливать с Flathub. Compass - родной GUI к MongoDB. В нашем случае он позволяет использовать high-perf-bio не вслепую: делает возможным просмотр создаваемых компонентами тулкита коллекций и индексов. Compass также даёт простор для ручного/продвинутого администрирования БД: составлять запросы, оценивать их производительность, экспериментировать с агрегациями и т.д.. MongoDB_Compass_high-perf-bio

GUI.

В июне 2022 года графический интерфейс на основе Streamlit достиг функционального паритета с консольным. Но, живя в мире биоинформатики, лучше привыкать к командной строке.

Запуск GUI.

Gnome.

Правой кнопкой по __run_gui_streamlit.sh --> Свойства --> Права --> галочку на Разрешить выполнение файла как программы. Это надо сделать лишь однократно. В контекстном меню __run_gui_streamlit.sh появится пункт Запустить как программу. Каждый раз после нажатия на Запустить как программу будет открываться вкладка браузера с GUI к high-perf-bio, а также, с чисто технической целью, окно Терминала. Как закончите работать с тулкитом, просто закройте и то, и другое.

KDE.

Правой кнопкой по __run_gui_streamlit.sh --> Свойства --> Права --> галочку на Является выполняемым --> OK. Правой кнопкой по __run_gui_streamlit.sh --> Открыть с помощью... --> Другое приложение... --> Система --> выделите Konsole --> OK. Вылезет консоль - закройте её. Левой кнопкой по __run_gui_streamlit.sh --> галочку на Больше не задавать этот вопрос --> Запустить. Вышеперечисленное надо сделать лишь один раз. После этих манипуляций GUI будет напрямую запускаться левым кликом по __run_gui_streamlit.sh.

Командная строка.

streamlit run _gui_streamlit.py

Тьюториал.

Основы и игрушечные примеры.

Управление через командный и графический интерфейсы схоже, и здесь для краткости я буду приводить только команды.

Чтобы не прописывать каждый раз путь к тому или иному компоненту, перейдём в основную папку инструментария.

cd $HOME/Биоинформатика/high-perf-bio-master

Примечание: корректируйте приведённые в примерах пути в зависимости от реального расположения ваших папок или файлов.

Для удобства можете вывести имена всех программ проекта high-perf-bio.

ls -R

Каждый инструмент содержит подробную справку. Пока с ней не ознакомились, реальные задачи не запускайте.

python3 create.py -h
python3 annotate.py -h

и т.д..

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

В папке high-perf-bio-master уже есть небольшие примеры данных. Но в качестве материала для базы лучше скачать что-то серьёзное, например, VCF-таблицу всех SNP. Переименуйте её, добавив расширение vcf перед gz.

Поскольку формат этого файла — VCF, программа создания БД сама отбросит метастроки, сконвертирует хромосомы из HGVS в обычный вид, преобразует ячейки определённых столбцов в BSON-структуры и проиндексирует поля с хромосомами, позициями и rsIDs. Файл — один, значит, распараллелить заливку не удастся. Таким образом, команда может получиться весьма минималистичной. Достаточно подать путь к папке со свежескачанным VCF.

python3 create.py -S $HOME/Биоинформатика/dbSNP

При работе с high_perf_bio высока вероятность возникновения двух проблем: вы забыли, как называется БД и из чего она состоит. В таких ситуациях перед запуском любого парсера требуется получить представление о парсимой БД. Для начала найдите её имя в списке всех имён.

python3 info.py

Затем выведите первостепенные характеристики конкретной базы.

python3 info.py -d dbSNP

Теперь мы многое знаем о БД и можем, к примеру, осуществить по ней запрос. Для разминки отберём снипы интервала chr14:105746998-105803861, упомянутого в статье "Локусы влияющие на экспрессию антигенов HLA в участке 14-й хромосомы, ассоциированном с развитием рассеянного склероза, и функции расположенных в них генов" (П.А. Быкадоров и др.). Создадим текстовый файл get_ms_locus.txt, впишем туда запрос и выполним его программой query. Помните, что операторы и строковые значения надо заключать в кавычки (неважно, одинарные или двойные), а целочисленные значения оставлять без ничего. Содержимое get_ms_locus.txt:

{'#CHROM': 14, 'POS': {'$gte': 105746998, '$lte': 105803861}}

Команда, которая прочитает запрос из файла и выполнит его:

python3 query.py -S $HOME/Биоинформатика/toy_queries/get_ms_locus.txt -D dbSNP -T $HOME/Биоинформатика/ms_locus

Кстати, инструмент query гораздо мощнее, чем могло показаться по примеру выше. Можно задавать не один, а сколько угодно запросов, а также тематически распределять их по разным файлам.

Ещё, чтобы лучше привыкнуть к high_perf_bio, найдём характеристики вариантов из папки high-perf-bio-master/test_data/TSV. С расположением rsID-столбца повезло — он первый. С единственной коллекцией базы dbSNP тоже всё ок — поскольку это ex-VCF, то там есть поле ID. Значит, обойдёмся минимумом аргументов: пути к папке с игрушечной таблицей и конечной папке, имя базы и игнор табличной шапки.

python3 annotate.py -S $HOME/Биоинформатика/high-perf-bio-master/test_data/TSV -D dbSNP -T $HOME/Биоинформатика/TSV_ann -m 1

Пример чуть посерьёзнее.

А вот задача явно ближе к реальной биоинформатике. Вы когда нибудь задумывались, есть ли на свете варианты, являющиеся eQTL сразу во всех тканях? Сейчас посмотрим. Скачаем таблицы пар вариант-ген. Нам понадобятся только *.signif_variant_gene_pairs.txt.gz-файлы. Таблицы так себе приспособлены к анализу: координаты, аллели и номер сборки запихнуты в один столбец, а rsIDs вовсе отсутствуют. Сконвертировать это непарсибельное хозяйство в VCF предлагаю скриптом, а обогатить VCFs сниповыми идентификаторами — high-perf-bio-плагином revitalize_id_column. Дарить rsIDs VCF-файлам готова ранее созданная нами MongoDB-версия dbSNP. Насыщение поля ID идентификаторами может оказаться неполным. Чтобы не оставалось строк с точками, добавим аргумент -r.

cd plugins
python3 revitalize_id_column.py -S $HOME/Биоинформатика/GTEx_Analysis_v8_eQTL_VCF -D dbSNP -T $HOME/Биоинформатика/GTEx_Analysis_v8_eQTL_VCF_rsIDs -r

Зальём удобоваримую версию GTEx-таблиц в базу. Коллекций — несколько десятков, поэтому лично я предпочёл бы распараллелить размещение и индексацию на 8 процессов.

python3 create.py -S $HOME/Биоинформатика/GTEx_Analysis_v8_eQTL_VCF_rsIDs -d GTEx -p 8

Напомню, мы хотели найти варианты, дающие о себе знать во всех тканях. Применим алгоритм пересечения. Что такое левые/правые коллекции и охват, подробно рассказано в справке ljoin — здесь я на этом останавливаться не буду. Указываем в явном виде одну левую коллекцию. Любую. Правыми коллекциями для решения задачи надо сделать все остальные, но перечислять их не нужно — в программе такой сценарий предусмотрен по умолчанию. То, что нам нужно пересекать, а не вычитать, опять же, явно прописывать не обязательно. С наличием идентификаторов у нас проблем нет, поэтому можно работать по ним, хотя координатные вычисления программой тоже поддерживаются. Подавать аргумент с именем пересекаемого поля не будем, т.к. по дефолту идёт ID. eQTL нам общие для всех тканей надо найти? Значит, задаём охват максимальный. Он равен количеству всех коллекций минус 1 (вычли одну левую). Программа, кстати, позволяет исследователю не тратить время на определение этого значения, а указать 0.

python3 ljoin.py -D GTEx -T $HOME/Биоинформатика/GTEx_common -l Adipose_Subcutaneous.v8.signif_variant_gene_pairs_rsIDs.vcf -c 0

Вездесущие eQTL обнаружились, и их 224045 штук. А сколько будет вариантов, влияющих на экспрессию генов только в крови и больше нигде? Это решается вычитанием.

python3 ljoin.py -D GTEx -T $HOME/Биоинформатика/GTEx_onlyWB -l Whole_Blood.v8.signif_variant_gene_pairs_rsIDs.vcf -a subtract -c 0

Таких уникальных вариантов — 86051 из 2399784 кровяных.

Найти что ли клинически значимые варианты среди эксклюзивно кровяных? Воспользуемся таблицей ассоциаций вариант-патология проекта DisGeNET. Эта таблица — небиоинформатического формата, поэтому при создании на её основе базы и дальнейшей работе с этой базой применимы не все удобные умолчания.

python3 create.py -S $HOME/Биоинформатика/DisGeNET -s semicolon -i snpId,chromosome,position

Можно, кстати, из этой БД создать новую, которая разбита на похромосомные коллекции, урезанные до 5 самых необходимых полей.

python3 split.py -D DisGeNET -T DisGeNET_chrs -f chromosome -k snpId,chromosome,position,diseaseName,score -i snpId,chromosome,position

Ой, кажется, забыли проиндексировать поле со скорами ассоциаций. С кем не бывает? Воспользуемся штатным редактором индексов, чтобы потом программы, отсекающие малоподтверждённые пары, не тратили силы на тупой перебор скоров.

python3 reindex.py -D DisGeNET_chrs -a score

Теперь, наконец, аннотируем уникальные Whole Blood eQTLs по DisGeNET. Результаты кидаем в свежую БД.

python3 annotate.py -S $HOME/Биоинформатика/GTEx_onlyWB -D DisGeNET_chrs -T GTEx_onlyWB_Clin -f snpId -i snpId,chromosome,position,score

Хоть какое-то упоминание в DisGeNET имеют 1914 кровяных вариантов. Отбирём наиболее клинически значимые. Как в DisGeNET присваиваются скоры? Если ассоциация упоминается в отобранных нейросетевым парсером BeFree статьях, ей прибавляют скор 0.01 * количество_статей, но не более 0.1. Если ассоциация взята из одного курируемого источника (UniProt, ClinVar, GWAS Catalog, GWASdb) - то за ней сразу бронируется аж 0.7. Два курируемых источника - 0.8, три и более - 0.9. Хочется, конечно, выставить 0.9, но тогда ничего не останется. Сделаем 0.81, чтобы связь подтверждалась, как минимум, двумя солидными базами и одной статьёй.

{'score': {'$gte': Decimal128('0.81')}}
python3 query.py -S $HOME/Биоинформатика/queries/parse_GTEx.txt -D GTEx_onlyWB_Clin -T GTEx_onlyWB_strongClin

Результаты разбросаны по хромосомам, а названия болезней иногда дублируются. Как тогда собрать список болезней? Очень просто: сконкатенировать коллекции, применив подходящие опции. Поскольку мы используем отбор полей, -u уникализирует не целиковые документы, а усечённые с помощью -k.

python3 concatenate.py -D GTEx_onlyWB_strongClin -T GTEx_onlyWB_081Dis -k diseaseName -u

Результат. Варианты, влияющие на экспрессию генов исключительно в крови, обуславливают развитие, вероятнее всего, перечисленных 10 патологий: "Psoriasis", "Inflammatory Bowel Diseases", "Vitiligo", "Crohn Disease", "Diabetes Mellitus, Non-Insulin-Dependent", "Age related macular degeneration", "ovarian neoplasm", "Behcet Syndrome", "IGA Glomerulonephritis", "Hirschsprung Disease".

Расскажите о применении high_perf_bio в своих исследованиях!

Механизм пересечения и вычитания.

Этот раздел — попытка максимально доступно объяснить, как работает программа ljoin.

Терминология.

Термин Аналог в реляционных БД Аналог в Python Пример Пояснение
Поле Столбец Ключ словаря 'f1' Может присутствовать в одних документах и при этом отсутствовать в других
Документ Строка Словарь {'f1': 1, 'f2': 'one', 'f3': 'c1.txt'} Объект из пар поле-значение. В качестве значений могут быть вложенные объекты
Коллекция Таблица Список словарей Набор, как правило, объединённых по смыслу документов

Тестовые данные.

В качестве примера будет БД из трёх игрушечных коллекций, которые легко пересекать и вычитать в уме. _id-поле здесь и далее опускается.

Коллекция c1.

{'f1': 1, 'f2': 'one', 'f3': 'c1.txt'}
{'f1': 2, 'f2': 'two', 'f3': 'c1.txt'}
{'f1': 3, 'f2': 'three', 'f3': 'c1.txt'}
{'f1': 4, 'f2': 'four', 'f3': 'c1.txt'}
{'f1': 5, 'f2': 'five', 'f3': 'c1.txt'}

Коллекция c2.

{'f1': 3, 'f2': 'three', 'f3': 'c2.txt'}
{'f1': 4, 'f2': 'four', 'f3': 'c2.txt'}
{'f1': 5, 'f2': 'five', 'f3': 'c2.txt'}
{'f1': 6, 'f2': 'six', 'f3': 'c2.txt'}
{'f1': 7, 'f2': 'seven', 'f3': 'c2.txt'}

Коллекция c3.

{'f1': 5, 'f2': 'five', 'f3': 'c3.txt'}
{'f1': 6, 'f2': 'six', 'f3': 'c3.txt'}
{'f1': 7, 'f2': 'seven', 'f3': 'c3.txt'}
{'f1': 8, 'f2': 'eight', 'f3': 'c3.txt'}
{'f1': 9, 'f2': 'nine', 'f3': 'c3.txt'}

Параметры отбора документов левой коллекции.

По справке к ljoin вы уже знаете, что каждая из коллекций, условно обозначенных как "левые", фильтруется по всем "правым". Левой пусть будет c1, а в качестве правых — c2 и c3. Непосредственно обрабатываемое поле — допустим, f1. Общее представление об условях отбора документов c1:

Действие Охват Документ из c1 сохраняется, если
пересечение 1 совпадение в c2 или c3
пересечение 2 совпадение в c2 и c3
вычитание 1 несовпадение в c2 или c3
вычитание 2 несовпадение в c2 и c3

Решение.

Стадии.

  1. Левостороннее внешнее объединение. Независимо от того, хотим мы пересекать или вычитать, выполнение $lookup-инструкции лишь объединит коллекции. Каждый документ, получающийся в результате объединения, содержит: 1. все поля документа левой коллекции; 2. поля с соответствиями из правых коллекций (далее — результирующие). Если для элемента выбранного поля данного левого документа не нашлось совпадений в одной из правых коллекций, то в результирующем поле появится пустой список (Python-представление Null-значения из мира баз данных). Если же совпадение имеется, то в качестве значения результирующего поля возвратится список с содержимым соответствующего правого документа. Если выявилось совпадение с несколькими документами какой-то одной правой коллекции, то они в полном составе поступят в результирующее поле.

Схема объединённого документа.

{левый документ, псевдоним правой коллекции 1: [правый документ 1.1, правый документ 1.2, ...],
                 псевдоним правой коллекции 2: [правый документ 2.1, правый документ 2.2, ...],
                 ...}

Объединённые документы тестовых коллекций (c1 vs c2, c3).

{'f1': 1, 'f2': 'one', 'f3': 'c1.txt', 'c2': [], 'c3': []}
{'f1': 2, 'f2': 'two', 'f3': 'c1.txt', 'c2': [], 'c3': []}
{'f1': 3, 'f2': 'three', 'f3': 'c1.txt', 'c2': [{'f1': 3, 'f2': 'three', 'f3': 'c2.txt'}], 'c3': []}
{'f1': 4, 'f2': 'four', 'f3': 'c1.txt', 'c2': [{'f1': 4, 'f2': 'four', 'f3': 'c2.txt'}], 'c3': []}
{'f1': 5, 'f2': 'five', 'f3': 'c1.txt', 'c2': [{'f1': 5, 'f2': 'five', 'f3': 'c2.txt'}], 'c3': [{'f1': 5, 'f2': 'five', 'f3': 'c3.txt'}]}

То же самое в схематичном виде (l — документ левой коллекции, r№a — alias (псевдоним) правой; 0 — в правой нет совпадений, 1 — есть одно).

l, r1a: 0, r2a: 0
l, r1a: 0, r2a: 0
l, r1a: 1, r2a: 0
l, r1a: 1, r2a: 0
l, r1a: 1, r2a: 1
  1. Фильтрация. Изучая предыдущий пункт, вы уже, наверное, догадались, что пересекать или вычитать — значит, фильтровать результаты объединения по количеству непустых или пустых результирующих списков соответственно.

Пересечение. Охват 1.

Чтобы значение левого поля f1 попало в результаты, такое же значение должно найтись хотя бы в одном правом f1.

Объединённые документы, отвечающие условию.

{'f1': 3, 'f2': 'three', 'f3': 'c1.txt', 'c2': [{'f1': 3, 'f2': 'three', 'f3': 'c2.txt'}], 'c3': []}
{'f1': 4, 'f2': 'four', 'f3': 'c1.txt', 'c2': [{'f1': 4, 'f2': 'four', 'f3': 'c2.txt'}], 'c3': []}
{'f1': 5, 'f2': 'five', 'f3': 'c1.txt', 'c2': [{'f1': 5, 'f2': 'five', 'f3': 'c2.txt'}], 'c3': [{'f1': 5, 'f2': 'five', 'f3': 'c3.txt'}]}

Схема соответствующих условию объединённых документов.

l, r1a: 1, r2a: 0
l, r1a: 1, r2a: 0
l, r1a: 1, r2a: 1

Окончательный результат в табличном виде.

f1 f2 f3
3 three c1.txt
4 four c1.txt
5 five c1.txt

Пересечение. Охват 2.

Чтобы значение левого поля f1 попало в результаты, такое же значение должно найтись в двух правых f1.

Объединённый документ, отвечающий условию.

{'f1': 5, 'f2': 'five', 'f3': 'c1.txt', 'c2': [{'f1': 5, 'f2': 'five', 'f3': 'c2.txt'}], 'c3': [{'f1': 5, 'f2': 'five', 'f3': 'c3.txt'}]}

Схема соответствующего условию объединённого документа.

l, r1a: 1, r2a: 1

Окончательный результат в табличном виде.

f1 f2 f3
5 five c1.txt

Вычитание. Охват 1.

Чтобы значение левого поля f1 попало в результаты, такого же значения не должно быть хотя бы в одном правом f1.

Объединённые документы, отвечающие условию.

{'f1': 1, 'f2': 'one', 'f3': 'c1.txt', 'c2': [], 'c3': []}
{'f1': 2, 'f2': 'two', 'f3': 'c1.txt', 'c2': [], 'c3': []}
{'f1': 3, 'f2': 'three', 'f3': 'c1.txt', 'c2': [{'f1': 3, 'f2': 'three', 'f3': 'c2.txt'}], 'c3': []}
{'f1': 4, 'f2': 'four', 'f3': 'c1.txt', 'c2': [{'f1': 4, 'f2': 'four', 'f3': 'c2.txt'}], 'c3': []}

Схема соответствующих условию объединённых документов.

l, r1a: 0, r2a: 0
l, r1a: 0, r2a: 0
l, r1a: 1, r2a: 0
l, r1a: 1, r2a: 0

Окончательный результат в табличном виде.

f1 f2 f3
1 one c1.txt
2 two c1.txt
3 three c1.txt
4 four c1.txt

Вычитание. Охват 2.

Чтобы значение левого поля f1 попало в результаты, такого же значения не должно быть в двух правых f1.

Объединённые документы, отвечающие условию.

{'f1': 1, 'f2': 'one', 'f3': 'c1.txt', 'c2': [], 'c3': []}
{'f1': 2, 'f2': 'two', 'f3': 'c1.txt', 'c2': [], 'c3': []}

Схема соответствующих условию объединённых документов.

l, r1a: 0, r2a: 0
l, r1a: 0, r2a: 0

Окончательный результат в табличном виде.

f1 f2 f3
1 one c1.txt
2 two c1.txt

Повторы.

В конечный файл направятся все дубли элемента, находящиеся в пределах поля левой коллекции, но от повторов правых элементов копийность результатов не зависит.

Задавайте вопросы!

Понятно, что реальные данные гораздо сложнее, чем раз-два-три-четыре-пять. Кидайте свои задачи на пересечение или вычитание в Issues — подскажу, как решить.

Полезные советы.

Метастроки.

У части тулзов high_perf_bio есть опция --meta-lines-quan/-m. Она даёт программе знать, сколько строк в начале файла не нужно никак трогать. Эти строки содержат т.н. метаинформацию, которая посвящает исследователя в различные особенности файла. Программам же они, как правило, только мешают. VCF — формат довольно строгий, и там метастроки чётко обозначены символами ##, что даёт возможность скипать их автоматически. Для остальных форматов требуется ручное указание количества игнорируемых строк. Как это количество узнать? Если файл маленький, просто откройте его в обычном блокноте. В хороших блокнотах, как, например, elementary OS Code, есть нумерация строк, что облегчает задачу. Если же файл большой, то так сделать не получится — в лучшем случае вылезет ошибка, в худшем — осложните себе работу зависанием. Вас раздражает командная строка? Можете считывать первые строки сжатого файла скриптом-предпросмотрщиком из моего репозитория bioinformatic-python-scripts. Уже привыкли к эмулятору терминала? Тогда юзайте командную утилиту zless.

zless -N $HOME/Биоинформатика/dbSNP_common/00-common_all.vcf.gz

Скроллить можно колесом мыши или двумя пальцами вверх-вниз по тачпаду. Закрыть предпросмотр: q.

Спецсимволы.

В ту или иную команду могут закрасться зарезервированные символы. Командная оболочка, на них наткнувшись, не сможет обработать вашу команду как единое целое. Я, к примеру, при тестировании high_perf_bio получил ошибку из-за решётки в начале имени хромосомного поля.

python3 create.py -S $HOME/Биоинформатика/dbSNP_common -i #CHROM,POS,ID
create.py: error: argument -i/--ind-col-names: expected one argument

Из-за наличия решётки интерпретатор подумал, что CHROM,POS,ID — комментарий, а не перечисление индексируемых полей. Получилось, что аргумент -i, строго требующий значения, остался без такового. Надо было просто взять набор имён полей в одинарные кавычки (заэкранировать).

python3 create.py -S $HOME/Биоинформатика/dbSNP_common -i '#CHROM,POS,ID'

А знаете, почему программы high_perf_bio заставляют перечислять что-либо через запятую без привычного нам по естественным языкам пробела? Не подумайте, что из вредности:). Пробел считается границей между аргументами. Т.е., при использовании запятых с пробелом каждый следующий элемент воспринимался бы шеллом в качестве нового аргумента.

python3 create.py -S $HOME/Биоинформатика/high-perf-bio-master/test_data/TSV -i SYMBOL, AF
create.py: error: unrecognized arguments: AF

Требование перечислять без пробела я ввёл, чтобы избавить вас от необходимости экранирования. Правильное оформление последней команды не подразумевает обязательного наличия надоедливых кавычек.

python3 create.py -S $HOME/Биоинформатика/high-perf-bio-master/test_data/TSV -i SYMBOL,AF

high-perf-bio's People

Contributors

platonb avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar

high-perf-bio's Issues

Мягкая проверка правильности ввода имени поля

Жёсткая проверка правильности ввода имён полей путём поиска по ключам первого документа первой коллекции иногда всё-таки приводит к недоразумениям. К примеру, в первом документе dbSNP ex-VCF, насколько я помню, нет GENEINFO (или PSEUDOGENEINFO - лень проверять), и если попытаться выполнить по нему запрос, вылезет ошибка, что, мол, этого поля нет во всей коллекции. От проверки отказываться не буду, но заменю исключение на предупреждение. Понадобится модуль warnings стандартной библиотеки.

Сохранять результаты в виде коллекций

Сейчас результаты выводятся в виде текстовых файлов. Для дальнейшего анализа их приходится вручную заливать в ещё одну БД. Надо бы сделать функциональность сохранения в базу напрямую.

Реализация:

  • Придумывать без особой необходимости новые опции я не люблю. Здесь есть прекрасная возможность обойтись имеющейся опцией -t. Любой путь (даже к домашней папке) содержит слэши. Сделаю так: значение -t без слэшей будет восприниматься как имя новой БД для результатов. Проработаю поведение на случай существования базы с тем же именем.
  • Сейчас в текстовые аутпуты добавляются строки с их историей. При перегонке из БД в БД нужно делать нечто аналогичное. Для метаданных можно смастерить либо отдельную коллекцию базы, либо особый документ каждой коллекции. При решении вопроса менеджмента метаданных хорошо бы заодно подумать над #9.
  • Обязательно сохранять в именах коллекций-потомков расширение файлов, по которым были получены родительские коллекции.
  • Индексировать новоиспечённые коллекции таким же образом, как это делает create_db.

[_gui_streamlit] Создать запускабельный двойным кликом .desktop-файл

Сабж представляю себе таким:

[Desktop Entry]
Type=Application
Version=1.0
Name=high-perf-bio
Comment=Open-source toolkit that simplifies and speeds up work with bioinformatics data

#Replace SPICIFY/PATH/TO with real path:
Path=SPICIFY/PATH/TO/high-perf-bio-master

Exec=bash -c 'streamlit run _gui_streamlit.py'
Icon=utilities-terminal
Terminal=true
Categories=Application;

#Next steps:
#cd SPICIFY/PATH/TO/high-perf-bio-master
#sudo cp _gui_streamlit.desktop /usr/share/applications/

Но, как минимум, в elementary OS двойной клик лишь открывает _gui_streamlit.desktop в блокноте. Если переместить файл в /usr/share/applications/, исполняться он от этого не научится.

Terminal=true, в Exec прописана работоспособная команда... Что ещё нужно для счастья?

[dock] Вывод каркаса запроса в конечную метастроку

Для каждой строки аннотируемой таблицы формируется свой запрос, но выделить общую часть всех ежестрочных запросов всё же возможно. С появлением функциональности дозапроса становится очевидной необходимость вывода метастроки с запросом, пусть и не до конца сформированным.

[create_db] Поддержка JSON Lines

Некоторые важные данные, например, источник соответствия rsIDs разных версий dbSNP, доступны в виде сабжа.

Реализиция:

  • Поскольку каждый элемент JSON Lines идёт отдельной строкой, то считывание будет обычное - построчное, с накоплением в списки заданного размера.
  • С учётом того, что объекты могут быть непредсказуемой вложенности, не понятно, как делать автоопределение типов данных. Скорее всего, придётся полагаться на правильность определения со стороны авторов JSON-файлов.
  • Квазирасширение надо давать, очевидно, TSV.

[Streamlit GUI] Выбор полей колекций через st.multiselect

Раз GUI-модули теперь подгружают MongoClient (#30), то открывается возможность сделать мышко-ориентированный выбор монгячих полей. Поскольку st.multiselect при невыборе выдаёт пустой список, а не пустую строку, понадобится учесть это в инитах соответствующих главных модулей.

[create] Поддержка CSQ-INFO в VCF

Сейчас может читаться и создаваться INFO-столбец, оформленный по стандарту ANN. Но существует и другой стандарт - CSQ (ConSeQuence annotations). В идеале надо, чтобы он тоже поддерживался.

[plugins] Программа пополнения ID-столбца src-VCF

Все компоненты high-perf-bio ориентированы на фильтрацию данных, размещённых в базы. Допущенное мною исключение - конкатенатор коллекций, да и там тоже возможна фильтрация ($project). Проги, сильно выходящие за рамки этого принципа, я буду размещать в отдельном каталоге - plugins. Планируемый обогатитель ID-столбца rs-идентификаторами, скорее всего, станет первым таким инструментом. Результаты будут формироваться не из коллекций, а из обычного VCF-файла. БД понадобится программе лишь как источник соответствия координат и аллелей с rsIDs.

Скрипт подсчёта строк

Уже есть такой продукт в старом скриптовом проекте. Но хочу переписать с нуля и переосмыслить функциональность.

  • CLI и GUI.
  • Стиль кода по стандартам high-perf-bio.
  • Количество строк для каждого файла, для всех файлов каждой папки и, возможно, для всех файлов дерева папок.
  • Опциональный подсчёт не полных строк, а значений столбца или нескольких столбцов.
  • Сохранение или удаление дублей.

[ljoin] $lookup не задействует индекс при работе с интервалами

Вчера, по-сути, завершил проект high-perf-bio, реализовав в left_join отбор геномных интервалов одной MongoDB-коллекции по их вхождению/невхождению в интервалы из других коллекций БД. Но, как выяснилось, не без ложки дёгтя.

Пояснение для IT-специалистов, не специализирующихся на работе с генетическими данными. Что такое геномный интервал? Есть хромосомы - надмолекулярные сткуктуры, включающие, помимо всего прочего, ДНК. Есть нуклеотиды - мономеры ДНК. Сослаться на тот или иной участок ДНК можно так: обозначение хромосомы, номер стартового нуклеотида, номер конечного нуклеотида. Эти два нуклеотида - границы как раз геномного интервала. В биоинформатических таблицах чаще всего мы видим интервалы с выявленным влиянием на организм.

Aggregation pipeline из программы left_join, построенный по официальным докам MongoDB, выдаёт правильные результаты, но работает бесконечно долго для больших коллекций. Ни compound, ни одиночные индексы не ускоряют вычисление. Вот основа пайплайна - код левостороннего внешнего объединения коллекций - источника описываемой проблемы:

pipeline = [{'$lookup': {'from': right_coll_name,
                         'let': {'chrom': '$chrom', 'start': '$start', 'end': '$end'},
                         'pipeline': [{'$match': {'$expr': {'$and': [{'$eq': ['$$chrom', '$chrom']},
                                                                     {'$lt': [{'$max': ['$$start', '$start']},
                                                                              {'$min': ['$$end', '$end']}]}]}}}],
                         'as': right_coll_name.replace('.', '_')}} for right_coll_name in right_coll_names]

Эта же конструкция, но чуть упрощённая (Ilya Vorontsov):

pipeline = [{'$lookup': {'from': right_coll_name,
                         'let': {'chrom': '$chrom', 'start': '$start', 'end': '$end'},
                         'pipeline': [{'$match': {'$expr': {'$and': [{'$eq': ['$$chrom', '$chrom']},
                                                                     {'$lt': ['$$start', '$end']},
                                                                     {'$lt': ['$start', '$$end']}]}}}],
                         'as': right_coll_name.replace('.', '_')}} for right_coll_name in right_coll_names]

Каждый объединённый документ представляет собой документ "левой" коллекции, в который вложены отвечающие запросу документы "правых" коллекций. Если в "правой" коллекции не нашлось соответствий, в объединённый документ попадает пустой список.

Кто-нибудь знает, как решить/обойти проблему игнора индексов?

P.S. Я уже создал аналогичную тему на официальном форуме MongoDB, но там далеко не всегда отвечают. Поэтому очень надеюсь на помощь коллег, знакомых и других посетителей репозитория.

@VorontsovIE, @yustinaivanova, может у вас будут какие-то идеи? Был бы очень благодарен.

Предполагаемые тормозящие факторы:

  • Вложенный пайплайн, какое бы выражение в нём ни было. UPD: точно нет: если попробовать объединять аналогичной конструкцией не по интервалам, а по одиночным полям, то всё посчитается почти мгновенно.
  • Подсчёт максимальной стартовой границы и минимальной конечной (Ilya Vorontsov).
  • Сочетание работы по хромосомам с работой по интеравалам.
  • Добавляемая к $lookup-пайплайну сортировка. Сразу говорю, что она не влияет.

Создать Dockerfile

Зачем?

  • Не на все дистрибутивы легко установить MongoDB. А для rpm-ostree-based систем вообще не понятно, как это делать.
  • Чтобы избавить вас от любой предэксплуатационной возни в принципе.

Недостаток использования high-perf-bio как Docker-контейнера: команды запуска компонентов распухнут из-за необходимости монтирования томов.

Превьюер MongoDB-коллекций

Для оценки, из чего состоит база, как-то недостаточно вывода имён полей и индексов. Надо иметь возможность демонстрации целых документов.

Варианты реализации:

  • Часть print_db_info. Эта фича, по-моему, неплохо сюда вписывается.
  • Часть make_request. Возможно, сделаю опцию верхнего порога количества отобранных документов. Выставление небольшого значения будет равносильно предпросмотру.
  • Отдельная программа. Вряд ли, т.к. я выступаю за минимализацию количества программ в пределах одного репозитория.

Программа развёртывания и свёртывания src-db-VCF по полю ALT

Раз задумка касается только VCF, то проге быть плагином.

Прототип:
MongoDB_unwind_group

Если с $unwind-стадией всё просто, то с формированием $group придётся чуть повозиться:

  1. Собрать имена всех полей (N штук) исходного VCF в список;
  2. Добавить в $group ID в качестве _id;
  3. Сгенерировать N $first-выражений, вставив туда вышеупомянутые имена полей;
  4. Добавить эти выражения в $group;
  5. Для ALT заменить $first на $push.

Вынести повторяющиеся от компонента к компоненту участки __init__ в отдельный модуль

Код инитов тогда существенно разгрузится.

Кандидаты на вынос:

  • Конфликт имён БД.
  • Количество процессов.
  • Дополнительный запрос.
  • Сортировка.
  • Отбор полей.
  • Вторичный разделитель.
  • Индексация.

Возможно, стоит также перенести в этот модуль функцию parse_nested_objs.

Переименовать args в inp

Пусть это будет единым для CLI и GUI именем объектов, накапливающих пользовательский ввод.

[create] Недостаточно глубоко парсится INFO dbSNP

К примеру, строка

RS=1570391677;dbSNPBuildID=154;SSR=0;PSEUDOGENEINFO=DDX11L1:100287102;VC=SNV;GNO;FREQ=KOREAN:0.9891,0.0109|SGDP_PRJ:0,1;COMMON

преобразуется в

    "INFO": [{
            "RS": 1570391677,
            "dbSNPBuildID": 154,
            "SSR": 0,
            "PSEUDOGENEINFO": "DDX11L1:100287102",
            "VC": "SNV",
            "FREQ": ["KOREAN:0.9891", "0.0109|SGDP_PRJ:0", 1]
        },
        ["GNO", "COMMON"]
    ]

Очевидно, что значения подполей вроде PSEUDOGENEINFO и FREQ должны делаться словарями.

Сортировка проиндексированных полей

Проиндексированные поля, вероятнее всего, являются наиболее значимыми для исследователя. Хорошо бы выдавать таблицы, отсортированные по соответствующим столбцам. Т.к. поля - с индексами, на производительность сортировка не должна влиять.

Игнорируются строки метаинформации VCF

Сейчас они выбрасываются ещё на этапе создания БД и не могут оказаться в результатах парсинга коллекций. Без строк метаинформации получающиеся таблицы не полностью соответствуют стандарту VCF. Надо бы эти строки протаскивать в базу и кидать в конечные файлы. Но остаются вопросы, требующие углубления в спецификации формата.

  • На сегодняшний день компоненты high_perf_bio прописывают в конечные файлы несколько закомментированных строк, хранящих в себе историю этих файлов. Позволительно ли их сочетать с генерируемыми коллерами метастроками VCF?
  • Программы, разумеется, выдают не полные VCF, а лишь строки, соответствующие запросам. Актуальны ли для фильтрованного VCF метастроки залитого в базу VCF?

При выводе результатов в другую БД нарушается порядок полей

Для VCF, BED и других биоинформатических форматов страшнее бага не придумать. Уже сообщил разработчиками MongoDB:
https://jira.mongodb.org/browse/SERVER-63853
Проголосуйте там, плз, откоментьте.

UPD.
Сейчас актуальнее прогососовать тут:
https://feedback.mongodb.com/forums/924280-database/suggestions/44875462-preserve-field-order-in-merge
(хотя по старой ссылке тоже не будет лишним)

Программа пересечения и вычитания BED-интервалов

Уже есть подобная функциональность у bedtools, но хочется, чтобы были:

  • выше скорость;
  • меньше потребление RAM;
  • настройка глубины.

Отличие от bedtools (UPD-13.04.2020):
Каждый конечный файл - BED по левой коллекции с немодифицированными интервалами. Т.е. пересечение и вычитание применяется исключительно для отбора данных. bedtools же по умолчанию выводит пересёкшиеся или непересёкшиеся области исходных интервалов, что не сочеталось бы с настройкой глубины intersect_subtract.

Реализация:

  • Как часть intersect_subtract. UPD-13.04.2020: опасения, что программа станет слишком нагромождённой, были напрасны - интервальная функциональность хорошо вписалась в имеющийся код и не потребовала создавать отдельные аргументы командной строки.
  • Как отдельная программа. UPD-13.04.2020: я всегда стараюсь этого подхода избегать.
  • MongoDB-инструкция будет такого типа: {'$lookup': {'from': ..., 'let': {...}, 'pipeline': ..., 'as': ...}}.

Разлочить проджекшен src-db-VCF

Сейчас отбирать поля ex-VCF запрещено, т.к. непонятно, как быть с INFO, FORMAT и полями с генотипами. Но есть ведь задачи типа выделения из VCF координат и/или rsIDs и/или аллелей. Так что разлочить надо, но с ворнингом, что сложные объекты тупо отсериализуются, а к вложенным полям доступа нет из-за вредности MongoDB.

Адаптировать к запуску из GUI

Планирую рефакторинг основных компонентов и плагинов, который позволит запускать их с помощью GUI или импортировать в другие Python-программы.

  • Переименовать класс PrepSingleProc в Main, т.к. его __init__ набивается всё большим количеством атрибутов, не связанных напрямую с распараллеливаемой функцией.
  • Перенести алгоритм определения оптимального количества процессов в __init__ вышеупомянутого Main.
  • Всё, что вне Main, сделать выполняемым только при прямом запуске модуля:
if __name__ == '__main__':
  • (UPD-2022.04.25) Вынести в отдельный набор модулей приветственные тексты всех CLI, чтобы их GUI тоже мог задействовать.
  • (UPD-2022.04.26) Принимать пустые значения не только в виде None, но и в виде ''. Первое любит возвращать argparse, второе - streamlit.

Возможность подачи произвольного (до)запроса

query потрясающ, но во многих ситуациях его функциональность избыточна. Нужно, чтобы как можно больше инструментов high-perf-bio, отличных от query, позволяли применить одиночный пользовательский запрос, дополняющий запрос основной. Принимать произвольное MQL/PyMongo-выражение программы должны будут через соответствующий аргумент CLI или виджет Streamlit-GUI.

[count] Ошибка горизонтальной конкатенации элементов, если часть из них - массивы

Удивительно, что я раньше этого не заметил. Если, к примеру, среди вариантов попались мультиаллельники, MongoDB выдаст Unsupported conversion from array to string. Пример куска агрегации, не рассчитанного на списки:
{'$group': {'_id': {'$concat': [{'$toString': '$ID'}, '_', {'$toString': '$REF'}, '_', {'$toString': '$ALT'}]}, <...>}}

GUI (Streamlit)

Задачи:

  • Создать классы, в инитах которых будет накапливаться инпут, полностью совместимый с уже привычным из argparse.
    • Знакомые по CLI приветственные принты реализовать в виде раскрывающихся блоков (st.expander).
    • Все виджеты кроме выбора компонента high-perf-bio завернуть в формы (st.form) с кнопкой запуска (st.form_submit_button). Если тупо, без формы, накидать виджетов, то этап сбора настроек будет проскакиваться.
  • Создать модуль (рабочее название - GUI.py):
    • Импорт GUI-классов, описанных чуть выше.
    • Импорт основных классов всех компонентов с подачей в иниты объектов, мимикрирующих под творения argparse.
    • Выпадающий список с тулзами.
    • Проверка, нажал ли исследователь кнопку запуска.
    • Спиннер (st.spinner), крутящийся, пока задача считается.
  • Вынести тексты, описывающие тулзы, в отдельные модули, и импортировать их как из CLI, так и из GUI (обновлённый meta-issue #15).
  • Ещё касательно #15: в инитах основных модулей принимать проигнорированную исследователем опцию не только в виде None, но и в виде ''.
  • Избавиться от интерактивного диалога пересоздания БД в пользу соответствующего элемента интерфейса (#25).

Проблемы:

  • Невозможно на этапе получения сигнала от виджета проверить, ввёл ли исследователь путь к папке / имя БД (апстрим-фичреквест: streamlit/streamlit#4645).
  • Для запуска веб-интерфейса придётся всё равно иметь дело с командной строкой (streamlit run GUI.py). Я-то потерплю, а каково будет небиоинформатикам? Возможное решение - создать .desktop-ярлык, но это не так уж и просто: #26.

[create] Поддержка SAM

  • Хедер придётся, как и в случае с BED, хардкодить. Он, как минимум, такой:
    QNAME, FLAG, RNAME, POS, MAPQ, CIGAR, RNEXT, PNEXT, TLEN, SEQ, QUAL
  • Теги будут посылаться в подструктуры ещё одного поля. Хоть я и против вмешательства в спеки форматов, придётся выдумать этому полю имя. Допустим, OPTIONAL.

[create_db] Возможность дозаписывать базы

Хочется, конечно, следовать ранее заявленному правилу "новый результат - в новый файл или базу", но практика показала, что фича очень нужная.

Чтобы исследователь осознавал, что модифицирует имеющуюся базу, думаю, стоит сделать отдельный разрешающий флажок.

Сортировка VCF и BED по хромосоме и позиции

Думаю над более конкретной фичей, чем #3.

Для VCF это не просто эстетическая штука, а обязательное требование. Что касается BED, там несортированность не смертельна, но снижает производительность некоторым программам (например, bedtools).

Реализация:

  1. create_db должен создавать компаундный индекс по первым двум (VCF) или трём (BED) полям. Без него парсеры при попытке сортировки по этим полям будут крайне тормозить и даже вылетать. Но от возможности создания раздельных индексов по упомянутым полям отказываться не хотелось бы. Я экспериментально убедился в том, что MongoDB разрешает совмещать несколько индексов для одного поля. Но предстоит ещё разобраться, как такое влияет на производительность запросов. Не запутается ли MongoDB при их выполнении? UPD-11.04.2020: не запутается.
  2. reindex_db. Возможно, будет полезен запрет на удаление обсуждаемого составного индекса.
  3. annotator, intersect_subtract, make_request. sort() для курсора и $sort в aggregation pipeline. UPD-18.04.2020: $sort надо держать в конце пайплайна, чтобы тот не мешал предыдущим звеньям использовать индекс. Сам он при этом индексом пользоваться не сможет, что должно частично компенсироваться режимом внешней сортировки.

Сохранять результаты в gzip-файлы

Многие компоненты принимают сжатые таблицы, но выводят несжатые. Перед тем, как подать вывод одного компонента на вход другому, приходится отдельно запускать башовый однострочник-сжиматель. При частой работе с high-perf-bio эта мелкая недоработка начинает бесить.

[query] Считывать запросы из файлов

Сейчас программа принимает одиночный запрос в качестве аргумента командной строки. В реальных исследованиях приходится для каждого набора данных выполнять много запросов.

Реализация:

  • Считывать запросы из JSON-lines-файлов. Почему "файлов" во множественном числе? На практике запросы то и дело нужно разносить по группам. Пример - получение подпороговых и надпороговых по популяционным частотам вариантов:
{'INFO.0.AF_AFR': {'$lt': Decimal128('0.02')}, 'INFO.0.AF_EUR': {'$gt': Decimal128('0.3')}, 'INFO.0.AF_EAS': {'$gt': Decimal128('0.3')}}
{'INFO.0.AF_AFR': {'$lt': Decimal128('0.02')}, 'INFO.0.AF_EUR': {'$gt': Decimal128('0.3')}, 'INFO.0.AF_EAS': {'$lt': Decimal128('0.02')}}
{'INFO.0.AF_AFR': {'$lt': Decimal128('0.02')}, 'INFO.0.AF_EUR': {'$lt': Decimal128('0.02')}, 'INFO.0.AF_EAS': {'$gt': Decimal128('0.3')}}
{'$or': [{'INFO.0.AF_AFR': {'$gte': Decimal128('0.02')}}, {'INFO.0.AF_EUR': {'$lte': Decimal128('0.3')}}, {'INFO.0.AF_EAS': {'$lte': Decimal128('0.3')}}]}
{'$or': [{'INFO.0.AF_AFR': {'$gte': Decimal128('0.02')}}, {'INFO.0.AF_EUR': {'$lte': Decimal128('0.3')}}, {'INFO.0.AF_EAS': {'$gte': Decimal128('0.02')}}]}
{'$or': [{'INFO.0.AF_AFR': {'$gte': Decimal128('0.02')}}, {'INFO.0.AF_EUR': {'$gte': Decimal128('0.02')}}, {'INFO.0.AF_EAS': {'$lte': Decimal128('0.3')}}]}
  • Продумать структуру папок и файлов, чтобы она не была слишком сложной.
  • Поскольку при наличии индексов запросы выполняются мгновенно, распараллеливание не имеет особого значения. Ради унификации с другими компонентами high-perf-bio, сделаю распараллеливание по исходным (запрососодержащим) файлам.

[left_join] Отбирать поля на уровне MongoDB aggregation

Сейчас поля отбираются блоком Python-кода при потрошении курсора. Если реализовать обычный MongoDB-project {'field_A': 1, 'field_C': 1, field_E: 1}, то пропадут замержденные правые документы, и невозможно будет фильтровать результаты по охвату. Думаю, надо действовать через исключение полей: {'field_B': 0, 'field_D': 0}. Это сложнее, но вполне реализуемо.

[create] Возможность задать произвольный набор имён полей trg-db-TSV

Многие TSV, увы, идут без шапки. А программе она нужна для формирования имён полей в создаваемой базе. Сейчас create подгружает имена полей из самой первой или первой после метаинфы табличной строки. В результате, в случае отсутствия в исходной таблице шапки, ключи документов получатся бессмысленными. Проще прописать шапку на уровне составления команды, чем расковыривать здоровенную таблицу. Кроме того, в ряде случаев, например, при работе с SAM, наличие шапки не приветствуется спецификациями формата. Кстати, о SAM: у меня в далеко идущих планах есть его полная поддержка (#21), но появление возможности добавления кастомного набора ключей было бы неплохим временным хаком, чтобы уже сейчас удавалось загружать SAM в БД хотя бы под видом TSV.

Реализация:
Как в CLI, так и в GUI надо запрашивать у исследователя строку, где элементы соединены с помощью \t. В CLI напоминать о необходимости заключать эту строку в кавычки. Что касается GUI, надеюсь, как можно скорее в @streamlit появится виждет каскадного ввода, и процесс построения шапки значительно упростится. Голосуйте за создание этого виджета!

Переформатировать весь тулкит с помощью autopep8

autopep8 я выбрал, ибо он - наименее агрессивный форматтер. Хочется сохранить собственный, выработанный годами код-стайл, но при этом по-максимуму соблюдать PEP-8.

  • autopep8 иногда довольно неэстетично разносит длинные по его мнению строки кода на несколько строк, поэтому зададим ему заведомо недостижимое значение максимальной длины строки:
"autopep8.args": [
  "--max-line-length=1000"
]

Подробная инструкция для начинающих любителей Visual Studio Code:

  1. CTRL+,
  2. Autopep8: Args
  3. Добавить элемент
  4. --max-line-length=1000
  5. OK
  • Несколько лет назад наличие __pycache__-папок портило результаты ld-tools. С тех пор я запрещаю их создание с помощью sys.dont_write_bytecode = True. Но эта строка должна идти в самом начале модуля, что, естественно, выбесит форматтер. Надо временно притупить его бдительность:
# autopep8: off
import sys; sys.dont_write_bytecode = True
# autopep8: on
  • Просто, чтобы не забыть: запуск переформатирования в vscode - CTRL+SHIFT+I

[create] Дробить генотипы VCF, получая списки

Для более гибкого запрашивания.

Как это будет выглядеть:
'1│0' --> [1, '│', 0]

Обратная совместимость.
Обновлять doc_to_line не придётся, т.к. он изначально способен работать со списками внутри GT.

[download_bio_data] Каждый процесс - скачивание данных одного портала

Сейчас программа параллельно качает несколько файлов из одного портала, потом то же самое - со следующего портала. Это создаёт чрезмерную нагрузку на сервер и вряд ли способствует значимому повышению скорости скачивания. Возможно, именно из-за такого поведения проги FTP Европейского Института Биоинформатики (EMBL-EBI) начинает вредничать, то и дело обрубая загрузку. Хочу переписать программу, чтобы параллелить на уровне порталов, а не файлов.

Аргумент/виджет разрешения пересоздать БД

Сейчас, если программа обнаружила конфликт имён, в диалоговом режиме предлагается стирание БД. Но этот способ слишком привязан к консоли, усложняя разработку GUI (#24). Сделаю для каждого компонента, способного создавать БД, универсальный CLI- и GUI-элемент erase-existing-db. Важные детали:

  • DbAlreadyExistsError оставлю на случай попытки стереть исходную базу, но переселю это исключение в common_errors.
  • В create будет по-прежнему append, но лучше его переименовать в нечто более человеко-понятное, например, replenish-existing-db.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.