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

1. Анализ качества чтений и их очистка

Мне была выдана 7-я хромосома человеческого генома (сборка версии hg19) chr7.fasta и, соответственно, файл с ридами для этой хромосомы chr7.fastq (без адаптеров).

В начале был сделан контроль качества ридов программой FastQC. В результате работы которой был получен архив (.zip) и отчет о результатах работы программы в виде html страницы [ chr7_fastqc.html ]

команда: fastqc chr7.fastq

Затем с помощью программы Trimmomatic была произведена очистка чтений:

команда: java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr7.fastq chr7_out.fastq TRAILING:20 MINLEN:50

Таким образом, в выходном файле chr7_out.fastq были оставлены только чтения длиной не меньше 50 нуклеотидов (MINLEN:50) и с конца каждого рида были удалены нуклеотиды с качеством ниже 20 (TRAILING:20).

Программой FastQC я оценила качество чтения, полученного после очистки:

команда: fastqc chr7_out.fastq

В результате работы программы был получен архив (.zip) и отчет о результатах работы программы в виде html страницы [ chr7_out_fastqc.html ]

ниже приведены основные характеристики изначального чтения и чтения, полученного после очистки программой Trimmomatic):

Basic Statistics

до чистки

после чистки

Per base sequence quality

до чистки

после чистки

Как видно, после чистки показатели качества в ридах стали намного лучше. Таким образом, были удалены те чтения, качество которых нас не устраивает и оставлены только те, с которыми можно работать дальше.

В результате чистки из 3752 ридов осталось 3650 (т.е. было удалено 102 рида, что составляет 2,7% от изначального количества)

Помимо приведенных в отчете изображений в архивах содержатся и другие.

(адекватные результаты обозначены — ✔; сомнительные результаты (warning) — !; когда все плохо (failure) — ✖):

✔Per sequence quality scores — распределение среднего качества среди ридов. В нашем случае все в порядке: пик на 38, т.е. качество у большинства ридов очень хорошее.

✔Per base sequence content — процент каждого нуклеотида в ридах. В случайной библиотеке процент разных нуклеотидов для всех позиций, по идее, должен быть примерно одинаковым, так что линии на графике должны быть параллельными и лежащими в одном диапазоне (чем параллельнее и ближе друг к другу, тем лучше (равнее) распределение). Если каких-то нуклеотидов в последовательности, в целом, больше, то графически это отразится в более высоком расположении кривой. Пики на графике обусловлены распределением нуклеотидов в ридах (в одном положении превалирует одно основание, в других – другое). В нашем случае неровности относительно несильны, а линии графика лежат в одном диапазоне, так что, можно считать, что все в порядке.

! Per sequence GC content — распределение GC-состава для всех ридов. Для нормальной случайной библиотеки распределение GC-состава должно быть нормальным, т.е. центральный пик должен совпадать с пиком для GC-состава всего генома. Т.к. исходного генома нет, вместо него используется геном, построенный с помощью расчетов из имеющихся данных. В нашем случае этот модуль отмечен как warning, т.к. отклонение от нормального (гауссовского) распределения существенно; и хотя в графике прослеживается относительная симметричность. (Притом медиана реального распределения не сильно отличается от медианы теоретического, если вообще отличается. То же самое можно сказать и о дисперсии: график реального распределения не намного выше графика распределения теоретического). Однако, на определенных участках ридов GC-состав относительно других положений значительно отличается (больше/меньше теоретической нормы: отсюда пики на графике).

✔Per base N content — график, отражающий позиции, в которых был неопределен нуклеотид (N). В нашем случае все хорошо

! Sequence Length Distribution — график, отображающий распределение длин ридов. В нашем случае модуль обозначен warning.

✔Sequence Duplication Levels

✔Overrepresented sequences

✔Adapter Content

✔Kmer Content

2. Картирование чтений. Анализ выравнивания

Далее надо было картировать чтения на имеющуюся референсную хромосому.

Была проиндексирована референсная последовательность

команда: bwa index chr7.fasta

