FASTQ->BAM

[Val_Metov]

Собрал геном в HG38 3 разными способами. BWA mem. Minimap2 и Strobealign. Последние 2 заметно быстрее в работе. Minimap2 хорошо справляется с длинными ридами. И на перспективу больше подходит.
Прогнал все 3 файла через deepvariant. VCF Y only все три объеденил. Добавил информацию по снипам по данным Yfull. Никаких особо критичных отличий нет, но bwa и minimap2 кажутся более надежными.
Ссылка на excel файл сведенный.
https://disk.yandex.ru/i/QQM4ol9MFCNtiA

[Val_Metov]

Подготовил скрипты для массовой обработки fastq, cram.
Желательно мощная система и достаточный объем Ram. На 64 гб + 32 потока все работает хорошо.
https://docs.google.com/document/d/1rgqlOd6JJXeK0nNqeSGFf-jKafqdY_R7G3Vmn9o9KX0/edit?usp=sharing

[Semargl]

Подготовил скрипты для массовой обработки fastq, cram.
Желательно мощная система и достаточный объем Ram. На 64 гб + 32 потока все работает хорошо.
https://docs.google.com/document/d/1rgqlOd6JJXeK0nNqeSGFf-jKafqdY_R7G3Vmn9o9KX0/edit?usp=sharing

Плюс.
Несколько вопросов.
"minimap2 -ax sr -t $(nproc) -I $memory -B $batch_size....."
Если настройки статичны, то почему не сохранить индекс референса один раз, а потом его использовать многократно? Конечно если не меняем опцию "-I"...
"-B" разве это не относится к установке пенальти при выравнивании и этот параметр по умолчанию 4 вроде.

-B INT       mismatch penalty (larger value for lower divergence) [4]

https://lh3.github.io/minimap2/minimap2.html

"samtools view -b $output chrY chrM > "${results_folder}/$(basename ${genome})_YM.bam"" и
"samtools view -T $ref_genome -C -o "${results_folder}/$(basename ${genome}).cram" $output'
можно бы добавить "-@ $(nproc)"

и в "@RG\\tPL:ILLUMINA\\tSM:Nebula" добавить бы уникальный ID. По типу "@RG\tID:${3}\tSM:${3}\tLB:library1"

Мне кажется проблема может быть только с "-B", в остальном это мелочи. Тоже пробую минимап для перемаппирования сотен научных. Важна скорость.

[Semargl]

Еще можно было бы улучшить конвейер, добавив fixmate и markdup. Конечно это увеличит время, но имхо так кошернее. Вопрос правда насколько это улучшение стоит потери производительности?)
Примерный алгоритм:
minimap2 | samtools fixmate | samtools sort | samtools markdup

[Semargl]

По перегону BAM/CRAM -> FASTQ.
Ранее пользовался пайпом samtools bamshuf | samtools bam2fq
но конструкция устарела и была выпилена из samtools
Теперь работает так:
samtools collate -u -O -@ 80 HG01494.final.cram col_HG01494 | samtools fastq -@ 80 -1 HG01494.1.fastq.gz -2 HG01494.2.fastq.gz -n

Вроде авторы samtools рекомендуют использовать collate вместо view

[Val_Metov]

Несколько вопросов.
"minimap2 -ax sr -t $(nproc) -I $memory -B $batch_size....."
Если настройки статичны, то почему не сохранить индекс референса один раз, а потом его использовать многократно? Конечно если не меняем опцию "-I"...

Для однотипных референсов я так и делал, но когда в дело вступает YSEQ там целый зоопарк разных версий начиная со 100 ридов и до 150. Плюс генерация референса занимает пару минут и я решил каждый раз создавать. В масштабах сборки всего генома это особо не влияет.

"samtools view -b $output chrY chrM > "${results_folder}/$(basename ${genome})_YM.bam"" и
"samtools view -T $ref_genome -C -o "${results_folder}/$(basename ${genome}).cram" $output'
можно бы добавить "-@ $(nproc)"

и в "@RG\\tPL:ILLUMINA\\tSM:Nebula" добавить бы уникальный ID. По типу "@RG\tID:${3}\tSM:${3}\tLB:library1"

Мне кажется проблема может быть только с "-B", в остальном это мелочи. Тоже пробую минимап для перемаппирования сотен научных. Важна скорость.

