За референс взята последовательность 10й хромосомы человека, которая находится в файле chr10.fna. Данный файл содержит нуклеотидную последовательность в .fasta формате, в его начале имеется строка >NC_000010.11 Homo sapiens chromosome 10, GRCh38.p13 Primary Assembly, где NC_000010.11 – это идентификатор последовательности, Homo sapiens chromosome 10 – 10-я хромосома человека, GRCh38.p13 Primary Assembly – первичная сборка генома (всё того же человека).
Сама последовательность в начале и в конце содержит неопределенные нуклеотиды N, что можно списать на то, что там находятся теломерные последовательности.
Для индексации референса была использована команда
bwa index -a bwtsw chr10.fna
где команда bwa index индексирует базу данных (наш референс) chr10.fna с помощью алгоритма bwtsw. Если с помощью параметра -a не указать конкретный алгоритм, по дефолту будет использоваться IS в силу своей простоты.
В результате были получены файлы chr10.fna.amb(1), chr10.fna.ann(2), chr10.fna.bwt(3), chr10.fna.pac(4) и chr10.fna.sa(5), причем 3 последних – бинарные.
Два файла с чтениями ДНК были скопированы в отдельную директорию, где далее изучались на предмет качества содержащихся в них чтений. Далее была введена команда
fastqc *.gz
в результате работы которой появились 4 файла, из которых было найдено, что количество прямых и обратных чтений совпадает и равно ожидаемому числу чтений (40990418).
Из Fig.1 видно, что качество прочтений высокое, однако к концу заметно его снижение в обоих случаях, что является довольно типичной ситуацией. Также можно отметить, что в среднем качество прямого прочтения выше, чем качество обратного.
На Fig.2 представлена зависимость количества чтений (ось ординат) от качества (ось абсцисс). Как можно видеть, пик максимальное количество ридов как прямых, так и обратных имеют качество, близкое к отметке 40. Если сравнивать графики для прямых и обратных прочтений между собой, то можно заметить, что у первых качество выше, поскольку во втором случае есть некоторое количество совершенно некачественных ридов.
Из данных, представленных на Fig.3, легко можно видеть, что в среднем длина прочтения в обоих случаях составила 75 bp.
Команда:
java -jar /usr/share/java/trimmomatic.jar PE -threads 10 -phred33 -trimlog trim_log SRR10720416_1.fastq.gz SRR10720416_2.fastq.gz paired_out_1.fq.gz unpaired_out_1.fq.gz paired_out_2.fq.gz unpaired_out_2.fq.gz TRAILING:20 MINLEN:50
Качество триммированных чтений проверялось уже знакомой командой fastqc:
fastqc *
* использовалась, поскольку в директории находились только нужные файлы (файл trim_log был удален, поскольку слишком много весил, а полезной информации содержал не столь много).
Количесвто оставшихся пар чтений как forward (1, fw), так и reverse (2, rev) равно 38977379, т.е. сохранилось 95.1% исходных чтений.
На Fig.4 явно видно, что после триммирования качество прочтений (paired) улучшилось как для fw, так и для rev, поскольку теперь даже ближе к концу усы (rev) лежат в зеленой зоне. Также очевидно, что качество paired на порядок выше, чем unpaired.
Если сравнивать результаты, представленные на Fig.5, и сравнить их с представленными на Fig.2, то видно, что глобально изменилось только значение качества на пике: было 38, стало 39, т.е. триммирование было не напрасным занятием. Количество ридов сильно не изменилось.
При сравнении результатов, представленных на Fig.6 и Fig.3, выясняется, что в среднем длина чтений после триммирования не изменилась и большинство до сих пор имеет длину 75 bp, однако в обоих случаях появилось некоторое количество коротких ридов с длиной от 49 bp. Интересно также, что, если присмотреться, на Fig.6 можно заметить в обоих случаях небольшие увеличения количества ридов с длиной 55, 60, 65 и 70, причем чем ближе к пику, тем относительно больше ридов с соответствующей длиной.