Практикум 11. Анализ данных NGS: референс, качество чтений, триммирование
Этот практикум - первая часть домашнего задания по анализу данных NGS. Здесь описаны процедуры индексирования референсной последовательности при помощи bwa, а также проверки качества прочтений и их триммирования.
Описание файла с референсом
В моём случае референсом будет служить последовательность четвёртой хромосомы человека. Ниже приведены первые несколько строк файла chr4.fna, содержащего референс.
>NC_000004.12 Homo sapiens chromosome 4, GRCh38.p13 Primary Assembly NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
По-видимому, это файл в формате fasta, в котором имеются заголовок с названием и описанием и собственно последовательность. Некоторое количество позиций в различных участках последовательности, в частности, в начале и в конце, заполнено буквами неопределённого нуклеотида - «N».
Индексация референса
Референсная последовательность была проиндексирована при помощи программы bwa. Была использована следующая команда:
bwa index -a bwtsw chr4.fna &> bwa_index.log
Здесь bwa index - команда для индексирования референса, -a bwtsw - опция, задающая алгоритм, используемый при индексировании (в описании аргумента bwtsw сказано, что «этот метод работает с цельным человеческим геномом»), а chr4.fna - файл с последовательностью, принимаемый на вход. Можно также видеть, что вся выдача программы в консоль (как stdout, так и stderr) перенаправлялась в log-файл. На выходе получилось 5 файлов - chr4.fna.amb, chr4.fna.ann, chr4.fna.bwt, chr4.fna.pac и chr4.fna.sa. В первом из них содержится информация об участках, состоящих из неизвестных нуклеотидов - «N», а именно: их число, координаты и длина; во втором файле - немного информации о входных данных: название последовательности, её длина и пр.; остальные файлы - бинарные.
Описание образца
Поиск SRR10720416 по базе данных NCBI SRA дал следующий результат: ссылка. Ниже приведена некоторая информация об эксперименте.
• Прибор: Illumina Genome Analyzer IIx
• Организм: Homo sapiens
• Стратегия секвенирования: экзомное
• Прочтения: парноконцевые
• Ожидаемое число прочтений: 40 990 418
Проверка качества исходных чтений
Для этой цели была использована программа fastqc. На вход ей были поданы файлы с прямым и обратным прочтениями - SRR10720416_1.fastq.gz и SRR10720416_2.fastq.gz соответственно. Конкретно, была выполнена команда:
fastqc SRR10720416_?.fastq.gz &> fastqc_run.log
С логами её выполнения можно ознакомиться здесь; результат работы программы, так удобно представленный в виде html-страниц, доступен по ссылкам:
SRR10720416_1_fastqc SRR10720416_2_fastqc
Далее - кратко о выдаче программы:
• Получилось 40 990 418 пар прочтений
• Это число одинаково для прямых и обратных прочтений
• Это число совпадает с полученным из базы данных в предыдущем задании
На рисунках 1 и 2 представлены графики качества прочтения в зависимости от позиции нуклеотида. Можно видеть, что как для прямого, так и для обратного прочтения качество несколько уменьшается ближе к концу последовательности, но всё же остаётся достаточно высоким (так, нижний квартиль всегда не менее 33).
На рисунках 3 и 4 изображены графики, иллюстрирующие количество прочтений с различным средним качеством. Больше всего прочтений имеют среднее качество Q=39, однако наблюдается также некоторое число обратных прочтений со средним качеством 2.
На рисунках 5 и 6 можно видеть графики распределения длины прочтений. В данном случае они крайне просты - все прочтения имеют длину 75 bp.
Фильтрация чтений
Для фильтрации прочтений была использована программа trimmomatic. Команда приведена далее:
java -jar /usr/share/java/trimmomatic.jar PE -phred33 -threads 10 SRR10720416_1.fastq.gz SRR10720416_2.fastq.gz trimmed_reads/SRR10720416_1_paired.fastq.gz trimmed_reads/SRR10720416_1_unpaired.fastq.gz trimmed_reads/SRR10720416_2_paired.fastq.gz trimmed_reads/SRR10720416_2_unpaired.fastq.gz TRAILING:20 MINLEN:50 &> trimmomatic_run.log
Здесь java -jar /usr/share/java/trimmomatic.jar - команда запуска исполняемого jar-файла trimmomatic с указанием его полного адреса, PE - аргумент для работы с парноконцевыми
чтениями, -phred33 - опция, указывающая формат Phred Score,
Проверка качества триммированных чтений
Программа fastqc была применена ещё раз, теперь - к полученным в предыдущем задании триммированным прочтениям. Команда приведена далее; логи, как обычно, лежат тут.
fastqc SRR10720416_?_*.fastq.gz &> fastqc_run_trimmed.log
По ссылкам ниже можно просмотреть результат работы программы:
SRR10720416_1_paired SRR10720416_1_unpaird SRR10720416_2_paired SRR10720416_2_unpaired
Краткое описание того, что получилось:
• Количество парных прочтений, оставшихся после триммирования: 38 977 379
• Это составляет 95.09% от исходного числа
На рисунках 7-10 изображены графики качества прочтений по нуклеотидам для триммированных прямых/обратных парных/непарных прочтений. Можно видеть, что качество парных прочтений (как прямых, так и обратных) после триммирования улучшилось по сравнению с исходными прочтениями (рисунки 1, 2). В то же время качество непарных прочтений значительно ниже, чем парных - так, в них встречаются позиции с Phred Score менее 28 (а в обратных непарных - менее 20).
На рисунках 11 и 12 приведены графики, отражающие количество парных триммированных прочтений с различным средним качеством. Графики несколько напоминают таковые для исходных прочтений (рисунки 3 и 4), но можно заметить, что качество всё же улучшилось, так как, судя по всему, отсутствуют триммированные прочтения с Q<14 (среди исходных они встречались).
На рисунках 13 и 14 можно видеть графики распределения по длине триммированных парных прочтений. Разница с рисунками 5 и 6 состоит в появлении некоторого количества прочтений длины менее 75 bp - это связано с удалением некачественных концевых участков некоторых прочтений, при этом прочтения длины менее 50 по-прежнему отсутствуют в связи с параметрами запуска trimmomatic. Любопытно, что маленькие дополнительные пики на графике соответствуют длинам прочтений в 70, 65, 60 и 55 bp - похоже, программа убирала нуклеотиды низкого качества преимущественно по 5 за раз.