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

1. Подготовка

Сначала я создала требуемую директорию (/nfs/srv/databases/ngs/e.catsrina) и скопировала в нее нужные файлы (chr10.fasta и chr10.fastq).

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

С помощью программы FastQC (установлена на kodomo) был проведен анализ качества чтений (результат - по ссылке):

| fastqc chr10.fastq

Затем с помощью программы Trimmomatic была проведена очистка чтений: с конца каждого чтения были отрезаны нуклеотиды с качеством ниже 20 (TRAILING:20) и оставлены только чтения длиной не меньше 50 нуклеотидов (MINLEN:50). После чего опять была использована программа FastQC (результат - по ссылке):

| java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr10.fastq chr10_2.fastq TRAILING:20 MINLEN:50
| fastqc chr10_2.fastq

На рисунках 1 и 2 можно увидеть некоторые графических выдачи программы FastQC - Per base sequence quality и длины ридов - до и после чистки соответственно.

Рисунок 1. До чистки	Рисунок 2. После чистки
Рисунок 3. Длины ридов до чистки	Рисунок 4. Длины ридов после чистки

Как видно по этим рисункам, качество ридов (синяя линия) повысилось, особенно в конце ридов, разброс по значениям ("усы") - уменьшился (особенно в конце). Также теперь все позиции попадают желтой областью (все, что между нижним и верхним квартилями) в зеленую зону - в зону с хорошим качеством. Сильнее всего изменился конец, т. к. ридов длиной меньше 50 н было мало, поэтому их выброс повлиял на картину не так сильно, как отрезание некачественных - т. е. как раз концевых - нуклеотидов.

Естественно, изменились длина и число ридов: до очистки было 10666 ридов с длинами 38-100, после - 10526 ридов с длинами 50-100. Изменение распределения длин можно увидеть на рис. 3 и 4: увеличилась минимальная длина (на графике изменилось начало оси абсцисс), из-за уменьшения числа ридов увеличилась доля ридов с длиной от 96 до 100 н (их присутствие стало заметно на графике).

3. Картирование чтений.

Дальнейшие действия - картирование очищенных чтений по референсной последовательности с помощью программы Hisat2. Действия проводились в директории /nfs/srv/databases/ngs/e.caterina (команда cd /nfs/srv/databases/ngs/e.caterina).
Сначала была извлечена необходимая программа:

| export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5

| hisat2-build chr10.fasta chr10 > hisat.txt

| hisat2 -x chr10 -U chr10_2.fastq --no-spliced-alignment --no-softclip -S chr10.sam

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

С помощью пакета samtools был произведен анализ выравнивания.
Сначала файл с выравниванием ридов по референсной последовательности (chr10.sam) был переведен в бинарный формат (chr10.bam, перевод в бинарный формат обозначается -b):

| samtools view chr10.sam -o chr10.bam -b

| samtools sort chr10.bam -T chr10_temp.txt -o chr10_sort.bam

| samtools index chr10_sort.bam

Так же данные для анализа были выданы командой hisat2:

10526 reads; of these:
10526 (100.00%) were unpaired; of these:
131 (1.24%) aligned 0 times
10393 (98.74%) aligned exactly 1 time
2 (0.02%) aligned >1 times
98.76% overall alignment rate

Файл chr10.sam тоже содержит много информации: для каждого рида в нем написано (помимо его идентификатора, координаты в хромосоме и т. д. и некоторых других параметров):

• AS:i: - счет выравнивания
• XN:i: - число неопределенных нуклеотидов в референсной последовательности, которую покрывает данное выравнивание
• XM:i: - число несоответствий
• XO:i: - число открытых гэпов в выравнивании и в соответствующей референсной последовательности
• XG:i: - длина гэпов в выравнивании и в соответствующей референсной последовательности
• NM:i: - минимальное число однонуклеотидных замен, необхожимое для полного совпадения выравнивания и референса
• YT:Z:UU - индикатор того, что этот рид не в паре

Покрытие экзона

| samtools depth chr10_sort.bam > depth.txt

