Учебная страница курса биоинформатики,
год поступления 2016
Практикум 11
Задача: Найти и описать полиморфизмы у пациента
Дано:
1. Чтения экзома, картирующиеся на участок хромосомы человека. Файлы с одноконцевыми чтениями в формате fastq лежат на kodomo в директории /P/y14/term3/block4/SNP/reads. Распределение файлов по студентам см. в табл..
2. Хромосомы человеческого генома (сборка версии hg19) лежат в директории /nfs/srv/databases/ngs/Human на kodomo.
В отчёт включите
- Табличку со всеми выполненными командами:
- команда (со всеми параметрами)
- что делает
- Отчёты по частям I – III
Часть I: подготовка чтений
0. Создание рабочей директории.
В директории /nfs/srv/databases/ngs/ создайте свою директорию и скопируйте в нее Ваши файлы с ридами (.fastq) и хромосомой (.fasta).
1. Анализ качества чтений.
Сделайте контроль качества Ваших чтений с помощью программы FastQC.
Комментарий: программа FastQC установлена на kodomo, её можно вызвать командой "fastqc file.fastq", где file.fastq — имя файла с чтениями. Версию с графическим интерфейсом можно поставить на свой компьютер. В результате работы программы Вы получите архив (.zip), который содержит отчет о программе в виде html файла.
2. Очистка чтений
Сделайте очистку чтений с помощью программы Trimmomatic. Отрежьте с конца каждого чтения нуклеотиды с качеством ниже 20, оставьте только чтения длиной не меньше 50 нуклеотидов.
Комментарий: программа Trimmomatic установлена на kodomo. Вызывать её можно так:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 infile.fastq outfile.fastq step
где infile.fastq и outfile.fastq — входной и выходной файлы с чтениями, а step — выражение, указывающее, какую операцию производить.
Например, для удаления участков плохого качества можно вместо "step" написать SLIDINGWINDOW:10:28, что означает пройти по прочтениям окном длиной 10 и удалить правый конец каждого прочтения после окна со средним качеством меньше 28 (если такое окно найдётся). Почитайте руководство пользователя и выясните, как удалить плохие буквы с конца и как оставить только прочтения длины не менее 50. В чтениях, с которыми Вы работаете, адаптеры уже удалены.
- Сделайте анализ качества очищенных чтений с помощью FastQC; сравните с прежней выдачей FastQC.
В отчет включите:
- команды;
- картинки из FastQC "Per base quality" до и после чистки;
- число чтений до и после чистки;
- объяснение того, какие чтения или их части отсеялись и почему;
- (*ДОПОЛНИТЕЛЬНО) прокомментируйте различия до и после чистки по другим картинкам из выдачи FastQC; объясните что представлено на них и почему они характеризуют качество чтений; см. примеры выдач FastQC, показанные по ссылке выше.
Часть II: картирование чтений
3. Картирование чтений.
Откартируйте очищенные чтения с помощью программы Hisat2. Этапы
- Сначала необходимо проиндексировать референсную последовательность; команда hisat2-build
- Затем построить выравнивание прочтений и референса в формате .sam. Запустите hisat2 с параметрами --no-spliced-alignment и --no-softclip
- Сохраните вывод программы hisat2 в отдельный файл
Комментарий: Все необходимые для запуска программы Hisat2 файлы лежат тут: /home/students/y06/anastaisha_w/hisat2-2.0.5
4. Анализ выравнивания
Переведите выравнивание чтений с референсом в бинарный формат .bam. Используйте пакет samtools, команда view: samtools view;
Отсортируйте выравнивание чтений с референсом (получившийся после картирования .bam файл) по координате в референсе начала чтения; команда samtools sort;
Проиндексируйте отсортированный .bam файл командой samtools index
- Выясните, сколько чтений откартировано на геном; загляните в вывод программы Hisat2.
- Какую еще информацию выдает Hisat2 о выравнивании чтений на геном?
Комментарий: программа samtools также стоит на kodomo. Чтобы правильно запускать ее, изучите руководство.
- (* ДОПОЛНИТЕЛЬНОЕ) Найдите один экзон и определите его среднее покрытие чтениями
Команда samtools depth вычислит покрытие для каждого из нуклеотидов. Для подсказки запустите эту подпрограмму без параметров.
- Выберите один из хорошо покрытых прочтениями нуклеотидов. Можно построить гистограмму, чтобы понять, что такое хорошо покрытые нуклеотиды.
Найдите в GenomeBrowser участок, содержащий этот нуклеотид. Надо выбрать версию генома hg19!
- Определить координаты экзона, содержащего этот нуклеотид. Внимание: не пугайтесь, если нуклеотид не попал в экзон! Всякое бывает...
Перезапустите samtools depth, указав границы экзона, и найдите среднее покрытие.
В отчете укажите
- команды
- число чтений, картированных на хромосому, и число чтений, не картированных на хромосому
- если выполнили дополнительное задание, то укажите границы экзона и среднее покрытие нуклеотида экзона чтениями
Часть III: Анализ SNP
5. Поиск SNP и инделей.
Создайте файл с полиморфизмами в формате .bcf; команда samtools mpileup -uf. Опции и формат описаны в руководстве.
Создайте файл со списком отличий между референсом и чтениями в формате .vcf. Используйте команду "bcftools call -cv" пакета bcftools.
- Найдите и опишите в отчете три полиморфизма из .vcf файла. Для каждого приведите:
- кординату;
- тип полиморфизма: замена, вставка или делеция;
- что было в референсе, что найдено в чтениях;
- глубина покрытия данного места;
- качество чтений в данном месте.
- (* – дополнительно) После выполнения пункта 6 (аннотации) визуализируйте эти три полиморфизма с помощью IGV. Приведите картинки в отчёте.
Комментарий: для работы с программой IGV ознакомьтесь с руководством. Помните, что Вы работаете со сборкой генома человека версии hg19. Загрузите в программу отсортированный .bam файл с выравниванием. Сначала Вы не увидите никаких чтений, т.к. на экране будет представлен сразу весь геном. После работы с annovar Вы уже знаете, в какие гены попали Ваши чтения. В строке поиска IGV укажите один из Ваших генов, чтобы приблизить выравнивание так, чтобы было удобно смотреть на чтения.
6. Аннотация SNP.
- С помощью программы annovar проаннотируйте только полученные snp (индели не надо!). Используйте базы данных: refgene, dbsnp, 1000 genomes, GWAS, Clinvar.
Комментарий: программа установлена на kodomo: /nfs/srv/databases/annovar. Для работы с программой annovar из .vcf файла необходимо получить файл, с которым умеет работать эта программа. Сделать это можно с помощью скрипта convert2annovar.pl. См. руководство. Для аннотации файла с snp с помощью предложенных баз данных используйте скрипт annotate_variation.pl. В руководстве можно найти всю необходимую информацию о работе с программой. Например, узнать, какие из Ваших snp имеют rs, можно с помощью команды:
annotate_variation.pl -filter -out outputfile -build hg19 -dbtype snp138 inputfile.human humandb/
где inputfile.human — входной файл, полученный после обработки .vcf с помощью convert2annovar.pl (расширение не имеет значения); outputfile — выходной файл; humandb/ — директория, в которой лежат базы данных (все необходимые базы данных уже есть на kodomo, пользоваться опцией -downdb не надо!); snp138 — база данных, с которой вы работаете. Базы данных в annovar часто обновляются, для корректного запуска программы всегда нужно знать, какая версия какой базы данных у Вас скачена. Для вас: refgene — refGene; dbsnp — snp138; 1000 genomes — 1000g2014oct; GWAS — gwasCatalog; Clinvar — clinvar_20150629. В Annovar существуют 3 типа аннотаций по базам данных, основанных на: генной разметке (gene-based annotation); разметке других регионов генома (region-based annotation); фильтрации (filter-based annotation). Команды, с помощью которых можно проаннотировать полиморфизмы по необходимым базам данных:
refgene - gene-based annotation
annotate_variation.pl -out ex1 -build hg19 example/ex1.avinput humandb/
dbsnp - filter-based annotation
annotate_variation.pl -filter -out ex1 -build hg19 -dbtype snp138 example/ex1.avinput humandb/
1000 genomes - filter-based annotation
annotate_variation.pl -filter -dbtype 1000g2014oct_all -buildver hg19 -out ex1 example/ex1.avinput humandb/
Gwas - region-based annotation
annotate_variation.pl -regionanno -build hg19 -out ex1 -dbtype gwasCatalog example/ex1.avinput humandb/
Clinvar - filter-based annotation
annotate_variation.pl example/ex1.avinput humandb/ -filter -dbtype clinvar_20140211 -buildver hg19 -out ex1
Не забывайте, пожалуйста, прописывать правильные пути до тех файлов, к которым Вы обращаетесь!!!
Отчет
В отчете укажите
- команды
- описание трех полиморфизмов из .vcf файла
- cколько snp и сколько инделей Вы получили?
- Хорошее ли покрытие и качество у найденных полиморфизмов?
- На какие категории делит snp база данных refseq в annovar? Сколько snp у Вас попало в каждую группу?
- В какие гены попали Ваши snp?
- К каким нуклеотидным и аминокислотным заменам привели snp?
- Сколько snp имеет rs?
- Что Вы можете сказать о частоте найденных snp?
- Что Вы можете сказать о клинической аннотации snp?
- Составьте сводную таблицу, в которую бы вошли все snp и их характеристики по использованным для аннотации базам данных. По желанию можно выделить и/или описать наиболее интересные на Ваш взгляд полиморфизмы.
- (*) Поместите изображение выравнивания из IGV. Как распределены Ваши чтения относительно экзонов гена?
Комментарий: программы samtools и bcftools стоят на kodomo. Вам помогут опции "samtools mpileup -ugf" и "bcftools view -vcg". Итоговый файл со списком будет иметь расширение vcf. Итоговый файл со списком snp и инделей будет иметь расширение vcf. Опции и формат описаны в руководстве. Обратите внимание на информацию, касающуюся глубины покрытия каждого полиморфизма и его качества (Phred).