NGS. De-novo сборка генома

В этом практикуме мы попытаемся de-novo собрать геном. Для этого потренируемся на данных секвенирования эндосимбионтной бактерии гороховой тли Buchnera aphidicola. Чтения были взяты здесь. Теперь проверим качество чтений:
fastqc SRR4240380.fastq.gz (отчёт)
Подготовка прочтений к сборке включала две стадии триммирования: обрезку адаптеров и отбор по качеству и длине (в две стадии, чтобы оценить долю адаптеров). Файл adapters.fasta получен объединением файлов адаптеров с kodomo.
java -jar /usr/share/java/trimmomatic.jar SE -threads 20 SRR4240380.fastq.gz trimadapt_SRR4240380.fastq.gz ILLUMINACLIP:adapters.fasta:2:7:7 2> trimadapt.err (лог)
fastqc trimadapt_SRR4240380.fastq.gz (отчёт)
Адаптерами являлось 98 174 (1.88%) последовательностей. Теперь почиститим по качеству с конца с порогом 20 и удалим последовательности короче 32, так как для сборки нужно, чтобы из прочтения можно было выделить хотя бы два 31-мера (слова длиной 31):
java -jar /usr/share/java/trimmomatic.jar SE -threads 20 trimadapt_SRR4240380.fastq.gz trim_SRR4240380.fastq.gz TRAILING:20 MINLEN:32 2> trim.err (лог)
fastqc trim_SRR4240380.fastq.gz (отчёт)
На этом этапе было удалено 253 785 (4.96%) прочтений. До очистки файл весил 108 Мб, а после - 99 Мб (в разархивированном виде было 526, стало 490).
Для сборки используем пакет программ velvet (velveth и velvetg):
velveth velvet/ 31 -fastq.gz -short trim_SRR4240380.fastq.gz
Здесь velvet/ — директория для сборки, 31 — размер k-мера для построения графа, -short — указание на короткие непарные прочтения.
velvetg velvet/
Проанализируем качество сборки. Из лог-файла узнаем N50 = 12 042. Найдем три самых длинных контига и их покрытие. Для этого откроем файл stats.txt в Excel и отсортируем по длине контига (см. таблицу 1)

Таблица 1. Характеристика самых длинных контигов

№	Длина	Покрытие
3	25 915	27,418676
20	23 850	24,763816
23	23 807	25,725921

Интересно, что есть и контиги с покрытием в пять раз большим покрытием (см. таблицу 2). Контиги длиной в несколько k-меров не учитывались, так как, вероятнее всего, это шум.

Таблица 2. Характеристика контигов с самым большим покрытием

№	Длина	Покрытие
56	9342	130,479657
11	2106	126,008547
75	5012	86,361277

Посмотрим, как наши контиги лягут на геном бактерии. Для этого проведем megablast двух последовательностей на геном с номером GenBank/EMBL AC — CP009253.

Таблица 3. Характеристика выравнивания самых длинных контигов на геном

№	Координаты участка хромосомы, соответствующие контигу	Число однонуклеотидных различий	Число гэпов
1	2004 : 11103 613 658 : 620 926 621 055 : 627 104	1992 (22%) 1535 (21%) 1492 (24%)	252 (2%) 190 (2%) 240 (3%)
2	236 859 : 232 358 232 057 : 229 411 252 161 : 248 967	1115 (24%) 529 (20%) 719 (22%)	130 (2%) 71 (2%) 94 (2%)
3	573 092 : 582 686 587 144 : 590 497	2610 (27%) 759 (22%)	471 (4%) 86 (2%)

Посмотрим также на карту локального сходства:

Видно, что контиги легли неплохо, однако контиг 2 выровнялся так, как будто у бактерии произошла инверсия. У контигов 1 и 3 такого нет. Прерывистость объясняется наличием вариабельных участков. Также стоит отметить, что все контиги переходят через нулевую точку, и в этом нет ничего странного, поскольку геном у бактерии кольцевой, а в fasta файле мы вынуждены его записывать в линейном виде, что неизбежно ведёт к разрыву в произвольном месте.

Дополнительные задания

В качестве бонуса, попробуем поиграть с данными и посмотреть, что влияет на качество сборки. В fastqc отчёте изначальных чтений видно, что качество неравномерно падает к концу, а проседает и в середине. Воспользуемся параметром очистки SLIDINGWINDOW, что означает пройтись по каждому чтению скользящим окном длины пять и, если среднее качество в окне ниже 20, удалить эти пять букв и все, что правее:
java -jar /usr/share/java/trimmomatic.jar SE -threads 20 trimadapt_SRR4240380.fastq.gz trim_slid_SRR4240380.fastq.gz SLIDINGWINDOW:5:20 MINLEN:32 2> trim_slid.err (файл)
fastqc trim_slid_SRR4240380.fastq.gz (отчёт)
velveth velvet_slid/ 31 -fastq.gz -short trim_slid_SRR4240380.fastq.gz
velvetg velvet_slid/ (stats)
Из лога видим, что N50 = 12 042, то есть качество сборки не изменилось. Возможно, оно изменилось бы, если бы прочтения имели качество хуже.
Теперь попробуем поменять длину k-мера с 31 на 27:
velveth velvet27/ 27 -fastq.gz -short trim_SRR4240380.fastq.gz
velvetg velvet27/ (stats)
Видим, что N50 = 13 525, то есть немного улучшился. Возможно, перед сборкой нужно подбирать оптимальную величину k-мера, например, с помощью JellyFish.
И, наконец, воспользуемся другим сборщиком. Я решил взять сборщик наших соотечественников - SPAdes. Запустим же его, указав одним из параметров аккуратную коррекцию ошибок:
spades -o spades -t 20 -s trim_SRR4240380.fastq.gz --careful
К сожалению, сам сборщик не выдает статистике о сборке, но его авторы предлагают пользоваться тулом Quast. Последуем их совету. Тул предлагает следующий отчёт. Из него видно, что N50 = 59 008, что в 5 раз выше, чем у Velvet.