Сборка de novo

Петренко Павел

Факультет биоинженерии и биоинформатики, Московский Государственный Университет имени М.В.Ломоносова

Введение

В этом практикуме попробуем собрать геном бактерии Buchnera aphidicola str. Tuc7 из коротких ридов длиной 39 нуклеотидов, полученных с помощью Illumina: SRR4240359.

Скачаем архив с чтениями в рабочую директориюс помощью команды wget:

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

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

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

cat /mnt/scratch/NGS/adapters/*.fa > united.fasta

Посчитаем количество отобранных адаптеров (всего получилось 27):

cat /mnt/scratch/NGS/adapters/*.fa > united.fasta

Наконец-то удалим адаптеры из ридов:

TrimmomaticSE -phred33 SRR4240359.fastq.gz clean.fastq.gz ILLUMINACLIP:united.fasta:2:7:7

Увидели, что у нас удалилось 55872 рида, что составляет 0.41% от общего количества:

Input Reads: 13557938 Surviving: 13502066 (99.59%) Dropped: 55872 (0.41%)

Теперь удалим с правых концов чтений нуклеотиды с качеством ниже 20, и оставим только такие чтения, длина которых не меньше 32 нуклеотидов:

cat /mnt/scratch/NGS/adapters/*.fa > delete.fastq.gz

Опять же посмотрим в выдачу и увидим, что теперь у нас отбросилось ещё 1317986 ридов, аж 9.76%:

Input Reads: 13502066 Surviving: 12184080 (90.24%) Dropped: 1317986 (9.76%)

Ради интереса посмотрим вес наших файлов и сравним их:

du -h file

Исходный файл весит 445M.

Файл с обрезанными адаптерами весит 443М.

Файл с ридами лучшего качества (без нуклеотидов низкого качества и коротких ридов) весит 385М.

Видим, что после нашей "чистки" размер файла уменьшился на 60М.

Работа с программой velveth

Velveth - это первая часть конвейера Velvet, программы для сборки геномов из коротких ридов. Воспользуемся velveth, чтобы разбить наши риды на k-меры:

velveth velvet 31 -short -fastq.gz delete.fastq.gz

Опции:
31 - длина k-мера
short - работаем с одноконцевыми чтениями
fastq.gz — указание формата файла, который подаём на вход

В результате в папке velvet появилось 3 файла:
Log - информация о запуске и ошибках
Roadmaps - бинарный файл, который содержит "карты" прохождения ридов через граф де Брёна
Sequences - нужные нам k-меры

Работа с программой velvetg

Velvetg - это вторая часть конвейера Velvet, velvetg берёт "сырой" граф k-mers от velveth, очищает его от ошибок, использует парную информацию для соединения фрагментов и выдаёт готовые контиги/скаффолды. Запускаем:

velvetg velvet

В результате в папку velvet появилось 5 файлов:
contigs.fa - fasta-файл, содержащий последовательности контигов длиной более двух слов из velveth
PreGraph, Graph, LastGraph - три последовательных состояния графа де Брёна в процессе работы Velvetg (каждый следующий файл чище предыдущего)
stats.txt - таблица с разной информацией

Посмотрим, сколько у нас получилось контигов:

grep -c "^>" velvet/contigs.fa

Видим, что всего получилось 285 контигов.

N50 = 70607 (10.3% от общей длины).

Final graph has 709 nodes and n50 of 70607, max 125674, total 682378, using 0/12184080 reads

Отберём 3 самых длинных контига и посчитаем их покрытие:

grep '^>' velvet/contigs.fa | awk -F'_' '{print $1, $2, $3, $4, $5, $6}' | sort -k4 -nr | head -n 3 | less

Находки:
>NODE 11 length 125674 cov 44.550949
>NODE 1 length 108447 cov 42.009186
>NODE 14 length 71403 cov 39.411552

Теперь посмотрим, есть ли у нас аномально большие или маленькие находки. Как я понял, мы берём медианный контиг (получается 143) и смотрим его покрытие, а затем сравниваем с ним все остальные. Найдём медианный контиг:

grep '^>' velvet/contigs.fa | awk -F'_' '{print $1, $2, $3, $4, $5, $6}' | sort -k6 -nr | tail -n143 | head -n1

Получили:

>NODE 342 length 68 cov 5.970588

Видим, что медиана очень маленькая, поэтому, скорее всего, мы найдём только контиги с аномально большим покрытием:

grep '^>' velvet/contigs.fa | awk -F'_' '{print $1, $2, $3, $4, $5, $6}' | awk '$6 > 29.85' | wc -l
grep '^>' velvet/contigs.fa | awk -F'_' '{print $1, $2, $3, $4, $5, $6}' | awk '$6 < 1.92' | wc -l

Как и ожидалось, контигов с аномально маленьким покрытий нет, зато нашлось 52 контига с аномально большим покрытием. Среди таких, например, >NODE_1_length_108447_cov_42.009186 (c покрытием ~ 42.01%, в 7 раз больше медианного) или >NODE_420_length_76_cov_54.263157 (c покрытием ~ 54.26%, в 9 раз больше медианного).

Анализ

Для работы в этом пункте задачи я выбрал рефференсный геном Buchnera aphidicola str. W106 (Myzus persicae) NZ_CP002699.

Контиг 1

Рис. 1. Карта локального сходства для первого контига.
Таблица 1. Информация о выравнивании.
Range 113933-140272
Score 19826
Identities 21463/26651(81%)
Gaps 701/26651(2%)

Контиг почти полностью совпадает с рефференсом.

Контиг 11

Рис. 2. Карта локального сходства для одиннадцатого контига.
Таблица 2. Информация о выравнивании.
Range 628956-643122
Score 12225
Identities 11780/14292(82%)
Gaps 272/14292(1%)

Видим, что на графике есть два хорошо выровнянных участка. Я предполагаю, что такая картина могла образоваться в связи с тем, что у бактерии кольцевая ДНК и начало секвенирования, совпадающее с концом полностью покрылось контигом (сначала захватилась первая область, потом вторая).

Контиг 14

Рис. 3. Карта локального сходства для четырнадцатого контига.
Таблица 3. Информация о выравнивании.
Range 223109-235768
Score 11457
Identities 10634/12789(83%)
Gaps 240/12789(1%)

Если сравним график с таковым для первого контига, то увидим, что он наклонён вниз, а не вверх. Это может говорить нам о том, что выравнивание контига пошло по минус-цепи рефференса.