После чего было построено выравнивание прочтений chr7_out.fastq и проиндексированного референса в формате (.sam) с использованием алгоритма mem:

команда: bwa mem chr7.fasta chr7_out.fastq > vyravnivanie.sam

Полученный в результате файл: vyravnivanie.sam

После этого с помощью пакета samtools файл с выравниванием был переведен в бинарный формат (.bam)

команда: samtools view vyravniv.sam -b -o vyravnivanie.bam

(-o - выходной файл; -b - бинарный формат)

В результате получился файл: vyravnivanie.bam

Получившийся файл был подан на вход другой команде, чтобы отсортировать выравнивание чтений с референсом по координате в референсе начала чтения:

команда: samtools sort vyravnivanie.bam -T vrem.txt -o sort.bam

(-T - создает временный файл(Параметр -T был задан для записи временных файлов в файл vrem.txt, т.к. иначе они выдавались в stdout); -o - выходной файл)

Полученный в результате файл: sort.bam

Далее полученный файл был проиндексирован:

команда: samtools index sort.bam

А затем было выяснено, сколько чтений откартировалось на геном:

команда: samtools idxstats sort.bam > map.out

В результате получился файл: map.out

полученный результат, записанный в файл map.out :

Таким образом, на геном откартировалось 3648 ридов из 3650 и не откартировались только 2.

3. Анализ SNP

Поиск SNP и инделей

Был создан файл с полиморфизмами в формате (.bcf)

команда: samtools mpileup -uf chr7.fasta sort.bam > snp.bcf

В результате получился файл: snp.bcf

Полученный файл snp.bcf подан на вход другой команде:

bcftools call -cv snp.bcf > snp.vcf

в результате чего был получен файл snp.vcf

Файл содержит список отличий между референсом и ридами. Всего было найдено 28 SNP и 3 вставки. Качество SNP варьирует ~ в пределах от 3 до 225, среднее значение — 115. При этом для 24 находок качество >20 (хорошее). Максимальное покрытие составляет 71, а для 5 SNP покрытие составляет 1 рид.

В следующей таблице приведено описание для трех полиморфизмов из этого файла.

Координата	Тип полиморфизма	В референсе	В ридах	Глубина покрытия	Качество ридов
100487721	Замена	G	T	4	38
120965718	Индель (вставка)	taa	tAAaa	2	17
120979525	индель (вставка)	CCT	CCTCTCT	4	144

4. Аннотация SNP

В этой части задания необходимо было с помощью программы annovar проаннотировать только полученные SNP (без инделей, индели были вручную удалены из файла snp.vcf , результат был сохранен как snp_nonindel.vcf).

Использовались следующие базы данных:

Название базы	Версия на kodomo
refgene	refGene
dbsnp	snp138
1000 genomes	1000g2014oct
GWAS	gwasCatalog
Clinvar	clinvar_20150629

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

Для начала с помощью скрипта convert2annovar.pl файл snp_nonindel.vcf был переведен в формат avinput, с которым может работать annovar

команда: perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 snp_nonindel.vcf -outfile snp.avinput

Полученный файл: snp.avinput

Далее полученный файл использовался для аннотации SNP с помощью скрипта annotate_variation.pl.

Аннотация по базе refgene:

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

команда: perl /nfs/srv/databases/annovar/annotate_variation.pl -out chr7_refgene -buildver hg19 snp.avinput /nfs/srv/databases/annovar/humandb/

В результате было получено 3 файла:

chr7_refgene.variant_function (содержит информацию обо всех найденных SNP),

chr7_refgene.exonic_variant_function (содержит описание SNP, попавших в экзоны)

chr7_refgene.log (содержит информацию о процессе работы программы).

*по какой-то причине файл chr7.refgene.log изначально не отобразился

В следующей таблице приведена общая информация о найденных SNP из файла chr7.refgene.variant_function. В нем SNP разделены на несколько групп в зависимости от положения в геноме (интроны, экзоны, UTR и т.д.). Также в этом файле указано, какой является замена: гетерозиготной или гомозиготной.

intronic	exonic	UTR3	hom	het
22	4	2	9	19

