NGS. Первичная обработка данных

В этом блоке нам, наконец, предстояло поработать с реальными данными секвенирования следующего поколения (Next generation sequencing). В этом практикуме нам нужно было одготовить необходимые файлы, изучить качество чтений, отфильтровать чтения и подготовить их для картирования на референс. Работали мы с геномом человека (мне досталась 9-ая хромосома).

1. Описание файла с референсом

Файл представляет собой запись в классическом fasta-формате, хотя и расширение у него .fna. Первое, что бросается в глаза - это большое количество неизвестных нуклеотидов (хотя казалось бы, это референсная сборка):

>NC_000009.12 Homo sapiens chromosome 9, GRCh38.p13 Primary Assembly
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

Это связано, как я полагаю, со сложностями в секвенировании отдельных участков генома вроде теломер и центромер, которые состоят из многокопийных повторов. Я попытался оценить количество таких нуклеотидов:
wc -m chr9.fna
140124720 - количество нуклеотидов (можно пренебречь парой десятков символов из заголовка)
grep 'N' chr9.fna | wc -m
16815099 - количество "N", которое составляет аж 12% от всех.

2. Индексация референса

Теперь из этой последовательности мы должны получить индекс, который в дальнейшем пригодится нам для маппирования на него чтений. Запускаем команду:
bwa index -a bwtsw chr9.fna, где параметр -a - алгоритм построения BWT индекса с опцией bwtsw (алгоритмом BWT-SW, адаптированным по скорости работы и занимаемой памяти к индексации генома человека). В результате в директории с файлом нашей последовательности появились файлы с расширениями: .amb, .ann, .bwt, .pac, .sa. Всё это вспомогательные файлы, которые будет использовать алгоритм для картирования.

3. Описание образца

Теперь нужно подготовить чтения секвенатора. Вот что удалось о них найти:
a. Ссылка на данные в базе NCBI;
b. Прибор - Illumina Genome Analyzer IIx;
c. Организм - Homo sapiens;
d. Стратегия секвенирования - Whole-exome;
e. Тип чтений - парноконцевые;
f. Ожидаемое количество чтений (spots) - 41 398 792.

4. Проверка качества исходных чтений

Теперь проанализируем качество исходных чтений с помощью программы FastQC.
a. Использованные команды:
fastqc SRR10720412_1.fastq.gz - отчёт
fastqc SRR10720412_2.fastq.gz - отчёт
b. Всего получилось 41 398 792 пары чтений;
c. У прямых и обратных чтений количество ридов совпадает;
d. Также количество пар чтений совпадает с числом ожидаемых из пункта 3f;
e.

f. Качество у чтений отменное, что неудивительно при количестве циклов элонгации в 75. В другой ситуации я бы даже их не чистил. Также наблюдаем типичную картину - качество постепенно уменьшается к концу. В среднем качество варьируется от 30 до 40 единиц;
g.

h. Отчёт о показателях качества для каждой последовательности позволяет нам увидеть, имеют ли подпоследовательности наших последовательностей низкие значения качества. Если бы они были, мы бы увидели заметный локальный пик слева на графике, что означало бы системную проблему, скорее всего, с проточной ячейкой;
i.

j. Этот график показывает распределение размеров фрагментов в анализируемом файле. На рисунке мы видим, что в файлах преобладают последовательности длиной 75 b.p., что говорит о хорошей пробоподготовке и качественной фрагментации DNA образца, хотя чаще фрагментируют на куски подлиннее, например, 300.

5. Фильтрация чтений

После проверки качества логичным этапом следует тримминг чтений. Для этого мы воспользуемся алгоритмом Trimmomatic в парноконцевом режиме. На вход подаём 2 файла с чтениями, а на выходе получаем 4 файла: 2 paired и 2 unpaired. Также отметим, что этот прибор кодирует качество в формате phred-score = 33. Для очистки удалим с конца чтений нуклеотиды с качеством ниже 20 (trailing = 20) и оставим среди них только такие чтения, длина которых не ниже 50 нуклеотидов (minlen = 50). Для ускорения работы возьмем побольше ядер (threads = 10). Прямые прочтения обозначены как f (forward), обратные - как r (reverse), спаренные как p (paired), а неспаренные как u (unpaired). Итоговая команда выглядит следующим образом:
java -jar /usr/share/java/trimmomatic.jar PE -threads 10 -phred33 -trimlog trimlog_dna ../dna_data/SRR10720412_1.fastq.gz ../dna_data/SRR10720412_2.fastq.gz dna_f_p.fastq.gz dna_f_u.fastq.gz dna_r_p.fastq.gz dna_r_u.fastq.gz TRAILING:20 MINLEN:50
Подав на вход 41 398 792 пар прочтений мы получили 39 853 177 (96.27%) парных, 678 011 (1.64%) непарных и 659625 (1.59%) обратных последовательностей. Всего было удалено 207 979 (0.50%) ридов.

6. Проверка качества триммированных чтений

Проверим теперь качество ридов после очистки. Ожидаем, что заметных изменений не будет.
a. Использованные команды:
fastqc dna_f_p.fastq.gz - отчёт
fastqc dna_f_u.fastq.gz - отчёт
fastqc dna_r_p.fastq.gz - отчёт
fastqc dna_r_u.fastq.gz - отчёт
b. В forward и reverse paired осталось 39 853 177 чтений, что составляет 96,2% от прежнего числа (см. пункт 4b). Отметим, что это совпадает с данными Trimmomatic;
c.

d. В глаза бросается различие по качеству между спаренными и неспаренными ридами. Вторые, очевидным образом, хуже. Я думаю, это связано с тем, что по у спаренных происходила взаимная коррекция ошибок и уточнение, примерно как мы делали в практикуме с секвенированием по Сэнгеру;
e. А вот спаренные, после очистки, по понятным причинам, выглядят лучше, чем до очистки;
f.

g. В этой графе после очистки уменьшился левый хвост, потому что мы убрали совсем короткие прочтения и почистили риды с конца;
h.

i. Эту картинку можно объяснить так: после очистки по качеству с конца появились следовые количества ридов, длиной до 75, которых не было до этого, а затем риды длиной меньше 50 (это видно по горизонтальной оси) были удалены.