Нахождение и описание полиморфизмов у человека

Часть I: подготовка чтений

Была создана директория /nfs/srv/databases/ngs/asya.kalashnikova, куда были скопированы файлы chr8.fastq и chr8.fasta - с ридами и хромосомой. Далее производился анализ качества чтений с помощью программы FastQC на Kodomo. Была использована команда (1). В итоге были получены 2 файла: chr8_fastqc.html и chr8_fastqc.zip . Проводилась очистка чтений с помощью программы Trimmomatic, причем необходимо оставить только чтения длиной не меньше 50 и отрезать с конца каждого чтения нуклеотиды с качеством ниже 20.
Использовались определенные операции:
1) TRAILING: отрезает основания в конце прочтения, качество которых будет ниже заданного значения;
2) MINLEN: убирает чтения, чья длина менее заданной.
(1) fastqc chr8.fastq
(2) java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr8.fastq 1.fastq TRAILING:20 MINLEN:50
(3) fastqc 1.fastq
Далее был проведен анализ качества очищенных чтений (файл: 1.fastq) с помощью команды (3). В итоге было получено 2 файла: 1_fastqc.html и 1_fastqc.zip .
Число чтений до и после чистки: 8367 (до), 8227 (после).

Рис. 1 - Per base sequence quality (до и после чистки) - увеличение при нажатии

Рис. 2 - Basic Statistics (до и после чистки) - увеличение при нажатии


Как можно увидеть на рис.1 очистка существенно улучшила показатель качества ридов, к тому же подняла ценз на длину - как и предполагалось (см. рис. 2). Показатель "GC Content" практически не изменяется. Показатель длины чтений до чистки составлял 34-100, после очистки - 50-100. Таким образом, отсев ридов происходил, как и предполагалось по качеству и длине рида.

Часть II: картирование чтений

1) Необходимо было откартировать чтения с помощью программы BWA. Для этого необходимо было проиндексировать референсную последовательность с помощью команды (1). Затем было построено выравнивание прочтений и референса в формате sam с помощью алгоритма mem команды (2).
2) Далее происходил анализ выравнивания: сначала выравнивание было переведено в бинарный формат bam с помощью пакета samtools, команды (3). Опция -o была использована для указания выходного файла, -b - для смены формата с sam на bam. Получившийся файл был отсортирован по координате в референсе начала чтения - команда (4). Далее файл sort.bam был проиндексирован командой (5) и исследован на предмет числа откартированных чтений на геном с помощью команды (6). Был получен файл chr8_car.txt .
Число не картированных/картированных чтений на хромосому: 3/8224 из 8227.
(1) bwa index chr8.fasta
(2) bwa mem chr8.fasta 1.fastq > align.sam
(3) samtools view align.sam -b -o align.bam
(4) samtools sort align.bam sort
(5) samtools index sort.bam
(6) samtools idxstats sort.bam > chr8_car.txt

Часть III: Анализ SNP

1) Поиск SNP и инделей: с помощью команды (1) был создан файл с полиморфизмами с формате bcf, а с помощью (2) был создан файл со списком отличий между референсом и чтениями в формате vcf. Был получен файл read.vcf. Число snp: 94; число инделей: 8. Далее необходимо было найти и описать 3 полиморфизма из vcf файла - таблица 1.

Таблица 1. Описание полиморфизмов

