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

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

Мне была дана 17 хромосома человека chr17.fasta и, соответственно, файл с ридами для той же хромосомы chr17.fastq (без адаптеров).

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

  java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr17.fastq chr17_out.fastq TRAILING:20 MINLEN:50

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

  fastqc *.fastq

где вместо * был вписан либо исходный файл chr17.fastq, либо файл с ридами после чистки chr17_out.fastq. В итоге было получено два архива, а также 2 .html файла с отчетами: chr17_fastqc.html и chr17_out_fastqc.html. На рис. 1-4 представлены основные характеристики для полученных файлов с чтениям.
Прим: конечно же, я сначала просмотрела характеристику чтений до очистки, чтобы понять, нужна ли очистка вообще, но обе характеристики приведены уже после очистки рядом друг с другом для удобства ☺

Как видно из рис. 1 и 2 после чистки действительно остались только риды с длиной больше 50 нуклеотидов. При этом было вырезано 178 чтений (1,6% от исходного количества), так что из 11046 чтений осталось только 10868.
Как видно из рис. 3 и 4 после чистки, показатели качества п.н. в ридах стали намного лучше. Таким образом, мы получили только те чтения, качество которых нас устраивает, и с которыми можно работать дальше.

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

• ✔Per sequence quality scores — распределение среднего качества среди ридов. В нашем случае все хорошо: пик на 38, т.е. качество у большинства ридов очень хорошее.
• ! Per base sequence quality — процент каждого нуклеотида в ридах. В рандомной библиотеке процент разных нуклеотидов для всех позиций, по идее, должен быть примерно одинаковым, так что линии на графике должны быть параллельными. Если каких-то нуклеотидов в последовательности больше, то это отразится в более высоко расположенной кривой, но параллельность все равно должна сохраняться. В нашем случае неровности относительно сильны, так что, возможно, это должно заставить задуматься о причинах.
• ✖Per sequence GC content — распределение GC-состава для всех ридов. Для нормальной рандомной библиотеки распределение GC-состава должно быть нормальным, т.е. центральный пик должен совпадать с пиком для GC-состава всего генома. Т.к. исходного генома нет, вместо него используется геном, построенный с помощью расчетов из имеющихся данных. В нашем случае этот модуль отмечен как failure, т.к. отклонение от нормального (гауссовского) распределения составляет более 30% ридов. Причиной может быть загрязненность библиотеки.
• ✔Per base N content — график, отражающий позиции, в которых был неопределен нуклеотид (N). В нашем случае все хорошо (для всех позиций 0 (< 5%)).
• ! Sequence Length Distribution — график, отображающий распределение длин ридов. В нашем случае модуль обозначен warning. Должен быть один пик для какой-то длины, а в нашем случае в библиотеке соотносимое число ридов с различающимися длинами (т.е. есть несколько вариантов преобладающих длин), что странно, т.к. обычно используются риды одной длины.
• ✔Sequence Duplication Levels
• ✔Overrepresented sequences
• ✔Adapter Content
• ✔Kmer Content

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

Для картирования чтений использовалась программа BWA. Сначала я проиндексировала референсную последовательность chr17.fasta:

  bwa index chr17.fasta

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

  bwa mem chr17.fasta chr17_out.fastq > align.sam

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

  samtools view align.sam -b -o align.bam

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

  samtools sort align.bam -T wtf.txt -o al_sort.bam

Параметр -T был задан для записи временных файлов в файл wtf.txt, т.к. иначе они выдавались в stdout. Далее полученный файл был проиндексирован:

  samtools index al_sort.bam

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

  samtools idxstats al_sort.bam > out.out

получив результат, записанный в файл out.out. Таким образом, на геном откартировалось 10867 ридов из 10868, т.е. не откартировался всего один. Что значат остальные цифры непонятно.

(*) Определение среднего покрытия для экзона
С помощью команды:

  samtools depth al_sort.bam > cover.out

был получен файл cover.out с вычисленными покрытиями для нуклеотидов. Я выбрала нуклеотид в позиции 44717428 (покрытие 105). Экзон, содержащий этот нуклеотид (5), был найден в GenomeBrowser. Он имеет границы 44,717,420-44,717,527 и принадлежит гену NSF.
Далее с помощью команды:

  samtools depth -r chr17:44717420-44717527 al_sort.bam > cover.exon

