Сборка генома de novo

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

Ссылка на исходные данные моего проекта SRR4240356. После подготовки файла с адптерами было проведено удаление остатков адаптеров. Отчет о выполнении доступен здесь. В первоначальном файле было 7511529 ридов, после чистки было удалено 153091 (2,04%), осталось 7358438 (97,96%) ридов.
После этого были удалены плохие буквы на концах чтений, а также чтения, с длиной меньше 30. На вход на этом этапе подавалось 7358438 ридов, удалено 73655 (1,00%), осталось 7284783 (99,00%). Исходный файл имел размер 793 999 866 байт, после чистки - 757 103 626 байт. Пройдя по ссылкам, можно увидеть результат выдачи программы fatqc (команда factqc 'имя файла'), для первоначального варианта и конечного. Можно заметить, что действительно были удалены риды, с рамахом неудовлетворительного (ниже 28) качества, а также удалены overrepresented sequences, то есть действительно адаптеры были удалены. Интересен тот файт, что содержание пар G-C составляет примерно 30%, а A-T -70%, то есть исследуемая бактерия не термофильна (поле per base sequence content. Причем видно, что в первоначальном файле содержание колеблется от рида к риду, а в очищенном- равномерное по всем ридам).

Подготовка k-меров и сборка

С помощью программы velveth были подготовлены k-меры нужной длины. Затем была проведена сборка программой velvetg[2], работающий на основе графов де Брайна для k-меров длины 25 и 29. Так как командой sort отсортировать данные не удалось, они были обработаны с помощью excel. Основные параметры приведены в таблице 2.

Таблица 2. Сравнение сборок

Параметр	k-mer length=25	k-mer length=29
N50	5 408	73 133
Общая длина, п.н.	695 924	664 866
Число контигов	3 319	608
Самые длинные контиги	a)ID:51, lgth 14 084, cov 80.41; b)ID:13,lgth 13 372, cov 81.19; c)ID:59,lgth 13 359, cov 76.87.	a)ID:7, lgth 115 468, cov 51.99 ; b)ID:19, lgth 106 076, cov 45.77; c)ID:8, lgth 75 082, cov 54.26.

Видно, что качество сборки для k-меров длины 29 лучше, так как больше N50, в 10 раз длиннее последовательность максимального контига а также в 5 раз меньше общего числа контигов. Поэтому дальше я буду работать только с этим вариантом.
Нашлись контиги с аномальным покрытием. Медиана покрытия составила 10,82. Контиг с максимальным (312 992) покрытием имеет номер 517 и длину 1 нуклеотид, поэтому понятно, почему у него такое большее покрытие. В файле stats.txt содержится информация о всех 610 узлах построенного графа, но в файле contigs.fa содержатся только поулченные контиги, с минимальной длиной, соответствующей длиной k-mer =29 нуклеотидов (таких контигов 16 штук).

Анализ полученных данных

Для анализа собранных чтений были извлечены .fasta файлы самых длинных контигов и сравнены с уже собранной хромосомой Buchnera aphidicola (GenBank/EMBL AC — CP009253) программой blast2seq (megablast). Результаты приведены в таблице 3. Согласно таблице видно, что процент несоответствий (миссматчей) и гэпов в выравниваниях примерно одинаковый между 3 -мя контигами. Расчеты производились в программе excel. Результаты доступны в том же файле, что и для k-меров excel.

Таблица 3. Анализ самых длинных контигов

Ссылка на fasta, длина выравнивания	E-value	Покрытие	Идентичность	Гэпы	Несоответствия
Контиг 7. 86 586 п.н.	0.0	73%	81%	2 205 (2,55%)	16 324 (18,85%)
Контиг 19. 73 984 п.н.	0.0	68%	79%	1 938 (2,62%)	14 770 (19,96%)
Контиг 8. 55 803 п.н.	0.0	73%	83%	1 172 (2,10%)	10 194 (18,27%)

На рисунке 1 представлена визуализация выравниваний для этих контигов. По оси Х расположены контиги , а по У - референсная хромосома. Контиг 7 соответствует участку от 480 660 до 584 329 на хромосоме (17 участков) , контиг 19 - от 248 967 до 349 674 (всего 20 участков. Кажется, что этот участок инвертирован, однако так как программа-сборщик не знает ориентации цепи, то она выбирается случаным образом, поэтому так как обращен контиг целиком, нельзя утверждать о том, что этот участок инвертирован.), а контиг 8 - от 11 103 до 35 124 и от 604 795 до 621 055 (13 участков). Такой большой "разрыв" связан с тем, что бактерии имеют кольцевую хромосому, и просто начало нумерации нуклеотидов в хромосоме совпало с серединой участка, выровненного с контигом, поэтому то мы и видим на карте два участка - один в начале, другой в конце. Анализируя таблицу, можно заметить, что у всех трех контигов процент несоответствий примерно одинаковый, что говорит о том, что контиги собраны правильно, а различия в позиционировании их на геном связаны с наличием делеций и вставок между двумя данными штаммами. Это и неудивительно, ведь референсный штамм BAg симбионт Соевой тли (Aphis glycines), а исследуемый Tuc7 - Гороховой тли (Acyrthosiphon pisum), поэтому их геномы и различаются. Если обратить внимание на общий вид карт, то можно заметить, что последовательности выровнены линейно, есть лишь какие-то делеции и инсерции, то есть общая (коровая) структура генома этих штаммов близка, а различия связаны с адаптациями к разным хозяевам.