Я с этими RG так и не разобрался. Хотел с fastq забирать часть заголовка и туда внедрять. Но в итоге забил.
Насчет научников, могу помочь тоже. Если скрипт поднастроить под эти конкретные задачи. Мои ресурсы сейчас, пока свежих геномов у меня нет, а свои ~ 40 я уже пересобрал, в почти всегда простаивают. Штук 15 ПГ в день можно обрабатывать. 2+2 tb m2 SSD должно хватить.

Там вроде еще более ускоренный вариант minamap2 есть. Но у меня пока с ним не особо заладилось, да и 1 час на сборку с обычным Minimap2 вполне себе хороший срок.
Спасибо за подсказки, сравним скорость и результаты.

Еще можно было бы улучшить конвейер, добавив fixmate и markdup. Конечно это увеличит время, но имхо так кошернее. Вопрос правда насколько это улучшение стоит потери производительности?)
Примерный алгоритм:
minimap2 | samtools fixmate | samtools sort | samtools markdup

Я немного с Рэнди общаюсь, одним с авторов WGSExtract. Они могут туда добавить скрипты, если какие-то рекомендации будут, к T2T итд. Он что-то не очень хвалил тоже YSEQ методы сборки.

[Semargl]

Для однотипных референсов я так и делал, но когда в дело вступает YSEQ там целый зоопарк разных версий начиная со 100 ридов и до 150. Плюс генерация референса занимает пару минут и я решил каждый раз создавать. В масштабах сборки всего генома это особо не влияет.

Согласен что пару минут не утянет, и можно генерить индекс каждый раз заново, с учетом возможного изменения входного параметра "-I".
Но при перегоне CRAM/BAM реф не нужен, если пользоваться рекомендованным алгоритмом (collate вместо view):
samtools collate -u -O -@ 80 HG01494.cram col_HG01494 | samtools fastq -@ 80 -1 HG01494.1.fastq.gz -2 HG01494.2.fastq.gz -s HG01494.S.fastq.gz -n
Конечно если не сильно влияет на скорость. Не уверен почему Рэнди использует при перегоне именно view и не знаю будет ли результат аналогичен использованию collate. Возможно у Рэнди свои соображения.

Я с этими RG так и не разобрался. Хотел с fastq забирать часть заголовка и туда внедрять. Но в итоге забил.

Использование @RG является хорошим тоном и используется для автоматического получения уникального имени семпла при автоматической обработке, а так же при формировании сравнительных файлов VCF значение из ID используется для обозначения столбцов. Кстати, Томас игнорирует заполнение этого поля.

"-B" разве это не относится к установке пенальти при выравнивании и этот параметр по умолчанию 4 вроде.

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

Насчет научников, могу помочь тоже. Если скрипт поднастроить под эти конкретные задачи.

Спасибо, напишу в личку. Скрипты накидал. Опубликую здесь, возможно кому пригодятся. Грубые, без циклов, но суть понятна. Команду перегонки CRAM/BAM->FASTQ привел выше в этом сообщении.

Он что-то не очень хвалил тоже YSEQ методы сборки.

Возможно из-за того что выравнивает без использования @RG, fixmate и markdup?

[Semargl]

Выравнивание

Код: [Выделить]

#!/bin/bash
reference=/path/to/file/chm13_v2.0.fna.gz
minimap2=/path/to/file/minimap2
samtools=/path/to/bin/samtools

set -o pipefail
$minimap2 -ax sr -t $4 -I 60G -a -R "@RG\tID:${3}\tSM:${3}\tLB:library1" $reference $1 $2 | \
$samtools fixmate -@$4 -O bam,level=0 -m - - | \
$samtools sort -O bam,level=0 -@$4 -T tmp - | \
$samtools markdup -@$4 --write-index --reference $reference - $3.dedup.bam && \
$samtools view -b $3.dedup.bam chrY chrM > $3.YM.bam && \
echo "${3} done!"

При желании можно закомментировать строки с fixmate и markdup
Вызывать с параметрами:
./my_script.sh FASTQ.R1.gz FASTQ.R2.gz IDNAME 80
Где IDNAME - то что будет прописано в @RG и используется для названия выходного файла. Только латиница и цифры. Короткое. Уникальное.
80 - количество потоков

[Val_Metov]

