Dzha_reseq

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

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

Даны следующие файлы:

Был проведен контроль качества чтений с помощью программы FastQC.

Данная программа установлена на kodomo. Она была запущена с помощью команды: 

	fastqc chr6.fastq
В итоге на выходе программа выдала архив (.zip) с отчетом [html] . (см. Рис.1)

Всего в файле 10289 последовательностей, их длина колеблется в диапазоне 33-100, процент содержания GC пар — 40.


align
Рис. 1. Качество нуклеотидов в последовательностях до очистки.

На горизонтальной оси отмечены позиции в риде (1-100). Вертикальная ось соответствует некоторому значению 'качества' прочтения. Жёлтый блок охватывает вторую и третью квартиль распределения качества для для каждой позиции. Красная линия - медиана, синяя - среднее значение. Черные 'уcы' обозначают границы 10 и 90 перцентилей.

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

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

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

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

Результаты работы:
	
	TrimmomaticSE: Started with arguments: -phred33 chr6.fastq chr6_clean.fastq TRA
	Automatically using 8 threads
	Input Reads: 10289 Surviving: 10123 (98,39%) Dropped: 166 (1,61%)
	TrimmomaticSE: Completed successfully

Оценка нуклеотидов после очистки — отчет html. (Рис.2)


align
Рис. 2. Качество нуклеотидов в последовательностях после очистки.

После обработки в ридах остались только те последовательности, в которых качество нуклеотидов >25.

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

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

	bwa index chr6.fasta

Очищенные риды были картированы на хромосому:

	bwa mem chr6.fasta chr6_clean.fastq > align.sam

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

	 samtools view -b align.sam > align.bam

Выравнивание чтений с референсом было отсортировано по координате в референсе начала чтения:

	 samtools sort align.bam align-sorted

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

	samtools index align-sorted.bam

Было подсчитано, сколько чтений откартировалось на геном:

	samtools idxstats align-sorted.bam > idxstats.out

В результате, получен файл idxstats.out; на геном успешно откартировалось 10123 ридов, не откартировалось — 1.

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

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

Cоздан файл с полиорфизмами в формате .bcf.

	samtools mpileup -uf chr17.fasta align-sorted.bam -o snp.bcf

Cоздан файл со списком отличий между референсом и чтениями в формате .vcf.

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

Был получен файл snp.vcf. 11 инделей, 83 полиморфизм.

Табл. 1. Описание полиморфизмов
Координата Тип полиморфизма Референс Рид Покрытие Качество ридов
106966857 Делеция cacacacacatacacacaca cacacacaca 58 217.468
107016838 Вставка CTTT CTTTTT 35 212.468
107009090 Замена G A 9 94.0077


Аннотация SNP.

Сначала из файла snp.vcf были удалены все индели -> snp_not_indels.vcf

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

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

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

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

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

База данных refgene делит снипы на следующие категории:


Распределение SNP по группам в файле:

Полиморфизмы попали в следующие гены:

Табл. 2. Описание полиморфизмов
Координата Замена Тип замены Ген Качество ридов
106974933 T > G het AIM1 81.013
106989212 G > A hom AIM1 221.999
106967833 G > A het AIM1 191.009

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

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

Получено три файла:

78 SNP имеют аннотацию в упомянутой базе данных, а 6 — не имеют.

Аннотация по базе GWAS: 
	perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out chr6.1000 -build hg19 -dbtype snp138 snp.avinput /nfs/srv/databases/annovar/humandb/

Аннотация по данной БД работает по такому же filter-based-принципу, как и в случае с dbsnp. Выходные файлы аналогичны (_dropped, _filered, .log). В этой базе аннотировано 46 SNP, а не аннотировано — 6

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

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

Было получено 2 файла: chr6.gwas.log и chr6.gwas.hg19_gwasCatalog (содержит SNP, для которых известно какое-то клиническое значение).

Название болезни Ген
Системная красная волчанка 138195723
Ревматоидный артирит 138195723
Коронарная болезнь OPRM1

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

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

В итоге было получено 3 файлa: chr6.clincar.hg19_clinvar_20150629_dropped (содержит найденные аннотированные SNP) и chr6.clincar.hg19_clinvar_20150629_filtered (содержит все остальные).