Анализ сборки, составленной из половины ридов

Из очищенного файла в отдельный файл была записана половина ридов (wc -l 2.fastq , tail -n 14569564 2.fastq > bad.fastq). Полученный файл был обработан аналогично тому, как описано выше (длина k-меров 29). Примечательно, что программа обрабатывала такой файл 18 секунд, а первоначальный 40-60 секунд. В таблице 4 приведено сравнение качества сборки файла, содержащего все риды и файла, содержащего только половину.

Таблица 4. Сравнение качества сборки

Параметр	Половина ридов	Все риды
N50	30 234	73 133
Общая длина, п.н.	658 385	664 866
Число контигов	339	608
Самые длинные контиги	a)ID11, lgth 42 588, cov 24.60; b)ID:6, lgth 41 368, cov 24.77; c)ID:8, lgth 36 360, cov 25.97.	a)ID:7, lgth 115 468, cov 51.99 ; b)ID:19, lgth 106 076, cov 45.77; c)ID:8, lgth 75 082, cov 54.26.

Видно, что геном собран полностью (общая длина примерно одинаковая), однако N50 в два раза меньше и число контигов тоже в два раза меньше (видимо мало коротких, так как значимых - длиной больше 100 п.н. в изначальном -56, а в урезанном - 64). Кроме того, меньше размер самых длинных контигов (раза в 3), а покрытие в среднем ниже в два раза. То есть можно сделать вывод, что геном собрать можно, однако достоверность будет ниже, но для получения чернового варианта такое количество ридов вполне сгодится.

Сборка с помощью SPAdes

Тот же самый файл 2.fastq был собран с помощью другой программы Spades [3]. Команда для запуска spades.py -s 2.fastq -k 29 -o spades. -s -флаг, указывающий, что исходный файл содержит простые риды, -k -длин k-меров, -o -директория, куда будут записаны выходные файлы. В моем случае этап Mismatch Correction был пропущен, так как риды одиночные, не парные, и нет референсных ридов.
В общем, программа SPAdes содержит гораздо больше опций для работы с разными типами ридов ( по сравнению с velvet), можно задавать файл с достоверными референсными ридами для сравнения на этапе корректировки и проверки ридов, собирать one-cell sequence (--sc), а также продолжать работу с места ошибки. также стоит отметить, что эта программа прездназначена конкретно для сборки бактериальных геномов. Но зато она работает гораздо дольше (вся процедура около 40 мин), чем velvet.
Выходные файлы содержат: contigs.fasta, contigs.fastg – последовательности контигов в разных форматах (fastg - похож на фаста, только содержит куски, записанные в виде графов, что полезно, не нужно записывать много фаста файлов, например для разных аллелей), scaffolds.fasta, scaffolds.fastg – последовательности скаффолдов, corrected/ – директория, содержащая файлы после коррекции ридов (актуально, если риды спаренные), configs/ – рабочие файлы для коррекции ридов, params.txt –параметры SPAdes в этой сессии, spades.log, dataset.info – общая информация о программе SPAdes.
В таблице 5 приведено сравнение работы программы velvet и SPAdes (статистика для программы SPAdes получена с помощью QUAST [4], информация о конкретных контигах получены из названия контигов в fasta-файле, обработано в excel). Интересно, что значимых контигов (длиной больше 100 п.о. в программе velvet получилось 56, а в SPAdes - 111), относительно общего количества гораздо больше, чем полученных программой velvet и они в среднем длиннее.

Таблица 5. Сравнение качества сборки разными программами

Параметр	SPAdes	Velvet
N50	13 841	73 133
Общая длина, п.н.	663 444	664 866
Число контигов	143	608
Самые длинные контиги	ID:1, lgth 32 944, cov 49,08 ; b)ID:3, lgth 31 521, cov 52,40; c)ID:5, lgth 31 329, cov 48,65.	a)ID:7, lgth 115 468, cov 51.99 ; b)ID:19, lgth 106 076, cov 45.77; c)ID:8, lgth 75 082, cov 54.26.

Список источников

[1] Программа Trimmomatic.
[2] Программа Velveth.
[3] Программа Spades.
[4] Программа QUAST.