С учетом рекомендаций Владимира модифицировал скрипт. Можно и для bwa-mem использовать. Нужно только раздел команд чуть поменять. Основная часть кода занимает поиск r1 и r2. Название файла берется из r1 и также добавляется в RG. + добавлена генерация лог файла. Хотел лог для каждой операции добавить, но корректно не получилось.

[Val_Metov]

https://docs.nvidia.com/clara/parabricks/latest/documentation/tooldocs/man_deepvariant_germline.html
Clara Parabricks v4.2.0
Работает на GPU Nvidia >= 16 gb gpu ram.
У меня к сожалению только 11 гб, проверить не получилось. out of memory

[Semargl]

https://docs.nvidia.com/clara/parabricks/latest/documentation/tooldocs/man_deepvariant_germline.html
Clara Parabricks v4.2.0
Работает на GPU Nvidia >= 16 gb gpu ram.
У меня к сожалению только 11 гб, проверить не получилось. out of memory

Странно. Может есть параметры, позволяющие ограничить требования к ресурсам?
Может это?

--bwa-cpu-thread-pool BWA_CPU_THREAD_POOL
Number of threads to devote to CPU thread pool per GPU (Argument only applies to --fq2bamfast) (default: 16)

К примеру попробовать 8 потоков. Ну или другие опции запуска.

Я не нашел в документации минимальных требований по железу.

[NathanS]

Я не нашел в документации минимальных требований по железу.

https://docs.nvidia.com/clara/parabricks/latest/gettingstarted.html#hardware-requirements

Hardware Requirements
Any NVIDIA GPU that supports CUDA architecture 70, 75, 80, 86, 89 or 90 and has at least 16GB of GPU RAM. NVIDIA Parabricks has been tested on the following NVIDIA GPUs: V100....

The fq2bam tool requires at least 24 GB of GPU memory by default; the --low-memory option will reduce this to 16 GB of GPU memory at the cost of slower processing. All other tools require at least 16 GB of GPU memory per GPU.

System Requirements:

A 2 GPU system should have at least 100GB CPU RAM and at least 24 CPU threads.

A 4 GPU system should have at least 196GB CPU RAM and at least 32 CPU threads.

A 8 GPU system should have at least 392GB CPU RAM and at least 48 CPU threads.

[Val_Metov]

Tesla p40 24 гб, аналог 1080ti. На бу рынке от 15 тыс. Надо еще охлаждение приколхозить. Интересно тренировкой нейросетей заняться, ищу по разумной цене 3090
Minimap2 и так довольно быстро собирает, но интересно попробовать разные технологии.

[NathanS]

Tesla p40 24 гб, аналог 1080ti.

Tesla P40 имеет устаревшую архитектуру Pascal SM60-61, а по требованиям там от SM70 (Volta: Tesla V100, GTX 1180 (GV104), Titan V, Quadro GV100). 3090 уже подходит, или SM89 or SM_89, compute_89 – NVIDIA GeForce RTX 4090, RTX 4080, RTX 6000, Tesla L40

Слишком дороговато для экспериментов, да и с требованиями по количеству памяти и процессорам достаточно будет вложится в хороший процессор. У меня выравнивание занимает около 1 суток, а если со стандартными параметрами, то 8-10 часов. Так как я не выравниваю дестками геномов в сутки, то в принципе приемлемо.

[Val_Metov]

Tesla p40 24 гб, аналог 1080ti.

Tesla P40 имеет устаревшую архитектуру Pascal SM60-61, а по требованиям там от SM70 (Volta: Tesla V100, GTX 1180 (GV104), Titan V, Quadro GV100). 3090 уже подходит, или SM89 or SM_89, compute_89 – NVIDIA GeForce RTX 4090, RTX 4080, RTX 6000, Tesla L40

Слишком дороговато для экспериментов, да и с требованиями по количеству памяти и процессорам достаточно будет вложится в хороший процессор. У меня выравнивание занимает около 1 суток, а если со стандартными параметрами, то 8-10 часов. Так как я не выравниваю дестками геномов в сутки, то в принципе приемлемо.

Железо у меня мощное. Выравнивание через minimap2 занимает 1-2 часа, в зависимости от доп обработок. В принципе этого более чем достаточно. Но интересно попробовать все что можно попробовать.

Насчет p40 спасибо за дополнения.
Я правда стоимость железа перекладываю на стоимость ПГ и этот аргумент обычно побеждает.