был получен файл cover.exon, содержащий покрытия нуклеотидов нашего экзона. В таблице есть полученные данные, а также вычисленное среднее покрытие, которое составило около 95.5 (много).

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

Сначала был создан файл с полиморфизмами в формате .bcf с помощью команды:

  samtools mpileup -uf chr17.fasta al_sort.bam -o snp.bcf

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

  bcftools call -cv snp.bcf -o snp.vcf

в результате чего был получен файл snp.vcf, содержащий список отличий между референсом и ридами в формате .vcf. Всего было найдено 52 SNP и 4 делеции. Качество SNP варьирует ~ в пределах от 3 до 225 (очень хорошее), среднее значение — 67. При этом для 56 находок качество >20 (хорошее). Максимальное покрытие составляет 63, а для целых 24 SNP покрытие составляет 1 рид [не совсем понимаю, какое покрытие мы условно считаем хорошим]. В табл. 1 приведено описание для трех полиморфизмов из этого файла.

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

Координата	Тип полиморфизма	В референсе	В ридах	Глубина покрытия	Качество ридов
44788310	Замена	G	A	74	221.999
79539193	Индель (делеция)	Ttttattta	Tttta	10	122.467
44833088	Индель (делеция)	ATTTTTTTTTTTTTTTTTTT	ATTTTTTTTTTTTTTTT, ATTTTTTTTTTTTTTTTT	101	5.78792

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

В этой части задания необходимо было с помощью программы annovar проаннотировать только полученные SNP (индели были вручную удалены в файле snp_only.vcf). Использовались базы данных: refgene, dbsnp, 1000 genomes, GWAS, Clinvar.

Сначала необходимо было из нашего файла snp_only.vcf получить файл, с которым может работать программа annovar. Для этого использовалась команда:

  perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 snp_only.vcf -outfile snp.avinput

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

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

  perl /nfs/srv/databases/annovar/annotate_variation.pl -out chr17.refgene -build hg19 snp.avinput /nfs/srv/databases/annovar/humandb/

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

Таблица 2. Общая информация о найденных SNP

intronic	exonic	UTR3	hom	het
44	3	5	34	18

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

Координата	Тип замены	Ген	В референсе -> в ридах	Замена аминокислоты	Глубина покрытия	Качество ридов	Изображение в IGV
62007498	синонимичная	CD79B	A -> G	C122C, C123C	63	221.999	тык
79589242	синонимичная	NPLOC4	G -> A	T53T	104	221.999	тык
79596811	синонимичная	NPLOC4	C -> T	P12P	48	221.999	тык

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

• CD79B (Cluster of Differentiation 79B) — ген, кодирующий белок Ig-β, входящий в состав белкового комплекса, образующего рецептор В-клетки, узнающий антиген¹;
• NPLOC4 (Nuclear protein localization protein 4 homolog) — ген, кодирующий белок NPL4 homolog), который вместе с двумя другими белками образует комплекс, связывающийся с убиквитинированными белками и обеспечивающий экспорт неправильно свернутых белков из ЭПС в цитоплазму, где они затем подвергаются деградации протеасомой². На NCBI также имеет второе название: ubiquitin recognition factor³.

(*) Визуализация полиморфизмов
С помощью программы IGV описанные полиморфизмы необходимо было визуализировать. Для этого в программе был открыт файл al_sort.bam. Затем в строке поиска были указаны позиции моих полиморфизмов. Так как изображения получились слишком большими, здесь я их выкладывать не буду; их можно открыть, нажав на соответствующую ссылку в таблице выше (откроются в новой вкладке).


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

  perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out chr17.dbsnp -build hg19 -dbtype snp138 snp.avinput 
  /nfs/srv/databases/annovar/humandb/

В результате было получено 3 файла: chr17.dbsnp.hg19_snp138_dropped (SNP, имеющие идентификатор rs, т.е. аннотированные в базе данных dbSNP), chr17.dbsnp.hg19_snp138_filtered (SNP, не имеющие rs) и такой же файл chr17.dbsnp.log, как в предыдущем случае.
Таким образом, 46 SNP имеют аннотацию в упомянутой базе данных, а 6 — не имеют. При этом все SNP без аннотации имеют низкое качество ридов, в то время как для SNP с rs качество ридов сильно варьирует.

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

  perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out chr17.1000 -build hg19 -dbtype snp138 snp.avinput 
  /nfs/srv/databases/annovar/humandb/

