Сборка генома de novo

В данном практикуме необходимо было работать с проектом по секвенированию бактерии Buchnera aphidicola.

Эмбрион тли (его клетки окрашены синим) содержит крупные специализированные клетки — бактериоциты, в которых живут симбиотические бактерии Buchnera (зеленые).

Бактерия Buchnera aphidicola относится к gamma-Proteobacteria, которые являются грам-отрицательными преимущественно патогенными и азотфиксирующими бактериями. Buchnera aphidicola также являюется первичным эндосимбионтом гороховых тлей Acyrthosiphon pisum. Основной источник пищи Acyrthosiphon pisum - соки растений, поэтому у них естественным образом возникает дефицит незаменимых аминокислот. Решением этой проблемы стал симбиоз c Buchnera, возникший примерно 200-150 миллионов лет назад. В результате такого симбиоза Buchnera aphidicola утратили значительную часть генома и многие ферменты, жизненно важные для свободноживущих бактерий. Геном довольно-таки мал (менее 1 Мб), имеется одна кольцевая хромосома и несколько плазмид.

Мне был дан код доступа SRR4240385.

По ссылке находится архив с короткими (длины 36) ридами, полученными по технологии Illumina.

На странице проекта я скачала fastq - файл в виде архива (.gz) и перенесла его в рабочую директорию (/nfs/srv/databases/ngs/mice.fatmer/pr15), после чего распаковала программой gunzip.

команда: gunzip SRR4240385.fastq.gz

Полученный файл: SRR4240385.fastq

Далее программой FastQC был проведен анализ качества чтений

команда: fastqc SRR4240385.fastq

В результате работы программы был получен архив (.zip) и отчет о результатах работы в виде html страницы [ SRR4240385.fastqc.html ]

Анализ качества ридов дал, мягко говоря, неутешительные результаты. Такое качество прочтений вполне можно назвать непригодным для работы. Отсюда возникает необходимость в дальнейшей очистке ридов.

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

Из файлов в папке с адаптерами (/P/y15/term3/block4/adapters) в рабочей директории (/nfs/srv/databases/ngs/mice.farmer) был создан один файл, содержащий все адаптеры
команда: cat *.fa > /nfs/srv/databases/ngs/mice.farmer/adapters.fasta

В результате получился файл adapters.fasta, который был перенесен в папку pr15 этой директории.

После чего из полученного ранее файла SRR4240385.fastq были удалены все адаптеры команда: java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240385.fastq SRR4240385_nonadapt.fastq ILLUMINACLIP:adapters.fasta:2:7:7

файл: SRR4240385_nonadapt.fastq

Комментарий: ILLUMINACLIP:adapters.fasta:2:7:7 - вырезает адаптеры, со значениями:

2 - отдельные несовпадения, 7 - порог для палиндромной шпильки, 7 - порог для простой шпильки

В результате работы программы было удалено очень малое количество ридов:

Input Reads: 14803535;

Surviving: 14797878 (99,96%)

Dropped: 5657 (0,04%).

Размер файла изменился с 1588 Mb до 1586 Мb.

После этого с концов ридов были удалены нуклеотиды с низким качеством, а также риды длиной менее 30 нуклеотидов

команда: java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240385_nonadapt.fastq SRR4240385_out.fastq TRAILING:20 MINLEN:30

Полученный файл: SRR4240385_out.fastq

Комментарий: MINLEN:30 - удаляет прочтения короче 30; TRAILING:20 - удаляют нуклеотиды ниже качества равного 20 с конца прочтения.

В результате работы программы было удалено значительное (огромное!) количество ридов (что еще раз заставляет усомниться в пригодности данных):

Input Reads: 14797878;

Surviving: 4286467 (28,97%)

Dropped: 10511411 (71,03%).

Размер файла изменился с 1586 Мb до 429 Mb.

Далее программой FastQC был проведен анализ качества чтений после чистки
команда: fastqc SRR4240385_out.fastq

В результате работы программы был получен архивы (.zip) и отчет о результатах работы в виде html страницы [ SRR4240385_out.fastq.html ], после чего было проведено сравнение анализов качества ридов после чистки и изначальных.

Basic statistics

до чистки

после чистки

Качество п.н. в ридах

до чистки

после чистки

Таким образом, после чистки из 14803535 ридов осталось 4286467 (28,95%), т.е. удалено было 10517068 ридов, что, составляет целых 71,05% от исходного числа. При этом качество п.н. повысилось, особенно на концевых нуклеотидах, но нельзя сказать, что изменения существенные, качество по-прежнему остается непригодным, что несложно заметить из графика.

Использование пакета velvet

Сначала были подготовлены k-меры. Подготовка k-меров была произведена с помощью программы velveth. Она предназначена для создания набора данных, которые далее могут обрабатываться программой velvetg.

Velveth принимает на вход несколько последовательностей, строит хэш-таблицу и создает в отдельной директории два файла - Sequences и Roadmaps, необходимые для работы velvetg.

команда: velveth velveth 29 -fastq -short SRR4240385_out.fastq

Комментарии: 29 - (hash length) - длина хэшируемых слов в парах оснований; -short - тип чтений (короткие или длинные, парные или непарные и т.д.); -fastq - формат чтений.

В нашем случае было необходимо подготовить k-меры длины 29 для коротких непарных чтений (-short) из файла в формате (-fastq). Выходные файлы записывались в папку velveth.

Cборка на основе k-меров была произведена программой velvetg с использованием данных, полученных на предыдущем этапе. Velvetg строит граф де Брёйна - ориентированный n-мерный граф из m символов, отражающий пересечения между последовательностями символов. Он имеет m^n вершин, состоящих из всех возможных последовательностей длины n из данных символов. Один и тот же символ может встречаться в последовательности несколько раз. Запуск программы без дополнительных параметров позволят получить fasta-файл с контигами и статистические данные в указанной директории.

команда: velvetg velveth

полученные файлы:

contigs.fa (содержит последовательности контигов)

stats.txt (содержащий статистику)

результат

В построенном программой графе оказалось 764 вершины. Информация по каждой из них отражена в файле stats.txt. Стоит отметить, что количество вершин не обязательно соответствует количеству контигов, так как "нормальными" являются только те контиги, длина которых не менее 29. Именно они прописываются в файле contigs.fa.

N50 составляет 70 п.н. (т.е. ридом с такой длиной и всеми ридами с большей длиной можно покрыть более половины генома).

В приведенной таблице описаны 3 самых длинных контига:

Среди собранных контигов было выбрано 5 контигов с аномальным покрытием (в 5 раз > «среднего»). В следующей таблице приведено их описание:

В принципе, такое большое покрытие неудивительно: длина "контигов" всего 1 k-мер. Не могу однозначно сказать, возможно ли избежать сборки таких контигов изначально.

Анализ

С помощью алгоритма megablast было проведено сравнение каждого из трех самых длинных контигов с хромосомой Buchnera aphidicola ( CP009253). Однако, сравнение не привело к положительному результату: никаких существенных находок обнаружено не было.
Вероятнее всего (наверняка) причиной этому послужило качество ридов.

К семестрам