Ресеквенирование

Таблица с командами

Этап	Команда	Что делала
1	fastqc chr8.fastq	Составляет отчёт о качестве
2	java -jar trimmomatic.jar SE -phred33 chr8.fastq chr8_1.fastq TRAILING:20 MINLEN:50	Отрезает нуклеотиды с нужными параметрами
3	hisat2-build chr8.fasta chr8_ind	Индексирует референсную последовательность
4	hisat2 -x chr8_ind -U chr8_1.fastq --no-spliced-alignment --no-sp > align.sam	Строит выранивание референса и прочтения
5	samtools view align.sam -bo align.bam	Переводит в бинарный формат
6	samtools sort align.bam -o align_sorted.bam	Сортирует по координате
7	samtools index align_sorted.bam	Индексирует выравнивание
8	samtools mpileup -uf chr8.fasta align_sorted.bam -o SNP.bcf	Создает файл с полиморфизмами
9	bcftools call -cv SNP.bcf -o SNP.vcf	Создает список отличий
10	perl convert2annovar.pl -format vcf4 SNP_1.vcf > SNP_1.av	Конвертирует в формат, необходимый для annovar
11	perl annotate_variation.pl -out snp_1.out -build hg19 SNP_1.av /nfs/srv/databases/annovar/humandb/	refgene - gene-based annotation
12	perl annotate_variation.pl -filter -dbtype snp138 -out snp_2.out -build hg19 SNP_1.av /nfs/srv/databases/annovar/humandb/	dbsnp - filter-based annotation
13	perl annotate_variation.pl -filter -dbtype 1000g2014oct_all -out snp_3.out -buildver hg19 SNP_1.av /nfs/srv/databases/annovar/humandb/	1000 genomes - filter-based annotation
14	perl annotate_variation.pl -regionanno -dbtype gwasCatalog -out snp_4.out -build hg19 SNP_1.av /nfs/srv/databases/annovar/humandb/	Gwas - region-based annotation
15	perl annotate_variation.pl -filter -dbtype clinvar_20150629 -out snp_5.out -buildver hg19 SNP_1.av /nfs/srv/databases/annovar/humandb/	Clinvar - filter-based annotation

Подготовка

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

Здесь лежит отчёт FastQC после.

Заметим, что после чистки Trimmomatic было оставлено 8227 чтений из 8367, которые удовлетворяли нашим параметрам и имели хорошее качество.

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

Всего было найдено 8227 рида. Из них 22 рида не было картировано, 7851 рида откартировано один раз, 354 рида откартировано более одного раза.

При картировании были использованны параметры: --no-softclip – запрет подрезания чтений и --no-spliced-alignment – картирование без разрывов, так как в геномной последовательности разрывов быть не должно.

Таблица с анализом некоторых SNP

Координата	Тип	Референс	Чтения	Глубина	Качество
116424270	Инсерция	A	AC	30	55.456
76468309	Делеция	GTTTTTTTTTTTT	GTTTTTTTTTT,GTTTTTTTTTTT	76	115.457
116423490	Замена	A	T	40	179.999

В дополнение к таблице можно сказать, что было найдено 39 SNP и 2 инделя. Один из этих инделей не до конца определён.

Анализ данных из annovar

Если быстро пройтись по вопросам, то можно сказать, что refseq делит SNP на две группы: синонимические и несинонимические. То есть приводящие к изменению кодируемой аминокислоты, или нет. Было получено 3 синонимических и 2 несинонимических замены в экзонах. Мутации попадали в гены CLU, HNF4G и TRPS1. Ну и ещё одно одно деление на экзонные, интронные и UTR регионы. В первом пункте я рассмотрел экзонные, так как лишь они приводят к значимым изменениям.

38 из 39 SNP было найдено в rs. При этом удивительно , что наибольшая частота SNP в популяции составляет 1. То есть все отсеквенированные люди её имеют. Не очень понятно, почему на референсе буква иная. Наименьшая частота тоже довольно далека от средней, 0.00499201. Но наибольшее удивление мне доставила клиническая аннотация SNP. Благо, их довольно мало. Очевидно, что не все SNP приводят к заболеваниям, да и аннотировать их довоьно сложно. Интересно, что они не совпали с экзонными SNP. Например, замена 27466315 T -> C Приводит к возникновению риска болезни Альцгеймера, а замена C -> T в 76478768 позиции приводит к изменению уровней мочевой кислоты. Файл с выходом clinvar оказался пуст. Это можно связать с отсутствием данных по этим мутациям в базе.