Тут все аналогично: опять три файла с таким же содержимым, но только для базы 1000 genomes. В итоге в этой базе аннотировано 46 SNP, а не аннотировано — 6, как и в предыдущем пункте. Частота варьирует от 0.000199681 (т.е. в крайней степени редко) до 0.979034 (очень часто). Примечательно, что SNP, входящие в экзоны (упоминались в аннотации по базе refgene) встречаются часто (частота > 0.5).

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

  perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -out chr17.gwas -build hg19 -dbtype gwasCatalog snp.avinput 
  /nfs/srv/databases/annovar/humandb/

Было получено 2 файла: chr17.gwas.log и chr17.gwas.hg19_gwasCatalog (содержит SNP, для которых известно какое-то клиническое значение). В нашем случае было найдено 3 полиморфизма с аннотацией в GWAS: первый вызывает предрасположенность к раку яичников и болезни Паркинсона, второй связан с ростом, а третий — с цветом глаз.

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

  perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out chr17.clincar -buildver hg19 -dbtype clinvar_20150629 snp.avinput 
  /nfs/srv/databases/annovar/humandb/

В итоге было получено 3 файлa. Не считая *.log, это файлы: chr17.clincar.hg19_clinvar_20150629_dropped (содержит найденные аннотированные SNP) и chr17.clincar.hg19_clinvar_20150629_filtered (содержит все остальные). В нашем случае первый файл оказался пустым, т.е. в данном БД ни один из изучаемых SNP не аннотирован.

Все найденные SNP (учитывались только файлы *dropped) можно найти в общей таблице Excel: pr13.xlsx.

Таблица 4. Сводная таблица со всеми командами

Команда	Что делает
java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr17.fastq chr17_out.fastq TRAILING:20 MINLEN:50	Принимает на вход файл с чтениями, в выходной файл записывает риды, длиной не менее 50 нуклеотидов (MINLEN:50), при этом с концов ридов отрезаны нуклеотиды с качеством ниже 20 (TRAILING:20)
bwa index chr17.fasta	Индексирует файл с референсной последовательностью
bwa mem chr17.fasta chr17_out.fastq > align.sam	Строит выравнивание ридов с референсной последовательностью и записывает их в файл .sam
samtools view align.sam -b -o align.bam	Переводит файл с выравниванием в формат .bam (-b) и записывает все в файл, указанный после -o
samtools sort align.bam -T wtf.txt -o al_sort.bam	Принимает на вход файл с выравниванием ридов с референсом и сортирует по координате в референсе начала рида. Выходной файл после -o, в файл после -T записывается все, что выдается в stdout
samtools index al_sort.bam	Индексирует файл с выравниванием
samtools idxstats al_sort.bam > out.out	Определяет, сколько ридов откартировалось на геном, записывая результаты в файл (out.out)
samtools depth al_sort.bam > cover.out	Вычисляет покрытия для нуклеотидов и записывает результат в файл (cover.out)
samtools depth -r chr17:44717420-44717527 al_sort.bam > cover.exon	Вычисляет покрытия для нуклеотидов из указанного региона (-r) и записывает в файл (cover.exon)
samtools mpileup -uf chr17.fasta al_sort.bam -o snp.bcf	Находит полиморфизмы и записывает их в файл (snp.bcf)
bcftools call -cv snp.bcf -o snp.vcf	Создает файл (после -o) с отличиями между референсом и рилами в формате .vsf
perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 snp_only.vcf -outfile snp.avinput	Создает из входного файла файл, ск оторым может работать программа annovar
perl /nfs/srv/databases/annovar/annotate_variation.pl -out chr17.* -build hg19 -dbtype *DATABASE* snp.avinput /nfs/srv/databases/annovar/humandb/	Создает файл с аннотацией SNP из входного файла (snp.avinout) по указанной базе данных (*DATABASE*).

Ссылки:

[1] CD79B molecule // NCBI. [URL]
[2] Q8TAT6 (NPL4_HUMAN) // UniProt. [URL] [3] NPL4 homolog, ubiquitin recognition factor // NCBI. [URL]