Ресеквенирование. Поиск полиморфизмов у человека
Таблица с командами
Программа | Что делает |
fastqc chr13.fastq | Анализ качества чтений |
java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr13.fastq chr13_aftertrim.fastq TRAILING:20; java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr13_aftertrim.fastq chr13_edit.fastq MINLEN:50 | Удаляет нуклеотиды с качеством меньше 20 с конца чтений, оставляет только чтения, чья длина больше 50 |
hisat2-build chr13.fasta chr13 | Индексирование референса (так быстрее потом картировать) |
hisat2 -x chr13 -U chr13_edit.fastq -S chr13.sam --nosoftclip --no-spliced-alignment &> hisat2_output | Кладёт прочтения на референс (запрещаем подрезать прочтения с концов и разрывать прочтения) |
samtools view chr13.sam -b -o chr13.bam | Преобразование полученного на предыдущем шаге sam - файла в бинарный bam - файл |
samtools sort chr13.bam chr13_sorted | Сортировка по координате прочтения (относительно референса) |
samtools index chr13_sorted.bam | Индексирование bam - файла |
samtools mpileup -u -g -f chr13.fasta -o chr13.bcf chr13_sorted.bam | Поиск отличий между тем, что мы отсеквенировали, и референсным геномом |
bcftools call -cv chr13.bcf -o chr13.vcf | Преобразование полученного на предыдущем шаге bcf - файла в vcf - файл (на него можно смотреть глазами) |
convert2annovar.pl -format vcf4 chr13.vcf -o chr13.avinput | Конвертирование vcf-файла в файл, с которым могут работать программа annovar |
annotate_variation.pl -out refg -build hg19 -dbtype refGene chr13.avinput /nfs/srv/databases/annovar/humandb/ | refgene - gene-based annotation |
annotate_variation.pl -filter -out snp -build hg19 -dbtype snp138 chr13.avinput /nfs/srv/databases/annovar/humandb.old/ | dbsnp - filter-based annotation |
annotate_variation.pl -filter -out 1000g -build hg19 -dbtype 1000g2014oct_all chr13.avinput /nfs/srv/databases/annovar/humandb.old/ | 1000 genomes - filter-based annotation |
annotate_variation.pl -regionanno -out gwas -build hg19 -dbtype gwasCatalog chr13.avinput /nfs/srv/databases/annovar/humandb.old/ | Gwas - region-based annotation |
annotate_variation.pl -filter -out clinvar -build hg19 -dbtype clinvar_20150629 chr13.avinput /nfs/srv/databases/annovar/humandb.old/ | Clinvar - filter-based annotation |
Часть 1: подготовка чтений
Выдача программы fastqc
До выполнения программы trimmomatic:
Количество чтений - 12155, длина - от 39 до 100.
После trimmomatic:
Количество чтений - 11933, длина - от 50 до 100. Видно, что качество чтений хуже всего в начале и в конце. Отрезая плохо прочтённые нуклеотиды с конца, мы укорачиваем все чтения, а затем убираем те из них, которые в результате этой процедуры стали слишком короткими. Благодаря этому, мы добились пристойного качества на концах прочтений и уменьшили разброс по качеству, начиная примерно с 60 нуклеотида (просто за счёт уменьшения количества длинных, некачественных в конце ридов).
Часть 2: картирование чтений
Выдача программы hisat2 выглядит следующим образом:
11933 reads; of these: 11933 (100.00%) were unpaired; of these: 117 (0.98%) aligned 0 times 11112 (93.12%) aligned exactly 1 time 704 (5.90%) aligned >1 times 99.02% overall alignment rate
Иными словами, из 119333 на хромосому однозначно откартировалось 11112 чтений. Ещё 704 откартировались неоднозначно и 117 вообще никуда не откартировались.
Часть 3: aнализ SNP
Рассмотрим несколько строк файла vcf.
Позиция Референс Что мы наблюдаем Качество Глубина покрытия chr13 25517089 TCCC TCC 217.468 103 chr13 20006620 C T 221.999 19 chr13 25518199 A G 5.46383 1
В верхней строке описана делеция, а в двух других - замены. Посдедний SNP приведён как иллюстрация ошибки: низкое качество и глубина в одно прочтение говорят в пользу того, что это не настоящий полиморфизм, а просто ошибка секвенирования, тогда как, например, делеция из первой строки, по всей видимости, действительно имеет место (причём отмечено, что она встречается примерно в половине прочтений, IDV=51;IMF=0.495146, что говорит о гетерозиготности в этом месте).
Всего в файле vcf 163 снипа (115 transitions, 48 transversions) и 15 инделей. Transition - мутация, при которой пуриновое основание переходит в пуриновое или пиримидиновое в пиримидиновое, transversion - наоборот.
На какие категории делит snp база данных refseq в annovar? Сколько snp у Вас попало в каждую группу?
exonic(9), intronic(90), intergenic(2), upstream(1), downstream(1), UTR5(1), ncRNA_exonic(12), ncRNA_intronic(61)
В какие гены попали Ваши snp?
TPTE2P1, COL4A1(исходя из выдачи refgene, учитываются только полиморфизмы, попавшие в экзоны)
К каким нуклеотидным и аминокислотным заменам привели snp?
6 синонимичных, 3 несинонимичных(refgene)
Сколько snp имеет rs?
140(по базе данных dbsnp)
Что Вы можете сказать о частоте найденных snp?
Частоту можно узнать из выдачи 1000genomes
Что Вы можете сказать о клинической аннотации snp?
В выдаче gwas:
Name=Obesity-related traits chr13 25533831 25533831 A C het 225.009 78 Name=Cardiac hypertrophy chr13 50080847 50080847 A G hom 221.999 45 Name=Arterial stiffness chr13 110818598 110818598 T G het 212.009 49
То есть найдены SNP, ассоциированные с ожирением, сердечной гипертрофией и артериальной ригидностью.
© Быкова Даша, 2018