Сборка de novo

В задании необходимо было поработать с данными проекта по секвенированию бактерии Buchnera aphidicola.
Buchnera aphidicola - вид гамма-протеобактерий, являющихся первичными эндосимбионтами гороховых тлей Acyrthosiphon pisum.
Основной источник пищи Acyrthosiphon pisum - соки растений, поэтому у них естественным образом возникает дефицит незаменимых
аминокислот. Решением этой проблемы стал симбиоз c Buchnera, возникший примерно 200 миллионов лет назад. Buchnera живут в
специальных клетках тли, бактериоцитах, и передаются вертикально, через яйцеклетки матери. Они синтезируют незаменимые
аминокислоты, а в ответ получают от тли множество питательных веществ и среду, богатую азотом. В результате такого симбиоза
Buchnera aphidicola утратили значительную часть генома и многие ферменты, жизненно важные для свободноживущих бактерий [1].
Работать было нужно с короткими (длины 36) чтениями, полученными по технологии Illumina. Мне был выдан код доступа SRR4240356.
На странице проекта я скачала fastq - файл в виде архива .gz и перенесла его в рабочую
директорию (/nfs/srv/databases/ngs/p.avdiunina/pr15), где распаковала программой gunzip.

Использованная команда: gunzip SRR4240356.fastq.gz

Был получен файл с чтениями SRR4240356.fastq.

Очистка чтений

С помощью программы Trimmomatic была проведена очистка чтений, а именно: удаление остатков адаптеров и плохих букв с концов. Для начала все адаптеры для Illumina были собраны в единый файл adapters.fasta. Затем были исполнены следующие команды:

Команда	Назначение	Выходной файл
java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240357.fastq SRR4240357_noad.fastq ILLUMINACLIP:adapters.fasta:2:7:7	Удаление остатков адаптеров	SRR4240357_noad.fastq
java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240357_noad.fastq SRR4240357_trim.fastq TRAILING:20 MINLEN:30	Обрезка с концов чтений нуклеотидов с качеством ниже 20 и отбор чтений длины не менее 30	SRR4240357_trim.fastq

Выдача программы на этапе удаления адаптеров: Input Reads: 7511579 Surviving: 7358438 (97,96%) Dropped: 153091 (2,04%).
Выдача программы на этапе обрезки концов и отбора по длине: Input Reads: 7358438 Surviving: 7075381 (96,15%) Dropped: 283057 (3,85%).

Подготовка k-меров

Подготовка k-меров была произведена с помощью программы velveth. Она предназначена для создания набора данных, которые далее могут обрабатываться программой velvetg. Velveth принимает на вход несколько последовательностей, строит хэш-таблицу и создает в отдельной директории два файла - Sequences и Roadmaps, необходимые для velvetg.
При запуске необходимо указывать длину k-мера (hash length) - длину хэшируемых слов в парах оснований, а также тип чтений (короткие или длинные, парные или непарные и т.д.). Возможные форматы входных файлов - fasta (default), fastq, fasta.gz, fastq.gz, sam, bam, eland, gerald. В нашем случае было необходимо подготовить k-меры длины 29 для коротких непарных чтений (-short) из файла в формате fastq (-fastq). Выходные файлы записывались в папку velveth. Использованная команда: velveth velveth 29 -fastq -short SRR4240357_trim.fastq.

Cборка на основе k-меров

Cборка на основе k-меров была произведена программой velvetg с использованием данных, полученных на предыдущем этапе. Velvetg строит граф де Брёйна - ориентированный n-мерный граф из m символов, отражающий пересечения между последовательностями символов. Он имеет m^n вершин, состоящих из всех возможных последовательностей длины n из данных символов. Один и тот же символ может встречаться в последовательности несколько раз. Запуск программы без дополнительных параметров позволят получить fasta-файл с контигами и статистические данные в указанной директории. Использованная команда: velvetg velveth. В постороенном программой графе оказалось 719 вершин. Информация по каждой из них отражена в файле stats.txt. Стоит отметить, что количество вершин не обязательно соответствует количеству контигов, так как "нормальными" являются только контиги длины не менее 29. Именно они прописываются в файле contigs.fa.

Длины и покрытия самых больших контигов
ID	Длина	Покрытие	Контиги
7	115468	52,223586	contig_7
20	106076	45,974914	contig_20
8	75082	54,512946	contig_8

N50 = 46003 Типичные значения покрытий выделить довольно сложно, так как они достаточно сильно разнятся. Чтобы найти отклоняющиеся показатели, были вычислены наиболее употребимые средние значения, а также найдены максимальные и минимальные покрытия.

Среднее арифметическое - 571.069
Медиана - 10
Максимальные значения покрытий - 599; 595,2; 591
Минимальные значения покрытий - 1; 1,142857; 1,166667

Контиги с аномально большим покрытием
ID	Длина	Покрытие
395	1	365,608
637	1	1660
709	1	735
103	127	664,598425

Видно, что в нашей выборке присутствуют контиги с аномально большим покрытием. Информация по ним представлена в таблице слева. Можно сказать, что для данных контигов характерна сравнительно небольшая длина, что, в общем-то, не удивительно. В случае контигов с минимальным покрытием я не уверена, можно ли их назвать "аномальными". Их покрытия действительно значительно меньше тех показателей, которые хотелось бы назвать типичными, однако в выборке их слишком много (141 контигов с покрытием < 3).

Анализ

С помощью алгоритма megablastn было проведено сравнение каждого из трех самых длинных контигов с хромосомой Buchnera aphidicola (CP009253). Результаты работы megablastn можно увидеть в таблице ниже.

Сравнение самых длинных контигов с хромосомой Buchnera aphidicola
ID	Координаты в геноме	Max score	Total score	Query cover	E-value	Ident	Alignment length	Gaps
7	528794-550219	17304	50759	73%	0.0	81%	21721	545/21721(2%)
20	266073-275551	6137	38814	68%	0.0	79%	9630	363/9661(3%)
8	35124-44693	8521	36733	73%	0.0	83%	9630	125/9630(1%)

Построенные выравнивания можно назвать достаточно неплохими. Далее требовалось выполнить аналогичный анализ для двух контигов с аномально большим покрытием.
При запуске megablastn с дефолтными параметрами ни для одного из 4-х описанных ранее "аномальных" контигов выравнивание построено не было, выдавалось сообщение "No significant similarity found". В FAQ программы указано, что одной из возможных причин такой ошибки могут являться слишком короткие query-sequences, что в нашем случае действительно так. Таким образом, ни для одного из контигов с аномально большим покрытием построить выравнивание не удалось.

Источники: