Практикум 15. De novo сборка генома

Основное задание

Задачей данного практикума является de novo сборка генома бактерии Buchnera aphidicola из чтений, имеющих ID: SRR4240361.

Чтения были скачаны из базы EMBL коммандой :

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/001/SRR4240361/SRR4240361.fastq.gz

Прежде чем начать сборку, чтения необходимо триммировать - для этого был использован trimmomatic с опицей -SE (опция для работы с одноконцевыми чтениями) .Сначала я избавился от адаптеров командной :

java -jar /usr/share/java/trimmomatic.jar SE -phred33 -threads 15 SRR4240361.fastq.gz SRR4240361_no_adapters.fastq.gz ILLUMINACLIP:adapters.fasta:2:7:7

В результате из 7272621 ридов осталось 7238089 , т.е. было выброшено 0.47% ридов. Затем я удалил с правых концов чтений нуклеотиды с качеством ниже 20 и удалил все риды длинной ниже 32. Команда :

java -jar /usr/share/java/trimmomatic.jar SE -phred33 -threads 15 SRR4240361_no_adapters.fastq.gz SRR4240361_good.fastq.gz TRAILING:20 MINLEN:32

В результате из 7238089 ридов осталось 6834335, т.е. было потеряно 5.58% ридов.

Далее из тримированных чтений были получены k-меры длинной 31, которые были использованы для сборки генома. Команды :

velveth ./SRR4240361_velvet 31 -fastq.gz -short SRR4240361_good.fastq.gz

velvetg ./SRR4240361_velvet

В результате сборки было получено 474 контигов. N50 сборки составил 25713 bp (то, что программа выдает длины в количестве k-меров учтено). Аномальных покрытий не наблюдается. Самые длинные контиги представлены в Таб.1 .

Таб.1Самые длинные контиги
Контиг Длина контига Покрытие
NODE_6_length_49238_cov_26.660851 49238 26.66
NODE_2_length_45555_cov_26.450466 45555 26.45
NODE_2_length_43866_cov_23.514977 45555 23.514977

Затем последовательности этих контигов были выравнены на геном нашего организма с помощью megablast (wordsize= 28) и blastn (wordzise=7). ID генома: CP009253

Письма мастера дзен Рис1. Dot plot полученный megablast выравниванием генома с NODE_6_length_49238_cov_26.660851
Письма мастера дзен Рис2.Dot plot полученный blastn выравниванием генома с NODE_6_length_49238_cov_26.660851

На Рис1 и Рис2 показаны Dot Plot выравнивания NODE_6_length_49238_cov_26.660851 с геномом алгоритмами megablast и blastn. На dot plot полученном megablast последовательности довольно хорошо выравнены. Имеется лишь немного пустых участков. Но особо крупных инверсий и делеций не наблюдается. Если посмотреть на карту, полученную blastn, то некоторые из этих пустых участков выравнивались, что говорит о том, что там есть мисмэтчи.

Информацию о выравнивании контига можно увидеть в Hit Table файле, так как у меня получилось многовато(45) выравниваний и делать из них таблицу затруднительно.

Письма мастера дзен Рис3. Dot plot полученный megablast выравниванием генома с NODE_2_length_45555_cov_26.450466
Письма мастера дзен Рис4.Dot plot полученный blastn выравниванием генома с NODE_2_length_45555_cov_26.450466

Этот контиг тоже довольно хорошо выравнился. Аналогично предыдущей карте локальных сходств имеюся пустые участки, которые выравниваются в blastn. Здесь линия выравнивания идет сверху вниз, что говорит о выравнивания с комплементарной цепью. Либо же этот участок подвергся инверсии. Учитывая длину контига такой вариант вполне возможен. Если смотреть на выдачу blastn, то есть участок между 26k-28k по оси абсцисс, который не выравнялся. Пристально посмотрев у меня появилось ощущение, что там есть сдвиг графика, т.е произошла делеция. Если посмотреть на Hit Table данные о выравнивании и на карту локального сходства, то можно увидеть, что выравнивание заканчивается на координатах [25826,460011], а продолжается с [28447,458588]. 28447-25826=2621. Получается относительно оси ординат график должен упасть на 2621, прежде чем заново начать выравниваться.Но у меня выходит |458588-460011|=|1423|.Получается, что случился сдвиг, которые, возможно, является делецией. (я надеюсь это так работает)

Письма мастера дзен Рис5. Dot plot полученный megablast выравниванием генома с NODE_34_length_43866_cov_23.514977
Письма мастера дзен Рис6.Dot plot полученный blastn выравниванием генома с NODE_34_length_43866_cov_23.514977

Файл с информацией о выравнивании. Контиг очень хорошо выравнился с участком генома, хоть и имеются некоторые пустые участки. Цепи были выбрано одинаково. При просмотре карты выданной blastn многие пустые участки выравнились, что говорит о наличии там мисмачтей.

Дополнительное заданин. k-меры 27

В рамках первого доп.задания я из тримированных чтений получил k-меры длинной 27 и использовал их для сборки. В результате сборки N50 оказалось ниже - 15564. Теперь самыми длинными контигами являются :

Таб.2Самые длинные контиги (k-меры длинной 27)
Контиг Длина контига Покрытие
NODE_14_length_50427_cov_46.11 50427 46.114819
NODE_4_length_39324_cov_41.74 39324 41.749313
NODE_1_length_34494_cov_45.15 34494 45.153882

Можно заметить, что в сравнении с прошлой сборкой эти контиги короче (кроме самого длинного), а покрытие у них больше. В целом, контигов стало больше. Если при прошлой сборке их было 474, то сейчас 2531. Можно сделать вывод, что при понижении длины k-меров сборка стала хуже.

Дополнительное заданин. Другие параметры триммирования

В рамках данного задания я пробовал иначе триммировать свои чтений (уже избавленные от адаптеров). Я использовал параметр SLIDINGWINDOW:5:20, который проходит по каждому чтению окном длины 5 нуклеотидов, и если среднее качество какого-либо участка чтения оказывается ниже 20, то удаляются эти 5 букв и все что, правее. Команда:

java -jar /usr/share/java/trimmomatic.jar SE -phred33 -threads 15 SRR4240361_no_adapters.fastq.gz SRR4240361_good.fastq.gz SLIDINGWINDOW:5:20

В результате запуска из 7238089 чтений осталось 7065728, т.е. было потеряно 2.38%.

Теперь из этих чтений я получаю k-меры длинной 31 и собираю геном. В результате сборки N50 получился 25713 (точно такой же, как и при первой сборке).В целом результаты сборки очень схожи с первой сборкой. Общее число образовавшихся контигов - 406. Двойка самых длинных контигов имеют точно такую жу длину, как и при первой сборке. Лишь немного изменилось их покрытие. (Таб.3)

Таб.3Самые длинные контиги (SLIDINGWINDOW:5:20)
Контиг Длина контига Покрытие
NODE_6_length_49238_cov_26.203481 49238 26.20
NODE_2_length_45555_cov_26.018329 45555 26.01
NODE_32_length_43866_cov_23.069667 43866 23.06

В целом можно сделать, вывод что качество сборки не сильно изменилось при использовании другого пути триммирования чтений.

Дополнительное задание. Сборка программой SPAdes

В рамках этого задания для сборки генома я пользуюсь прораммой SPAdes. При этом для сборки я использую k-меры длинной 31 (как и в основном задании) и соответсвенно сравнивать я буду сборку из основного задания и ту, которую сделаю здесь. Команда для сборки:

spades.py -s SRR4240361_good.fastq.gz -k 31 --only-assembler -o SRR4240361_SPAdes

В таблице 4 показаны некоторые характеристики сборок, как число образовавшихся контигов, N50 (учтено, что теперь программа выдает длину в количестве нуклеотидов), длины самых длинных контигов и их покрытие. В результате сборки программой SPAdes получается значительно меньше контигов и более высокое значение N50, чем в результате сборкой программой velvet. Из этого я делаю вывод, что сборка полученная SPAdes превосходит сборку полученную velvet.

Таб.4Сравнение сборок программами SPAdes и velvet
Программа N50 Длина самого длинного контига Покрытие самого длинного контига Число контигов
Velvet 25713 49238 26.66 474
SPAdes 30583 49269 21.87 122