Учебный сайт Птицыной Елены

Cтудентки первого курса факультета биоинженерии и биоинформатики Московского государственного университета имени М.В. Ломоносова

Семестр 3, практикум 11

Назад на учебную страницу Птицыной Елены

Ресеквенирование. Поиск полиморфизмов у человека.

Задача этого практикума - найти и описать полиморфизмы у пациента (в хромосоме 16 в нашем случае).

Команды

$ cp /nfs/srv/databases/ngs/human/reads/chr16.fastq/nfs/srv/databases/ngs/elena-pt

$ cp /nfs/srv/databases/ngs/human/chr16.fasta/nfs/srv/databases/ngs/elena-pt

Загрузка файлов: чтений и референса

$ fastqc chr16.fastq

Анализ качества чтений до триммирования

$ java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar se -phred33 /nfs/srv/databases/ngs/human/reads/chr16.fastq chr16_trimmed.fastq TRAILING:20 MINLEN:50

Очистка чтений=триммирование с параметрами: TRAILING:20 (удаление нуклеотидов с конца с качеством менее 20), MINLEN:50 (установка минимальной длины чтения в 50 нуклеотидов)

результат: input reads: 3965 surviving: 3798 (95,79%) dropped: 167 (4,21%)

$ fastqc chr16_trimmed.fastq

Анализ качества чтений после триммирования

$ hisat2-build chr16.fasta mychr

Индексирование референса

$ hisat2 -x mychr -u chr16_trimmed.fastq -s chr16_aligntoref.sam --no-spliced-alignment --no-softclip

Картирование чтений на проиндексированную рефернсную последовательность с параметрами: --nosoftclip (выравнивание не должно допускать частичных несовпадений, то есть по всему риду).--no-spliced-alignment (выравнивание не должно допускать гэпов и инделей в риде).

результат:

3798 reads; of these:
3798 (100.00%) were unpaired; of these:
32 (0.84%) aligned 0 times
3763 (99.08%) aligned exactly 1 time
3 (0.08%) aligned >1 times
99.16% overall alignment rate

$ samtools view -b chr16_aligntoref.sam -o chr16_aligntoref.bam

Перевод из sam в bam (бинарный формат)

$ samtools sort chr16_aligntoref.bam chr16_aligntoref_sorted

Сортировка выравнивания по координате в референсе

$ samtools index chr16_aligntoref_sorted.bam

Индексирование отсортированного файла.

$ samtools flagstat chr16_aligntoref_sorted.bam

Статистика.

результат:

