В качестве исходного файла был дан архив (.gz) .fastq-файла, который необходимо было скачать отсюда. Вообще эти данные представляют собой данные по проекту секвенирования генома бактерии Buchnera aphidicola.
Эта бактерия является симбионтом тлей рода Acyrthosiphon, синтезируя незаменимые аминокислоты и получая питательные вещества и среду, обогащенную азотом.
После скачивания архива, он был перемещен в папку /nfs/srv/databases/ngs/emadevich/pr15 и распакован gunzip.

Команда: gunzip SRR4240357.fastq.gz. На выходе получаем .fastq-файл SRR4240357.fastq.

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

Следующим шагом стала непосредственная очистка чтений (удаление адаптеров (1) и букв плохого качества с концов чтений (2)) программой trimmomatic. Для этого был создан файл, содержащий адаптеры. Следующие команды были использованы:

Команда для (1): java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240357.fastq SRR4240357_wad.fastq ILLUMINACLIP:adapters.fasta:2:7:7
На выходе: 7937705 чтений записаны в SRR4240357_wad.fastq, 161274 чтений было удалено.
Изменение размера файла: размер файла уменьшился с 863 Мбайт до 845 Мбайт.

Команда для (2): java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240357_wad.fastq SRR4240357_trim.fastq TRAILING:20 MINLEN:30
На выходе: 7248100 чтений записаны в SRR4240357_trim.fastq, 689605 чтений были удалены.
Изменение размера файла: размер файла уменьшился с 845 Мбайт до 748 Мбайт.

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

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

Velveth принимает на вход файл с ридами, строит хэш-таблицу и создает в отдельной директории два файла - Sequences и Roadmaps, необходимые для velvetg и один, Log, в котором содержится справочная информация. При запуске необходимо указывать hash length - длину хэшируемых слов (bp), а также тип чтений.
В нашем случае было необходимо подготовить k-меры длины 29 для коротких непарных чтений (-short) из файла в формате fastq (-fastq).

Команда: velveth velveth 29 -fastq -short SRR4240357_trim.fastq.

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

После подготовки всех необходимых файлов была произведена сборка. Программа использует алгоритм, строящий граф де Брёйна - граф, отражающий пересечения всех возможных последовательностей длины n.

Команда: velvetg velveth.

На выходе несколько файлов, из которых нам будут нужны два: .txt и .fa. Первый содержит данные о всех возможных контигах (даже те, которые не удовлетворяют нашим условиям: длина не менее 29 bp). Второй - fasta последовательности контигов.
Из первого файла были удалены все контиги с длиной меньше 29 bp, в итоге из 914 контигов осталось 178. В этой таблице приведены основные данные по оставшимся контигам (среднее покрытий - 81.03, медиана покрытий - 14.57, N50 длины - 28739). Ниже составлена таблица для трех самых длинных контигов:


К слову, контигов, отличающихся от среднего значения/медианы покрытия, довольно много. Ниже приведена сводная таблица по трем контигам, покрытия которых очень велики. Такой же таблицы по контигам с небольшим покрытием (1-3) к сожалению привести не представляется возможным, так как количество таких контигов 24.


Видно, что приведенные контиги имеют малую длину, что достаточно логично (например максимальная длина контига при покрытии не меньше 100 - 212).

Анализ

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


Видно, что все построенные выравнивания достаточно неплохие.

Что же касается выравнивания контигов с аномально большим покрытием, то тут возникла проблема. Ни при каких изменениях дефолтных параметров (и при самих дефолтных параметрах) ни одного выравнивания потсроено не было. На каждом выходе алгоритм выдавал следующее: 'No significant similarity found'. К слову ни одна из приведенных ниже tips не помогла:


Сборка генома при k = 25

Были проделаны те же шаги и получена эта таблица, что и для случая с k = 29. Как и ожидалось, количество контигов возросло: 1025 при k = 25 против 178 при k = 29 (дальше случай с k = 29 будет писаться в скобках). Что касается основных статистических показателей, то тут тоже все ожидаемо: N50 длины - 2067 (28739), среднее покрытий - 36.41 (81.03), медиана покрытий - 47.1 (14.57). Интересно, что в случае k = 29 большой процент контигов имеют малое покрытие, в отличие от таковых в при k = 25 (видно из значений медианы относительно среднего).
Ниже приведена таблица для трех контигов, обладающих самой большой длиной:


Видно, что по сравнению с длинами контигов при k = 29, максимальная длина контигов при k = 25, куда меньше (примерно в 5 раз), что логично, ведь увеличилось количество самих контигов (примерно в 9 раз!), а длина сборки осталась приблизительно такой же.

Сборка при помощи SPAdes

SPAdes является еще одним сборщиком геномов, выдачу которого необходимо было сравнить с выдачей velvet.
Изучив мануал, была использована следующая команда:
spades.py -k 29 --only-assembler -s [file.fastq] -o [folder name], где флажок -k означает длину k-мера, --only-assembler - включен только режим сборки, -s - обработка одноконцевых чтений, которые лежат в файле [file.fastq], -o - запись в папку [folder name]

В итоге программа работала 100 минут (!). Несмотря на такую длительную работу, выдача программы была куда лучше, чем выдача velvet.
Так, количество контигов - 46, N50 длины - 77869, среднее покрытий - 9093.68, медиана покрытий - 30.1. Видно, что значение медианы кардинально сильно отличается от значения среднего, что говорит о наличии нескольких (в данном случае двух) контигов с аномально большим покрытием. Ниже приведена таблица с описанием нескольких контигов с аномально большим покрытием:


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

Сборка генома при отсутствии половины контигов

Для того, чтобы удалить половину чтений, можно было использовать команду head:
head SRR4240357_trim.fastq -n 14496200 > ilov.fastq, где флажок -n означает количество строчек, которые нужно записать

В итоге, после использования программы velvet (ну просто потому что она работает быстрее, чем SPAdes) и использования Excel, была получена таблица. Следующие статистические характеристики были посчитаны: количество контигов - 169, N50 длины - 10700 (против 28739), среднее покрытий - 26.66 (против 81.03) и медиана покрытий - 17.34 (против 14.57).
Видно, что значения медианы и среднего почти одинаковы, что говорит о том, что значения покрытий распределены более-менее равномерно, то есть нет аномальных значений (на самом деле есть, но их буквально 1-2 штуки).

Что касается значений длин, то тут наблюдается похожая ситуация, что и с покрытием: значения длин распределены достаточно равномерно. Еще одним подтверждением являются значения самых длинных контигов:


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

Ниже приведена таблица выдачи алгоритма megablast для трех самых длинных контигов:


Сразу бросается в глаза то, что значения покрытия (query cover) стали меньше, при этом остальные параметры не особо изменились.

Если говорить в целом, то стоит сказать, что такие результаты в ДАННОМ конкретном случае объяснимы. В среднем покрытие в начальном случае достаточно большое, поэтому даже удаление половины ридов не несёт катастрофических изменений (хотя понятно, что будь покрытие в начальном случае от 1 до 10, то после удаления половины чтений, результаты были бы плачевными).
Резюмируя, можно лишь сказать, что в выдаче стоит ожидать худшего, а как оно получится на самом деле, зависит от данных.