Практикум 15. Сборка de novo

Этот практикум посвящён основным этапам обработки данных при сборке генома de novo - подготовке данных секвенирования, сборке контигов методом построения графа де Брёйна и минимальному анализу результатов.

Подготовка прочтений

Архив с чтениями (полученными при секвенировании генома бактерии Buchnera aphidicola) был взят из базы данных ENA при помощи команды:

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/009/SRR4240379/SRR4240379.fastq.gz

Далее следовало очистить прочтения от возможных остатков последовательностей адаптеров и удалить нечитаемые правые концы. Чтобы отслеживать изменения в количестве и качестве прочтений в ходе очистки, была использована программа fastqc. Результат её применения к исходному файлу с прочтениями (команда ниже) доступен по ссылке. В файле содержится 7 400 155 одиночных прочтений, все имеют длину 36 нуклеотидов, Phred score нуклеотидов ближе к правому концу иногда ниже 20, а в разделе Overrepresented sequences отмечены некоторые адаптеры.

fastqc SRR4240379.fastq.gz

Итак, соберём все предложенные последовательности адаптеров в один файл - adapters.fasta - и применим программу trimmomatic с целью избавиться от них:

java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240379.fastq.gz no_adapters_reads.fastq.gz ILLUMINACLIP:adapters.fasta:2:7:7

В ходе работы программы на консоль была выведена информация о том, что после обработки осталось 7 269 852 прочтений, а 130 303 (1.76% от исходного числа) было отброшено. Результат анализа fastqc можно просмотреть здесь, видно, что последовательности адаптеров исчезли и появилось некоторое количество прочтений малой длины. Ещё раз применим trimmomatic, теперь с параметрами, позволяющими убрать концевые участки прочтений с качеством менее 20 и оставить только прочтения длины не менее 32:

java -jar /usr/share/java/trimmomatic.jar SE -phred33 no_adapters_reads.fastq.gz trimmed_reads.fastq.gz TRAILING:20 MINLEN:32

На этот раз осталось 6 974 267 прочтений, удалено - 295 585 (4.07% от поданных на вход). Анализ fastqc показывает, что, действительно, прочтений длины менее 32 больше нет, а среднее качество концевых участков заметно улучшилось. Размер исходного файла с прочтениями (SRR4240379.fastq.gz) составлял 167 Мб, а файла после фильтрации прочтений - 156 Мб.

Построение k-меров

Для этой задачи предлагалось использовать программу velveth. Она была запущена следующей командой:

velveth assembly 31 -fastq.gz -short trimmed_reads.fastq.gz

Здесь assembly - имя директории для выходных файлов, 31 - длина k-меров, -fastq.gz - опция, указывающая формат входного файла, -short - опция, указывающая, что во входном файле находятся короткие одноконцевые прочтения, trimmed_reads.fastq.gz - имя входного файла. На выходе была получена директория assembly с тремя файлами: Log, Roadmaps, и Sequences.

Сборка

Построение графа де Брёйна на основе полученных k-меров и сборка генома были осуществлены с помощью программы velvetg. Далее приведена команда её запуска:

velvetg assembly

В результате её работы в директории assembly появились файлы: contigs.fa, Graph, LastGraph, PreGraph и stats.txt. Также программа в числе прочего вывела на консоль значение N50 получившейся сборки, которое оказалось равно 25 646. Файл stats.txt содержал любопытную информацию о собранных контигах, поэтому он был импортирован в Excel для дальнейшей работы. Самыми длинными контигами были: ID 6 (длина 49 912, покрытие 35.907237), ID 9 (длина 49 262, покрытие 34.772177), ID 5 (длина 33 085, покрытие 36.25903). Было подсчитано медианное значение покрытия контигов, оно составило 12.703704 и было принято за типичное. Число контигов, покрытие которых отличалось от типичного более чем в 5 раз, оказалось достаточно велико - 103 при общем количестве 440, и отличия были как в большую, так и в меньшую сторону. Самым большим покрытием - 474 299 - обладал контиг ID 133 длины 1 (самым большим, если не считать «вырожденный» контиг ID 292 длины 0 с покрытием infinity), одним из самымх маленьких - 1 - контиг ID 437 длиной 1. Среди контигов неединичной длины наибольшее покрытие - 181.739865, наименьшее - также 1.

Анализ

Для выполнения этого задания три наиболее длинных контига (ID 6, 9, 5) были при помощи megablast выровнены с хромосомой Buchnera aphidicola (GenBank AC: CP009253).

На рисунке 1 представлена карта локального сходства контига 6 и хромосомы Buchnera aphidicola.

Контигу 6 соответствует участок хромосомы 127 825 bp - 173 180 bp. Всего было построено 5 локальных выравниваний, в них - 7 990 несовпадений, 1 027 инделей.

На рисунке 2 представлена карта локального сходства контига 9 и хромосомы Buchnera aphidicola.

Контигу 9 соответствует участок хромосомы 480 874 bp - 529 211 bp. Всего было построено 10 локальных выравниваний, в них - 7 396 несовпадений, 1 073 инделей.

На рисунке 3 представлена карта локального сходства контига 5 и хромосомы Buchnera aphidicola.

Контигу 5 соответствует участок хромосомы 451 729 bp - 480 660 bp. Всего было построено 4 локальных выравнивания, в них - 4 245 несовпадений, 561 инделей.