3804 + 0 in total (qc-passed reads + qc-failed reads)
6 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
3772 + 0 mapped (99.16%:-nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (-nan%:-nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (-nan%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapq>=5)

$ samtools mpileup -uf chr16.fasta -o chr16_polym.bcf chr16_aligntoref_sorted.bam

Получение файла с полиморфизмами.

$ bcftools call -cv -o chr16_polym.vcf chr16_polym.bcf

Изменение формата, вывод только изменившихся позиций.

$ convert2annovar.pl -format vcf4 chr16_polym.vcf -outfile chr16_polym.avinput

Конвертация в формат, распознаваемый командами дальнейших действий.

$ annotate_variation.pl -out ann_chr16_refgene -build hg19 -dbtyperefgene chr16_polym.avinput /nfs/srv/databases/annovar/humandb

Аннотация refgene

$ annotate_variation.pl -filter -out ann_chr16_dbsnp -build hg19 -dbtype snp138 chr16_polym.avinput /nfs/srv/databases/annovar/humandb.old/

Аннотация dbsnp

$ annotate_variation.pl -filter -dbtype 1000g2014oct_all -buildver hg19 -out ann_chr16_1000g chr16_polym.avinput /nfs/srv/databases/annovar/humandb.old/

Аннотация 1000 genomes

$ annotate_variation.pl -regionanno
-build hg19 -out ann_chr16_gwas -dbtype gwascatalog chr16_polym.avinput /nfs/srv/databases/annovar/humandb.old/

Аннотация gwas

$ annotate_variation.pl -filter -db
type clinvar_20150629 -buildver hg19 -out ann_chr16_clinvar chr16_polym.avinput /nfs/srv/databases/annovar/humandb.old

Аннотация clinvar

Исходное количество чтений

3965 (см. вывод Trimmomatic: input reads: 3965 surviving: 3798 (95,79%) dropped: 167 (4,21%))

До и после очистки

Программа trimmomatic отбросила 167 (4,21%) последовательностей.

Графики качества чтения оснований были сгненерированы fastqc (рис.1) На этих графиках красной линий обозначается медиана, желтый блок означает межквартильный диапазон, верхние и нижние усы - 10 и 90%, синяя линия - среднее значение. По оси y отложены quality scores. До очистки качество чтения было хорошее, но наблюдалось и выпадение в область плохих quality scores. После очистки ни один нуклеотид не имеет плохого качества прочтения, значит, очистка была оправдана.

Пример 1
Рисунок 1. Графики Per base sequence quality до и после очистки.

К сожалению, очистка trimmomatic не позволяет исправить gc content. Распределение в обычной случайной библиотеке должно быть нормальным, в нашем случае это не совсем так, много острых пиков. Это может свидетельствовать о загрязнении библиотеки.

Пример 2
Рисунок 2. GC content до и после очистки.

Картирование

3763 последовательностей aligned exactly 1 time. Можно сказать, что качество хорошее, так как это составляет 99.08% от всех последовательностей. Кроме того, только 3 последовательности были выровнены более 1 раза. Cм. вывод hisat2-build: 

3798 reads; of these:
3798 (100.00%) were unpaired; of these:
32 (0.84%) aligned 0 times
3763 (99.08%) aligned exactly 1 time
3 (0.08%) aligned >1 times
99.16% overall alignment rate

Описание полиморфизмов из .vcf

Всего в хромосоме 16 пациента 65 полиморфизмов, из них 64 замены и 1 вставка, делеций не было.

Примеры. В первой колонке таблицы ниже указано название хромосомы, далее позиция полиморфизма, его суть, качество чтения и глубина покрытия. 


chr16	11348273	T		A                       4.12853		DP=10
chr16	11361895	C		G                      	225.009		DP=50
chr16	11444454	gaaaaaaaaaaa	gaaaaaaaaaa,gaaaaaaaaa	9.03477		DP=55

Для качества чтений медиана 19,8099, среднее арифметическое 75,826144, для глубины покрытия медиана 3,00, среднее арифметическое 17,56.

Категории snp по база данных refseq в annovar

Эту информацию можно получить из файла ann_chr16_refgene.variant_function: downstream (3), exonic (9), intergenic (17), intronic (24), ncRNA_intronic (1), upstream (5), UTR3 (3), UTR5 (3).

Гены с snp, нуклеотидные и аминокислотные замены

Эту информацию можно получить из того же файла. Интересно посмотреть на уточнённый файл, охватывающий экзонные мутации - ann_chr16_refgene.exonic_variant_function: TNR2, PRM3, PRM2, PRM1, RMI2, PRSS53, HERPUD1. В ann_chr16_refgene.exonic_variant_function можно увидеть обозначение мутации в нуклеотидах и в аминокислотах (например, R131W означает замену 131-ого аргинина на триптофан. 

			
nonsynonymous SNV	TNP2:NM_005425:exon1:c.C391T:p.R131W,		G	A
stopgain		PRM3:NM_021247:exon1:c.C310T:p.R104X,		G	A
nonsynonymous SNV	PRM3:NM_021247:exon1:c.G299A:p.R100Q,		C	T
nonsynonymous SNV	PRM2:NM_001286359:exon1:c.C290T:p.P97L,		G	A
nonsynonymous SNV	PRM2:NM_001286359:exon1:c.C281T:p.A94V,		G	A
synonymous SNV		PRM1:NM_002761:exon2:c.C139A:p.R47R,		G	T
synonymous SNV		RMI2:NM_152308:exon2:c.A369C:p.T123T,		A	C
nonsynonymous SNV	PRSS53:NM_001039503:exon8:c.C1216G:p.P406A,	G	C
nonsynonymous SNV	HERPUD1:NM_001010989:exon2:c.G149A:p.R50H,
			HERPUD1:NM_001272103:exon2:c.G149A:p.R50H,
			HERPUD1:NM_014685:exon2:c.G149A:p.R50H,		G	A

snp c rs

Эту информацию можно получить из файла ann_chr16_dbsnp.hg19_snp138_dropped: 58 snp.

Средняя частота snp

Эту информацию можно получить из файла ann_chr16_1000g.hg19_ALL.sites.2014_10_dropped: медиана 0,4715455, среднее арифметическое 0,464222799 .

Клиническая аннотация snp

Эту информацию можно получить из ann_chr16_gwas.hg19_gwasCatalog: 3 SNP имеют клиническое значение. Они связаны с Obesity-related traits (экзонная синонимичная мутация гена PRM1), Parkinson's disease (интронная мутация PRSS53), Metabolic syndrome chr16 (экзонная несинонимичная мутация гена HERPUD1).

Аннотация clinvar

Эту информацию можно получить из ann_chr16_clinvar.hg19_clinvar_20150629_dropped. Ни один из изучаемых SNP не аннотирован.