На сайте ENA нашли проект по секвенированию бактерии Buchnera aphidicola, на нём копировали ссылку для скачивания SRR4240359.fastq.gz. Командой wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/009/SRR4240359/SRR4240359.fastq.gz скопировали файл в рабочую директорию.

Проанализируем качество чтений: fastqc SRR4240359.fastq.gz 2> 1.log. В выходном файле SRR4240359_fastqc.html видно, что длина чтений соответствует ожидаемой (лежит в диапазоне от 35 до 37, медиана — 36), в графе Overrepresented sequences только последовательности с большим количеством А или N, идущих подряд. Последнее объясняется крайне низким качеством концов ридов. Таким образом, адаптеры изначально удалены, всё хорошо.

На всякий случай (в т.ч. чтобы формально выполнить задание) удалим возможные остатки адаптеров. Создадим файл со списком всех адаптеров с помощью команды cat /mnt/scratch/NGS/adapters/*.fa >> /mnt/scratch/NGS/starlingsden/15pr/adapters.fa. Теперь можно убрать адаптеры: java -jar /usr/share/java/trimmomatic.jar SE -phred33 -threads 10 SRR4240359.fastq.gz 15.fastq.gz ILLUMINACLIP:adapters.fa:2:7:7 2> 2.log. Из лог-файла 2.log узнаём, что было удалено 55872 (0.41%) чтений. После триммирования файл с чтениями сократился с 445 до 443 мб.

Теперь удалим с правых концов чтений нуклеотиды с качеством ниже 20, оставив только чтения с длиной не меньше 32 нуклеотидов: java -jar /usr/share/java/trimmomatic.jar SE -phred33 -threads 10 15.fastq.gz 15_1.fatsq.gz TRAILING:20 MINLEN:32 2> 3.log. Из файла 3.log видно, что удалилось 1317986 (9.76%) чтений. В этот раз файл уменьшился до 385 мб. ПОсле этих манипуляций ситуация стала лучше: 15_1.fatsq_fastqc.html.

Подготовим k-меры длины 31: velveth velveth 31 -fastq.gz -short 15_1.fatsq.gz 2> 4.log. Во второй позиции стоит название директории, в которую сохранятся созданные файлы. 31 — длина k-меров, -fastq.gz задаёт формат входного файла, -short обозначает, что мы работаем с короткими одноконцевыми чтениями.

Воспользуемся программой velvetg, использующей граф де Брёйна (de Bruijn graph):velvetg velveth -cov_cutoff auto. Опция -cov_cutoff auto позволяет сразу убрать элементы с низким покрытием и получить оптимальные контиги.

В директории velveth получает следущие файлы: contigs.fa, Graph, LastGraph, Log, PreGraph, Roadmaps, Sequences, stats.txt. Из Log узнаём, что N50 имеет значение 125956. Импорт в гугл-таблицы stats.txt позволяет анализировать информацию о контигах.

Также есть два контига, выделяющиеся по покрытию относительно остальных. Это 39 и 71 контиги с покрытиями 4111220 и 1395 соответственно при стандартном покрытии в пределах 20-150.

Сравним программой megablast каждый из трёх самых длинных контигов с хромосомой Buchnera aphidicola. Для этого на странице blastn NCBI отметим "Align two or more sequences". В верхнее окно поместим AC генома, к нижнему прикрепим файл с контигом, выберем в окне Program Selection megablast (Изображение 1). Изначально файлы с отдельными контигами были получены командой seqretsplit contigs.fa -auto, из них выбрали файлы с самыми длинными контигами (Таблица 1).

Часть характеристик выравнивания бралась из раздела Descriptions а странице с выдачей blastn. Часть из скачанного в том же разделе (Download -> Hit Table (text)) файла, который был импортирован в таблицу (страницы 2, 8 и 3, по ID контига), чтобы легче было считать суммарное число однонуклеотидных различий и гэпов, а также сразу увидеть координаты начала и конца выравненных участков.

Теперь можно приступить непостредственно к анализу. В Таблице 2 указаны основные характеристики выравнивания. Видно, что 3 самых длинных контига суммарно соответствуют почти половине хромосомы. Несмотря на высокую идентичность, на Изображениях 2-4 видно заметное количество инделей, что может быть связано либо с параметрами запуска blastn, либо с тем фактом, что контиги Buchnera aphidicola str. Tuc7 (Acyrthosiphon pisum) сравнивались с хромосомой другого штамма (Buchnera aphidicola str. BAg (Aphis glycines)).