Полиморфизм Координата Качество чтений Глубина покрытия Было в референсе Найдено в чтениях Тип полиморфизма
№1 27466315 154.998 29 T C Замена
№2 116424270 59.456 31 A AC Вставка
№3 76452313 221.999 30 G A Замена
Покрытие: min - 1; max - 134; среднее - 14,8.
Качество: min - 3,01618; max - 225,009; среднее - 68,1076.
(Я не совсем представляю себе критерии "хороших" покрытия и качества, но предполагаю, что среднее значение качества являет довольно неплохим, в то время как среднее значение глубины покрытия я бы назвала "средним").
SNP привели к данным нуклеотидным заменам:
T - (C, A, G); C - (T, G, A); G - (A, T); A - (G, T, C).
2) Аннотация SNP:
с помощью annovar необходимо было проаннотировать полученные snp по базам данных refgene, dbsnp, 1000 genomes, GWAS, Clinvar. Сначала был получен файл без инделей: No_indels.vcf. Потом было необходимо получить файл формата, с которым умеет работать annovar с помощью скрипта convert2annovar.pl, который был скопирован в директорию с остальными файлами. Это было сделано командой (3).
Dbsnp:
аннотация была получена с помощью команды (4). Данная база используется для того, чтобы узнать, какие из snp имеют rs. Аннотация по этой базе данных основана на фильтрации (filter-based annotation). На выходе были получены 3 файла: dbsnp.hg19_snp138_dropped (файл с snp, имеющими rs), dbsnp.hg19_snp138_filtered (файл с snp, не имеющих rs), dbsnp.log (файл, содержащий информацию по запуску). Таким образом, 76 snp из 94 имеют rs (18 не имеют).
Refgene:
аннотация была получена с помощью команды (5). Аннотация по этой базе данных основана на генной разметке (gene-based annotation). На выходе 3 файла: refgene.exonic (файл содержит snp, найденные в экзонах), refgene.log (файл с информацией о запуске и работе программы), refgene.variant_function (файл, содержащий все snp).
Категории: 1) Exonic: 5; 2) Intronic: 58; 3) Intergenic: 18; 4) UTR3: 13.
Гены, в которые попали snp: TRPS1 (exonic, UTR3, intronic), CLU (UTR3, intronic, exonic), CASC9 (intergenic), HNF4G (intergenic, intronic, exonic, UTR3).
1000 genomes:
аннотация была получена с помощью команды (6). Аннотация по этой базе данных основана на фильтрации (filter-based annotation). Выход похож на выход dbsnp, но по базе 1000 genomes. На выходе получены 3 файла: 1000gen.hg19_ALL.sites.2014_10_filtered, 1000gen.log (файл с информацией по запуску), 1000gen.hg19_ALL.sites.2014_10_dropped. Аналогично, не было аннотировано 18 snp.
Частота snp: максимальная - 1; минимальная - 0,00499201; средняя - 0,466446.
Gwas:
аннотация была получена с помощью команды (7). Аннотация по этой базе данных основана на разметке других регионов генома (region-based annotation). Данная база аннотирует snp, ассоциируемые с определенными заболеваниями, то есть база отвечает за клиническую аннотацию. На выходе было получено 2 файла: gwas.log (файл с информацией по запуску), gwas.hg19_gwasCatalog (файл с информацией по болезням). Было аннотировано 4 snp: 2 ассоциируемы с болезнью Альцгеймера, 1 - с "хорошим" холестерином (HDL cholesterol), 1 - с уровнем мочевой кислоты.
Clinvar:
аннотация была получена с помощью команды (8). Аннотация по этой базе данных основана на фильтрации (filter-based annotation). На выход получаем 3 файла, аналогичные выходу dbsnp, но в данном случае аннотированных snp не наблюдается. Файл с параметрами запуска: clinvar.log; файл с неанотированными snp (всеми): clinvar_filtered.
Таблица snp с аннотациями по базам данных
(1) samtools mpileup -uf chr8.fasta sort.bam -o polymorf.bcf
(2) bcftools call -cv polymorf.bcf -o read.vcf
(3) perl convert2annovar.pl -format vcf4 No_indels.vcf -outfile No_indels.avinput
(4) perl annotate_variation.pl -filter -out dbsnp -build hg19 -dbtype snp138 No_indels.avinput /nfs/srv/databases/annovar/humandb
(5) perl annotate_variation.pl -out refgene -build hg19 No_indels.avinput /nfs/srv/databases/annovar/humandb
(6) perl annotate_variation.pl -filter -dbtype 1000g2014oct_all -buildver hg19 -out 1000gen No_indels.avinput /nfs/srv/databases/annovar/humandb
(7) perl annotate_variation.pl -regionanno -build hg19 -out gwas -dbtype gwasCatalog No_indels.avinput /nfs/srv/databases/annovar/humandb
(8) perl annotate_variation.pl -filter -dbtype clinvar_20150629 -buildver hg19 -out clinvar No_indels.avinput /nfs/srv/databases/annovar/humandb

© Kalashnikova Anastasia, 2016