Как видно из таблицы, больше всего замен было найдено в интронах, а меньше всего — в нетранслируемых областях, в экзонах их тоже немного, что, в принципе, логично: замены в кодирующих последовательностях могут привести к замене аминокислоты в соответствующей белковой последовательности, что может привести к нарушению функции белка.

В область экзонов попало 4 SNP, краткая информация о которых приведена в следующей таблице

Координата	Тип замены	Ген	В референсе -> в ридах	Замена аминокислоты	Глубина покрытия	Качество ридов
100490077	синонимичная	ACHE	G -> A	C1431T:p.P477P C1167T:p.P389P	34	221.999
120969769	несинонимичная	WNT16	G -> A	G214A:p.G72R G244A:p.G82R	14	222.974
120979089	несинонимичная	WNT16	C -> T	C758T:p.T253I C788T:p.T263I	13	225.2
134264286	синонимичная	AKR1B15	C -> T	C1020T:p.F340F	42	187.009

При этом половина всех замен нуклеотидов в экзонах синонимична, т.е. они не привели к замене аминокислоты в белке на другую.

информация об упомянутых генах и белках, которые они кодируют:

ACHE - ген, кодирующий фермент ацетилхолинэстеразу. Ацетилхолинэстераза (сокр. АХЭ, англ. AChE) — гидролитический фермент , из семейства эстераз, который содержится в синапсах и катализирует гидролиз нейромедиатора ацетилхолина до холина и остатка уксусной кислоты . Реакция, катализируемая ацетилхолинэстеразой, необходима для дезактивации ацетилхолина в синаптической щели и перехода клетки-мишени в состояние покоя (например, для расслабления мышечной клетки). Известно, что АХЭ играет важную роль в патогенезе таких заболеваний как болезнь Альцгеймера, черепно-мозговые травмы, сосудисто-мозговые патологии, миастения Гравис и т.д.

WNT16 – ген, кодирующий белок Wnt-16. Семейство генов WNT состоит из структурно родственных генов, кодирующих сигнальные белки. Эти белки участвуют в онкогенезе и в некоторых процессах развития, включая дифференциацию клеток в процессе эмбриогенеза. Предполагается, что стимуляция экспрессии WNT16 в близлежащих нормальных клетках отвечает за развитие устойчивости к химиотерапии у раковых клеток.

AKR1B15 (Aldo-Keto Reductase Family 1 Member B15) – белоккодирующий ген. Одним из процессов, связанных с этим геном, является метаболизм стероидных гормонов (кодируемый им белок действует, в основном, в качестве фермента, катализирующего восстановление андрогенов и эстрогенов).

Аннотация по базе dbsnp

dbsnp - база данных уже аннотированных полиморфизмов. Командва выдает 3 файла:

команда:perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out chr7.dbsnp -build hg19 -dbtype snp138 snp.avinput /nfs/srv/databases/annovar/humandb/

В результате было получено 3 файла:

chr7.dbsnp.hg19_snp138_dropped (SNP, имеющие идентификатор rs, т.е. аннотированные в базе данных dbSNP)

chr7.dbsnp.hg19_snp138_filtered (SNP, не имеющие rs)

chr7.dbsnp.log (содержит информацию о работе программы)

Было выяснено, что 25 SNP имеют аннотацию в упомянутой базе данных и 3 — не имеют. При этом все SNP без аннотации имеют низкое качество ридов (не более 8), в то время как для SNP с rs качество ридов сильно варьирует (от 1 до 71).

Аннотация по базе 1000 genomes

1000genomes - база данных полиморфизмов из 1000 геномов людей.

команда: perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out chr7_1000genomes -build hg19 -dbtype snp138 snp.avinput /nfs/srv/databases/annovar/humandb/

аналогично: опять 3 файла с таким же содержимым, но уже для базы 1000 genomes:

chr7_1000genomes.hg19_snp138_dropped (SNP, имеющие идентификатор rs, т.е. аннотированные в базе данных dbSNP)

chr7_1000genomes.hg19_snp138_filtered (SNP, не имеющие rs)

