Сборка de novo
Задание 1: Подготовка чтений
Для этого я использовал программу
$ seqret '*.fa' 'fasta::adapters.fasta'
И вот только после этого я могу воспользоваться программой
$ java -jar /usr/share/java/trimmomatic.jar SE -threads 12 -phred33 SRR4240358.fastq.gz SRR4240358_clear.fastq.gz
ILLUMINACLIP:adapters.fasta:2:7:7 TRAILING:20 MINLEN:32
SE - параметр, поясняющий программе, что у нас чтения одноконцевые (single-ended);-threads 12 - этот параметр говорит программе использовать 12 ядер. Нужно для того, чтобы ускорить процесс;-phred33 - поясняем, какую форму записи качества прочтения нуклеотида имеют наши чтения;SRR4240358.fastq.gz - подаём на вход исходный файл с чтениями;SRR4240358_clear.fastq.gz - имя выходного файла с очищенными чтениями;ILLUMINACLIP:adapters.fasta:2:7:7 - указываем, от каких адаптеров чистить чтения и насколько строго;TRAILING:20 - говорим удалять с правого конца чтений нуклеотиды с качеством прочтения ниже 20;MINLEN:32 - "отставляем в живых" только те чтения, длина которых получилась больше или равна 32 нуклеотидам.
Таким образом получается, что из исходного файла с чтениями было удалено 2527404 последовательности (24%), а размер файла уменьшился на 129 Mb с 470 до 341 Mb.
Задание 2: Подготовка k-меров
Подготавливать k-меры я буду с помощью программы
$ velveth kmers 31 -short -fastq.gz SRR4240358_clear.fastq.gz
velveth - обращение к программеvelveth ;kmers - output-директория, куда будут сложены мои готовые k-меры (созданная предварительно);31 - указываем длину k-меров;-short - указываем программе, что наши чтения короткие и непарные;-fastq.gz - поясняем формат входного файла с чтениями;SRR4240358_clear.fastq.gz - собственно входной файл с чтениями.
Задание 3: Сборка на основе k-меров
Собственно сборку я буду проводить с помощью программы
$ velvetg kmers
velvetg - обращение к программеvelvetg ;kmers - подача на вход программе директории с подготовленными на основе наших чтений k-мерами.
N50 сборки: 8600;
| Таблица 1. Характеристики 3-х самых длинных контигов | |
|---|---|
| Длина, k-меры | Покрытие |
| 19821 | 29.475859 |
| 18714 | 29.922678 |
| 16436 | 30.793624 |
| Таблица 2. Контиги с аномальным покрытием | |
|---|---|
| Длина, k-меры | Покрытие |
| 1 | 111576 |
| 11 | 552.636364 |
| 1 | 499 |
Собственно, насчёт покрытия лично я наблюдаю некую закономерность: чем меньше длина контига, тем больше у него покрытие и наоборот. Это можно объяснить тем, что чем меньше у нас длина нуклеотидной последовательности, тем чаще она будет совпадать со случайной последовательностью той же длины (нуклеотидов всего 4 в наборе, а не 20, как аминокислот в белках), т. е. с уменьшением длины падает уровень уникальности (если можно так сказать) последовательности. Например, самый длинный контиг (см. Таблицу 1) имеет покрытие ~29.5, в то время как контиг с самым большим покрытием (см. Таблицу 2) имеет длину 1 k-мер (по сути никуда не определённый несчастный k-мер). Но при этом стоит отметить, что есть и исключения из моей тенденции: есть в списке контигов индивиды с длиной в районе 50-70 k-меров, которые имеют покрытие 1-2. Это, я думаю, можно объяснить тем, что последовательность, собранная из каких-либо 50-70 k-меров, получилась настолько уникальной, что нашла себе почти однозначное место.
Задание 4: Анализ
В этом задании я буду использовать
| Таблица 3. Характеристики выравниваний контигов и хромосомы ( |
|||||||||
|---|---|---|---|---|---|---|---|---|---|
| № | Длина контига, k-меры | E-value | Кол-во идентичных нуклеотидов, шт. (%) | Кол-во гэпов, шт. (%) | Длина выравнивания, нукл. | Координаты контига | Координаты хромосомы | ||
| Начало | Конец | Начало | Конец | ||||||
| 1 | 19821 | 0.0 | 3256 (75%) | 154 (3%) | 4324 | 948 | 5226 | 496111 | 500325 |
| 6516 (76%) | 351 (4%) | 8617 | 5342 | 13787 | 500370 | 508806 | |||
| 3577 (81%) | 77 (1%) | 4393 | 15478 | 19851 | 510438 | 514772 | |||
| 2 | 18714 | 0.0 | 1977 (78%) | 50 (1%) | 2525 | 1 | 2495 | 8599 | 11103 |
| 2e-110 | 392 (82%) | 9 (1%) | 478 | 5505 | 5979 | 13994 | 14465 | ||
| 0.0 | 2450 (76%) | 86 (2%) | 3225 | 6139 | 9309 | 14727 | 17919 | ||
| 1896 (85%) | 30 (1%) | 2220 | 9387 | 11586 | 17962 | 20171 | |||
| 1509 (82%) | 51 (2%) | 1851 | 12176 | 14000 | 20358 | 22183 | |||
| 2933 (78%) | 140 (3%) | 3779 | 15025 | 18744 | 23067 | 26764 | |||
| 3 | 16436 | 0.0 | 3861 (77%) | 162 (3%) | 5015 | 6919 | 11860 | 467421 | 462496 |
| 5343 (77%) | 204 (2%) | 6961 | 3 | 6889 | 474242 | 467412 | |||
Какие мы можем сделать выводы из выше приведённой информации: контиги легли на хромосому частями и не полностью (покрытия следуют из Таблицы 3: 1 - 87.45%, 2 - 75.23%, 3 - 72.86%), судя по длинам выравнивания, но те куски, что легли, легли прилично, с хорошими показателями выравнивания (см. Таблицу 3). Стоит отметить, что выравнивания разных частей контигов и хромосом не пересеклись, и по координатам хромосомы/контига их можно разместить упорядоченно (в таблице это сделано по координатам хромосомы). Также на всех картах локального сходства мы можем заметить разрывы или участки низкой гомологии (особенно на рис. 2). Достойно внимания и то, что диагональ на карте локального сходства 3-го по длине контига и данной нам хромосомы инвертировалась.
Далее имеет смысл проверить, действительно ли разрывы следуют из того, что там участки низкой гомологии, или же всё-таки это ошибки алгоритма поиска. Для того, чтобы это проверить, я также запустил
На первый взгляд, количество разрывов только увеличилось, но если внимательно сравнить карты локального сходства, то становится понятно, что крупных разрывов стало гораздо меньше, а в целом покрытия выравниваний увеличились (об этом говорят и получившиеся выравнивания: покрытия 93%, 89%, 78% у выравниваний 1-го по длине, 2-го по длине и 3-го по длине контигов с хромосомой соответственно). Это говорит нам о том, что получившееся на предыдущих картах локального сходства (рисунки 1, 2 и 3) разрывы получились в результате неточности алгоритма.