Сборка de novo
Мне достался проект SRR4240356 по секвенированию бактерии Buchnera aphidicola. В качестве подготовки необходимых файлов были скачаны сами чтения, а также сохранены в рабочую папку и объединены в один файл остатки адаптеров командами:
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/006/SRR4240356/SRR4240356.fastq.gz
cp ../../adapters/* adapters/
cat *.fa >> adapters_all.fa
Подготовка чтений программой trimmomatic
Теперь, имея все необходимое для дальнейшей работы, стоит начать отбирать только самые лучшие чтения. Для начала были удалены остатки адаптеров командой:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 -trimlog log_trimm.txt -threads 12 SRR4240356.fastq.gz SRR4240356_trimm.fastq ILLUMINA CLIP:adapters_all.fa:2:7:7
В полученном файле осталось 97.96% чтнений. Лог программы тут.
Теперь стоит удалить с правых концов чтений нуклеотиды с качеством ниже 20, оставьте только такие чтения, длина которых не меньше 32 нуклеотидов. Я реализовал это командой:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 -threads 12 SRR4240356_trimm.fastq trimmed_reads.fatsq TRAILING:20 MINLEN:32 2> log_tr imm2.txt
В этот раз сохранились только 95.85% чтений файла без адаптера (305092 чтений было удалено на этом этапе). Файл же стал меньше на 31,72 МБ (777308708Б было и 744045190Б стало). Лог программы тут.
Подготовка k-меров
Теперь пора готовить k-меры для сборки, в данном случае их длина будет 31. Для этого использовалась программа velveth после пары минут разбора параметров ее запуска:
velveth . 31 -short -fastq trimmed_reads.fatsq
Тут "." - команда на запись в текущую директорию, "31" - длина k-мера, "-short" сообщает о коротких непарных чтениях во вносимом файле, "-fastq" показывает формат входного файла а "trimmed_reads.fatsq" - файл с триммированными чтениями. Лог программы тут. Результатом работы послужили три новых файла: Log, Roadmaps, Sequences.
Сборка
Пользуясь полученными данными, можно уже и начать сборку. Для этого использовалась программа velvetg:
velvetg . -cov_cutoff 30
Было выбрано значение "-cov_cutoff" 30 как, по идее, наиболее соответствующее самой идее графов де Брёйна - длина k-мера минус один. Судя по логу выполнения программы, N50 полученных контигов - 65554. Информация о трех самых длинных контигах находится в таблице 1 и получена при помощи команды:
sort -n -r -k 2 stats.txt
| ID контига | 7 | 5 | 9 |
|---|---|---|---|
| Длина контига | 111962 | 107488 | 80939 |
| Покрытие контига | 38.66 | 34.17 | 37.52 |
"Типичное" покрытие, ввиду того, что два контига длины один имеют покрытии по 266951 и 1134, что в разы больше следующих за ними контигов (проанализировано по выдаче
sort -n -r -k 6 stats.txt), стоит считать по медиане:
awk {'print $6'} stats.txt | sort -n | awk '{a[i++]=$1;} END {print a[int(i/2)];}'
Итак, медиана - 40,211. Очевидно, уже описанные контиги (их ID - 33 и 39) длины 1 и с зашкаливающими порытиями очень сильно от нее отличаются. Впрочем, толькоо они - похоже, из-за своей ничтожной длины.
Анализ
Далее каждый из трех самых длинных контигов (ID 7, 5, 9) был выровнен в megablast с уже имеющейся в GenBank хромосомой Buchnera aphidicola при стандартных настройках. Результатом этих выравниваний может служить таблица 2:
| ID контига | 7 | 5 | 9 |
|---|---|---|---|
| Выдача | Тут | Тут | Тут |
| Координаты в хромосоме | 451729-555905 | 220869-323043 | 195400-126623 |
| Координаты в контиге | 2390-109238 | 146-107252 | 37-69033 |
| Идентичность | 81.46% | 78.76% | 74.88% |
| Покрытие | 75% | 74% | 65% |
| Число гэпов | 2092 | 2171 | 1202 |
| Число различий | 16054 | 16294 | 10414 |
| Число выравниваний | 15 | 18 | 11 |
| Dot plot |  |
 |
 |
Сразу стоит отметить, что контиги наложились на известную последовательность далеко не по всей длине, причем не выровнены куски не только с краев контигов - достаточно большие куски несоответствия находятся и внутри последовательности контигов, так что в итоге покрытие во всех трех случаях не превышает трех четвертей длины контига. Кроме того, в отличие от 5 и 7 контигов, 9 контиг наложился на референсный геном в обратной последовательности, что можно объяснить не зависящиим от направления хромосомы секвенированием, в результате которого и были получены исходные риды - часть прямых, часть обратных. В целом такой случай был вполне ожидаем.
Что до довольно обширных зон несоответствия - мне не верится, что они действительно столь велики. Я предположил, что такое их обилие вызвано эвристичностью алгоритма blast, и, чтобы проверить это, повторил выравнивание, понизив длину слова до минимально возможной - 16. И да, на примере 9 контига пробелы между участками заметно сократились (рис. 4), так что, похоже, при выравнивании, например, needle, в выдаче окажутся хоть сколько-то длинные группы гэпов только в начале и конце, но такое выравнивание провести было бы крайне ресурсозатратно. В общем, судя по хорошему выравнивнию, сборка произошла вполне успешно!