| samtools depth chr10_sort.bam -r chr10:63952845-63958202 > depth2.txt 

Рисунок 5. Окрестность нуклеотида с максимальным покрытием

Рисунок 6. Экзон 12

Рисунок 7.Экзон 6

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

Первым делом был создан файл с полиморфизмами (chr10_5.bcf):

samtools mpileup -uf chr10.fasta chr10_sort.bam -o chr10_5.bcf

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

	ID	REF		ALT		QUAL	INFO (тип - по умолчанию - замена, DP - глубина прочтения, прочее)			FORMAT	chr10_sort.bam
chr10	5766152	AGTAT		AGTATGTAT	92.4668	INDEL;DP=8;IDV=3;IMF=0.375;VDB=0.00746206;SGB=-0.511536;MQ0F=0;AF1=0.5;AC1=1;
							DP4=5,0,3,0;MQ=60;FQ=95.4767;PV4=1,1,1,1						GT:PL	0/1:130,0,194
chr10	5766337	G		A		225.009	DP=93;VDB=0.00111352;SGB=-0.693145;RPB=0.992146;MQB=1;MQSB=1;BQB=0.489669;MQ0F=0;
							AF1=0.5;AC1=1;DP4=38,11,35,6;MQ=60;FQ=225.007;PV4=0.423127,0.251369,1,1			GT:PL	0/1:255,0,255
chr10	5766572	Catatatatatat	CATatatatatatat	10.8393	INDEL;DP=5;IDV=1;IMF=0.2;VDB=0.267404;SGB=-0.379885;MQ0F=0;AF1=0.499996;AC1=1;
							DP4=0,2,0,1;MQ=60;FQ=13.6582;PV4=1,1,1,1						GT:PL	0/1:48,0,75

Файл chr10_5.vcf также был визуализирован с помощью программы IGV (см. ниже).

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

Аннотирование только что полученных SNP производилось с помощью программы annovar. Сначала вручную были удалены индели из файла chr10_5.vcf (получился файл chr10_6.vcf). Затем этот файл был переведен в нужный формат (perl необходим для получения доступа):

| perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 chr10_6.vcf -o chr10.avinput

1. По базе данных snp:

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

   Результат: файлы chr10.snp.hg19_snp138_dropped, chr10.snp.hg19_snp138_filtered, chr10.snp.log. В аннтотациях по типу -filter всегда файл _dropped включает в себя аннотированные snp, _filtered - не аннотированные (не имеющие аннотаций в этой базе данных), .log - сведения о команде. Было получено 54 аннотированных snp (из 57) - т. е. 54 snp имеют rs.

2. По базе данных 1000 genomes:

   | perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -dbtype 1000g2014oct_all -out chr10_1000g -buildver hg19 chr10.avinput /nfs/srv/databases/annovar/humandb/

   Результат: файлы chr10_1000g.hg19_ALL.sites.2014_10_dropped, chr10_1000g.hg19_ALL.sites.2014_10_filtered, chr10_1000g.log. Частоту встречаемости см в общей таблице ниже: частоты были найдены для всех snp, кроме трех, минимальная частота - 0.000798722, максимальная - 0.997604, средняя - 0.61, медиана - 0.72.

3. По базе данных Clinvar:

   | perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -dbtype clinvar_20150629 -out chr10_clinvar -buildver hg19 chr10.avinput /nfs/srv/databases/annovar/humandb/

   Результат: файлы chr10_clinvar.hg19_clinvar_20150629_dropped, chr10_clinvar.hg19_clinvar_20150629_filtered, chr10_clinvar.log. Все аннотированные по этой базе данных snp далжны были быть записаны в файл .hg19_clinvar_20150629_dropped, однако этот файл пуст. Что говорит о том, что никаких известных этой базе данных клинических последствий данные замены не несут.

