Сборка генома de novo
Подготовка чтений
Для работы мне был выдан код доступа SRR4240361 проекта по секвенированию генома бактерии Buchnera aphidicola. По ссылке я скачала архив SRR4240361.fastq.gz с короткими чтениями (длины 36), полученными по технологии Illumina. Перенеся архив в рабочую директорию, я распаковала его с помощью программы:
gunzip SRR4240361.fastq.gz
В результате я получила файл SRR4240361.fastq с чтениями. Затем было необходимо выполнить подготовку чтений программой trimmomatic. Сперва я собрала все адаптеры для Illumina в один файл adapters.fasta, после чего удалила возможные остатки адаптеров из чтений:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240361.fastq SRR4240361_outi.fastq
ILLUMINACLIP:adapters.fasta:2:7:7
В результате имеем:
| Таблица 1. Вывод программы | ||||
| Исходное число чтений Input Reads |
Размер исходного файла | Число оставшихся чтений Surviving (%) |
Размер получившегося файла | Число отброшенных чтений Dropped (%) |
| 7272621 | 733 Мбайт | 7238089 (99.53%) |
729 Мбайт | 34532 (0.47%) |
Затем я удалила плохие буквы с концов чтений (с качеством ниже 20), оставив только чтения длиной не менее 30. Для этого я воспользовалась командой:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240361_outi.fastq SRR4240361_out.fastq
TRAILING:20 MINLEN:30
В результате имеем:
| Таблица 2. Вывод программы | ||||
| Исходное число чтений Input Reads |
Размер исходного файла | Число оставшихся чтений Surviving (%) |
Размер получившегося файла | Число отброшенных чтений Dropped (%) |
| 7238089 | 729 Мбайт | 6881704 (95.08%) |
690 Мбайт | 356385 (4.92%) |
После подготовки чтений я получила файл SRR4240361_out.fastq.
Подготовка k-меров
Программа velveth помогает создать набор данных для работы последующей программы velvetg и сообщает системе, что представляет собой каждый файл с последовательностью. Программа velveth принимает на вход несколько файлов с последовательностями, создаёт хэш-таблицу и на выходе даёт 2 файла (Sequences и Roadmaps) в отдельной директории, которые затем используются программой velvetg. Длина k-мера (hash length), указываемая при запуске программы velveth, соответствует длине хэшируемых слов в парах оснований.
Программа velveth способна обрабатывать файлы форматов: fasta (по умолчанию), fastq, fasta.gz, fastq.gz, sam, bam, raw, raw.gz.
Программа velveth способна обрабатывать чтения следующих категорий: short (по умолчанию, короткие непарные), shortPaired (короткие парные), short2, shortPaired2, long, longPaired, reference.
В нашем случае нужно было подготовить k-меры длины k=29 для коротких непарных чтений (short). Поэтому я воспользовалась программой:
velveth . 29 -fastq -short SRR4240361_out.fastq
В результате программа работала 103.60 секунд, и я получила в рабочей директории файлы Sequences, Roadmaps и Log (отчёт о работе программы).
Сборка на основе k-меров
Программа velvetg строит граф де Брёйна - компактную структуру, основанную на k-мерах. Для этого рассматривается 4-буквенный алфавит (A, T, G, C) и формируется множество всех возможных строк длины k над этим алфавитом (k = 30-1 = 29, так как минимальная длина чтений в нашем случае 30). Затем строится граф B(4,k+1), в котором в качестве вершин стоят все элементы построенного множества, а две вершины v1, v2 соединены направленным ребром тогда и только тогда, когда существует слово длины (k+1), чьим префиксом является строка из v1, а суффиксом строка из v2. Граф содержит n^k вершин. Затем по полученному графу получают контиги. Наиболее простая и популярная стратегия заключается в том, чтобы взять в качестве контигов последовательности, соответствующие отдельным ребрам графа. При этом про структуру графа обычно полностью забывают и больше её не используют. Затем происходит вычисление длин зазоров между контигами, и объединение их в скэффолды.
В командной строке я указала:
velvetg .
В результате я получила файлы: contigs.fa (последовательности контигов в формате fasta), PreGraph (перечень чтений с указанием их порядкового номера в графе), Graph (граф, состоящий из контигов), LastGraph (граф, состоящий из контигов), stats.txt (информация по каждой вершине графа в виде таблицы).
В построенном графе оказалось 1222 вершин, N50 = 49972, максимальная длина контига = 155850, суммарная длина всех контигов = 690757.
С помощью написанной мной программы fasta.py я получила файл contigs.txt, содержащий описание всех контигов (порядковый номер контига|номер вершины|длина|покрытие).
fasta.py contigs.fa contigs.txt
Ниже приведена таблица с характеристиками трёх самых длинных и трёх самых коротких контигов и контигов с максимальным и минимальным покрытиями.
| Таблица 3. Характеристики некоторых контигов | |||
| Порядковый номер | Номер вершины | Длина | Покрытие |
| Самые длинные контиги | |||
| 3 | 3 | 155850 | 33,079514 |
| 11 | 11 | 85024 | 34,670528 |
| 1 | 1 | 72780 | 35,516788 |
| Самые короткие контиги | |||
| 24 | 32 | 29 | 29,827587 |
| 30 | 43 | 29 | 17,310345 |
| 38 | 69 | 29 | 40,137932 |
| Контиги с максимальным покрытием | |||
| 41 | 72 | 36 | 130,750000 |
| 328 | 581 | 32 | 122,593750 |
| 317 | 505 | 38 | 91,473686 |
| Контиги с минимальным покрытием | |||
| 272 | 432 | 58 | 1,724138 |
| 359 | 840 | 29 | 1,827586 |
| 310 | 492 | 59 | 1,847458 |
Самых коротких контигов длины 29 получилось 35 штук из всех 363 контигов. Контигов с покрытием < 3 получилось 45 штук.
Среднее значение покрытия составляет 11,785589.
Медиана покрытий составляет 4,948718.
Мода покрытий составляет 3,000000.
Среднее значение длины составляет 1886,5.
Анализ
С помощью программы megablast я сравнила каждый из трёх самых длинных контигов с хромосомой Buchnera aphidicola (CP009253). Затем повторила запуск программы для двух контигов с аномально большим покрытием. Результаты сравнения приведены в таблице 4. Query - исследуемый контиг, Subject - геном бактерии длиной 628164 пн.
| Таблица 4. Сравнение контигов с хромосомой | |||||||||
| Номер вершины |
Координаты в геноме |
Max score |
Total score |
Query cover |
E value |
Ident | Gaps | Mismatch | Alignment length |
| Самые длинные контиги | |||||||||
| 3 | 353822 - 202390 | 6154 | 63791 | 75% | 0.0 | 91862 77% |
3889 3% |
23803 20% |
119554 |
| 11 | 454069 - 384182 | 3605 | 26639 | 54% | 0.0 | 36579 77% |
1495 3% |
9210 20% |
47284 |
| 1 | 531590 - 462496 | 4047 | 31952 | 78% | 0.0 | 44949 77% |
2022 4% |
11298 19% |
58269 |
Когда я запустила megablast с контигами с максимальным покрытием, то программа выдала "No significant similarity found" (Не найдено значительное сходство). Если посмотреть на вводимые контиги, то они довольно короткие и часто встречаются в геноме. Они сами состоят из повторяющихся фрагментов, из-за чего при сборке их нельзя правильно собрать в контиги большей длины. Поэтому невозможно построить достоверное выравнивание.
>NODE_72_length_36_cov_130.750000 TGATTGTATAATATTATTCGTGGGTACTTGAAACTTCTAAAGTATACTATTATATATCTA TGAT
>NODE_581_length_32_cov_122.593750 AATTTTCGTGATTTTTCCATATTTTGTCATTTTTTGAACTTTAAATGCTTATAAAAAAAA