chr7_1000genomes.log (содержит информацию о работе программы)

В итоге в этой базе аннотировано 25 SNP, а не аннотированы те же самые 3, как и в предыдущем пункте.

Аннотация по базе GWAS

GWAS - база данных аннотированных полиморфизмов, которые ассоциированы с заболеваниями или особыми последствиями.

команда: perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -out chr7_gwas -build hg19 -dbtype gwasCatalog snp.avinput /nfs/srv/databases/annovar/humandb/

Было получено 2 файла:

chr7_gwas.hg19_gwasCatalog (с SNP, с известными клиническими значениями)

chr7_gwas.log (содержит информацию о работе программы)

В нашем случае было найдено 4 полиморфизма с аннотацией в GWAS:

Клинически значимые полиморфизмы	Координата	Замена	Тип замены	Качество ридов	Глубина покрытия
Name=Type 2 diabetes	100490077	G -> A	hom	221.999	34
Name=Bone mineral density	120969769	G -> A	hom	222.974	14
Name=Cortical thickness	120979089	C -> T	het	225.2	13
Name=Longevity	134250322	A -> C	het	225.009	68

1) синонимичная замена из гена ACHE: вызывает предрасположенность к диабету II типа

2) несинонимичная замена из гена WNT16: влияет на интенсивность минерализации костной ткани

3) несинонимичная замена из гена WNT16: влияет на толщину коркового слоя

4) Последняя замена: влияет на продолжительность жизни.

Аннотация по базе Clinvar

Clinvar - база данных о связях между полиморфизмами и их влиянием на организм человека (патогенные, неизвестные, отвечающие на лечение и пр.)

команда: perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out chr7_clincar -buildver hg19 -dbtype clinvar_20150629 snp.avinput /nfs/srv/databases/annovar/humandb/

В итоге было получено 3 файлa:

chr7_clincar.hg19_clinvar_20150629_dropped (содержит найденные аннотированные SNP)

chr7_clincar.hg19_clinvar_20150629_filtered (содержит все остальные SNP)

chr7_clincar.log

У нас первый файл оказался пустым, т.е. в данной базе ни один из изучаемых SNP не аннотирован.

Общая таблица найденных SNP

локализация	координата	координата	замена ->		тип замены	качество ридов	глубина покрытия
chr7	134261302	134261302	T	C	het	138.008	26
chr7	134246036	134246036	G	A	het	3.54557	1
chr7	100490077	100490077	G	A	hom	221.999	34
chr7	134254029	134254029	G	A	het	212.009	47
chr7	134252718	134252718	C	T	het	78.0075	6
chr7	134248045	134248045	A	G	hom	11.3429	1
chr7	134262441	134262441	T	C	het	225.009	71
chr7	134261674	134261674	C	T	het	225.009	54
chr7	100487721	100487721	G	T	hom	37.5136	4
chr7	134237294	134237294	A	G	hom	6.20226	1
chr7	134255326	134255326	C	T	hom	7.79993	1
chr7	134254427	134254427	G	A	hom	185.999	23
chr7	134252691	134252691	C	T	het	26.0177	5
chr7	134250322	134250322	A	C	het	225.009	68
chr7	120969769	120969769	G	A	hom	222.974	14
chr7	134264546	134264546	C	G	het	3.0136	3
chr7	134262747	134262747	G	C	het	68.0074	8
chr7	134261097	134261097	C	G	het	225.009	66
chr7	134260464	134260464	G	A	het	225.009	71
chr7	120965652	120965652	C	T	het	62.0073	8
chr7	134239978	134239978	C	A	het	29.0123	10
chr7	134260106	134260106	C	T	het	135.008	38
chr7	134264286	134264286	C	T	het	187.009	42
chr7	134259951	134259951	T	C	het	87.0076	8
chr7	134237239	134237239	T	C	hom	6.20226	1
chr7	120979089	120979089	C	T	het	225.2	13
chr7	134253269	134253269	C	T	hom	130.015	9
chr7	134259962	134259962	T	C	het	105.008	12

К семестрам