Учебная страница курса биоинформатики,
год поступления 2014

Семестры Студенты Преподаватели

• 2024
• 2023
• 2022
• 2021
• 2020
• 2019
• 2018
• 2017

• ...

Практикум 13

Задача: Найти и описать полиморфизмы у пациента

Дано:

1. Чтения экзома, картирующиеся на участок хромосомы человека. Файлы с одноконцевыми чтениями в формате fastq лежат в директории /P/y14/term3/block4/SNP/reads. Распределение файлов по студентам см. в таблице.

2. Хромосомы человеческого генома (сборка версии hg19) лежат в директории /nfs/srv/databases/ngs/Human.

В отчёт включите

• Табличку со всеми выполненными командами:
  • команда (со всеми параметрами)
  • что делает
• Отчёты по частям 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. Картирование чтений.

• Откартируйте очищенные чтения с помощью программы BWA. Этапы

  • Сначала необходимо проиндексировать референсную последовательность; команда bwa index
  • Затем построить выравнивание прочтений и референса в формате .sam; используйте алгоритм mem; команда bwa mem.

Комментарий: программа BWA стоит на kodomo. Чтобы правильно запустить ее, изучите руководство. В разделе "Синопсис" приведены образцы подходящих команд.

Для интересующихся: программа BWA использует преобразование Барроуза — Уилера, про которое можно почитать, например, здесь.

4. Анализ выравнивания

• Переведите выравнивание чтений с референсом в бинарный формат .bam. Используйте пакет samtools, команда view: samtools view;

• Отсортируйте выравнивание чтений с референсом (получившийся после картирования .bam файл) по координате в референсе начала чтения; команда samtools sort;

• Проиндексируйте отсортированный .bam файл командой samtools index

• Выясните, сколько чтений откартировалось на геном; команда samtools idxstats.

Комментарий: программа 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. Как распределены Ваши чтения относительно экзонов гена?