4. По базе данных GWAS:

   | perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -dbtype gwasCatalog -out chr10_gwas -buildver hg19 chr10.avinput /nfs/srv/databases/annovar/humandb/

   Результат: Файлы chr10_gwas.hg19_gwasCatalog, chr10_gwas.log. Было найдено всего 3 замены (что составляет всего 5,3% от всех snp), ассоциированные с различными болезнями (остеосаркоме - раку кости, ревматоидному артриту и витилиго - исчезновению меланина с отдельных участков кожи) - см ниже (и в общей таблице). Интересно, что snp, приводящий к витилиго, назодится не в экзоне (по крайней данным refgen), в отличие от двух других snp:

   База данных	Проявление			Координата              REF	ALT	генотип	Качество рида	Глубина рида	
   gwasCatalog	Name=Osteosarcoma		chr10	5804531		G	A	het	225.009		70
   gwasCatalog	Name=Rheumatoid arthritis	chr10	63958112	T	C	het	225.009		75
   gwasCatalog	Name=Vitiligo			chr10	115481018	C	T	het	225.009		91
   

5. По базе данных refgene:

    
   | perl /nfs/srv/databases/annovar/annotate_variation.pl -out chr10_refg -build hg19 chr10.avinput /nfs/srv/databases/annovar/humandb/

   Результат: Файлы chr10_refg.exonic_variant_function, chr10_refg.variant_function, chr10_refg.log. В интронах содержится 45 замен, в экзонах 10 (причем из них 3 синонимичных и 10 несинонимичных), в UTR3 (нетранслируемый участок на 3'-конце зрелого транскрипта - после стоп-кодона) - 2. Информация о заменах в экзонах (включая получившуюся замену аминокислоты) сведена с таблицу в файле .exonic_variant_function:

   Тип			Ген (два ID):экзон:замена нукл:замена а. к.	Координата		REF	ALT	генотип Качество рида	Глубина рида
   nonsynonymous SNV	FAM208B:NM_017782:exon13:c.T1495G:p.C499G,	chr10	5781628		T	G	het	117.008		21
   nonsynonymous SNV	FAM208B:NM_017782:exon14:c.A2419G:p.I807V,	chr10	5784151		A	G	het	225.009		83
   nonsynonymous SNV	FAM208B:NM_017782:exon15:c.G3224C:p.R1075P,	chr10	5788608		G	C	het	108.008		14
   nonsynonymous SNV	FAM208B:NM_017782:exon15:c.T5036C:p.V1679A,	chr10	5790420		T	C	het	135.008		14
   nonsynonymous SNV	FAM208B:NM_017782:exon20:c.G7211A:p.S2404N,	chr10	5804531		G	A	het	225.009		70
   nonsynonymous SNV	RTKN2:NM_145307:exon12:c.A1385G:p.H462R,	chr10	63958112	T	C	het	225.009		75
   nonsynonymous SNV	CASP7:NM_001267057:exon2:c.T236G:p.I79S,
   			CASP7:NM_001227:exon2:c.T12G:p.D4E,
   			CASP7:NM_033338:exon3:c.T111G:p.D37E,
   			CASP7:NM_001267056:exon2:c.T12G:p.D4E,
   			CASP7:NM_033340:exon2:c.T12G:p.D4E,
   			CASP7:NM_033339:exon3:c.T12G:p.D4E,		chr10	115457264	T	G	het	225.009		71
   synonymous SNV		RTKN2:NM_145307:exon10:c.G1149A:p.K383K,	chr10	63964653	C	T	het	206.009		38
   synonymous SNV		RTKN2:NM_145307:exon6:c.G528A:p.V176V,		chr10	63995983	C	T	hom	221.999		154
   synonymous SNV		FAM208B:NM_017782:exon13:c.A1836T:p.T612T,	chr10	5781969		A	T	het	225.009		18
   

Ниже (в таблице 1) приведены все SNP с аннотациями их этих баз данных (базы Clinvar нет, потому что она не выдала ни одной аннотации).

Таблица 1. SNP и их аннотации

С помощью программы IGV был визуализирован участок, в котором есть индель и две замены (см рис. 5 и 6):

Рисунок 5.Визуализация

Рисунок 6.Визуализация

НАЗАД ➜

© <Рюмина Екатерина>, 2017