На главную

Практикум 11

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

Команда	Эффект от выполнения
fastqc chr21.fastq	Принимает чтения в формате fasctq и создает zip архив и html файл с их характеристиками (качество, длина, GC состав etc.)
java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr21.fastq res_trimm_21.fastq TRAILING:20 MINLEN:50	Принимает на вход чтения в формате fastq и удаляет с конца каждого из них нуклеотиды с качеством меньше 20, оставляя только чтения с длиной не менее 50 нуклеотидов. Phred 33 обозначает особенности записи чтений софтом секвенатора (для новых данных phred всегда 33). В результате работы программы создается файл в формате fastq с обработанными чтениями res_trimm_21.fastq
fastqc res_trimm_21.fastq	Принимает обработанные в программе trimmonatic чтения в формате fasctq и создает zip архив и html файл с их характеристиками (качество, длина, GC состав etc.)

При чистке с концов чтений были удалены нуклеотиды с качеством ниже 20, а затем из всех чтений были оставлены те, что имеют длину не менее 50 нуклеотидов. Таким образом были отсеяны части чтений с низким качеством, и среди оставшихся чтений были отсеяны короткие. Отрезать плохие нуклеотиды с конца важно, потому что при секвенировании Illumina с увеличением числа раундов секвенирования увеличивается рассинхронизация между отдельными молекулами внутри "пятна из клонов". На box-plot результаты фильтрации заметны как исчезновение заезжающих в красную зону качества (ниже 20) усов (крайних значений качества) на концах ридов.
До чистки - 8158 чтений, после чистки - 7858 чтений.

box plot распределения качества секвенирования нуклеотидов для всех ридов из неотфильтрованного набора

box plot распределения качества секвенирования нуклеотидов для всех ридов из отфильтрованного набора

На изображенных ниже графиках вдоль оси х отмечено среднее качество прочтения нуклеотида на рид а по оси у - количество ридов с таким средним качеством. После чистки удаляются части ридов с плохим качеством и нижний порог среднего качества прочтения на нуклеотид сильно увеличивается (с 2 до 48-49), а мода среднего качества увеличивается с 48 до 101. (фактически не может быть меньше 20)

Распределение числа ридов, обладающих определенным значением среднего качества для неочищенных ридов

Распределение числа ридов, обладающих определенным значением среднего качества для очищенных ридов

Часть 2 : Картирование чтений

Команда	Эффект от выполнения
hisat2-build chr21.fasta chr21_index_base	Программа принимает на вход последовательность хромосомы в формате fasta а затем создает проиндексированную базу для работы программ картировщиков ридов на геном.
hisat2 —no-spliced-alignment —no-softclip -x chr21_index_base -U res_trimm_21.fastq -S mapped_reads.sam	Программа принимает на вход файлы с проиндексированной последовательностью хромосомы и отфильтрованные риды в формате fastq, затем программа картирует риды на геном, запрещая разделять рид на части и картировать его на разные части генома (считая рид при этом местом соединения двух экзонов лежащих в геноме в разных местах) (--no-spliced-alignment) и запрещая подрезать риды с концов (--no-softclip). Результатом работы программы является файл в формате sam с описанием ридов и тех мест в геноме, куда они картировались.
samtools view -b -o mapped_reads.bam mapped_reads.sam	Программа принимает откартированные риды в формате sam в формат bam для дальнейшей обработки программами.
samtools sort -T ./tmp/sorted.tmp -o sorted_reads.bam -O bam mapped_reads.bam	Сортирует файл с картированными ридами по тому месту, куда они откартировались в геноме в формате bam и записывает в файл.
samtools index sorted_reads.bam	Индексирует отсортированнные риды и записывает в файл sorted_reads.bam.bai

46 ридов из 7858 не было откартировано, 7811 откартировались на геном один раз, 1 откартировался больше одного раза.
Выдача hisat2 после картирования чтений на геном содержит информацию для чтения с каждым ID: было ли чтение откартировано и произошло ли это один или более раз, координату места, куда откартировалось чтение, качество картирования, количество совпадений вставок и делеций относительно референса при картировании в формате CIGAR и другое.

Команда	Эффект от выполнения
samtools depth sorted_reads.bam > depth_bases.txt	Программа принимает на вход сортированные чтения в формате bam и выдает глубину покрытия для каждого нуклеотида хромосомы
Скрипт	При помощи скрипта на python (принимает на вход файл с глубиной покрытия по каждому нуклеотиду) быпо посчитано количество нуклеотидов обладающих определенным значением глубины покрытия и построена гистограмма

Распределение числа нуклеотидов имеющих определенную глубину покрытия (по оси х - глубина покрытия по у - число нуклеотидов)

Из гистограммы видно, что большее число нуклеотидов покрыто один раз, максимальная глубина покрытия - 177 попробуем найти тот участок генома, в котором он лежит и понадемся что он окажется экзоном.
Нуклеотид с таким покрытием имеет координату на 21 хромосоме версии сборки генома hg:19 16336354. Он попадает в экзон nuclear receptor interacting protein или NRIP1. Согласно браузеру ensembl на выбранном мной варианте транскрипта NRIP1-001 или ENST00000400199.1 этот экзон имеет координаты 16,333,556 - 16,340,847.

depth -r chr21:16333556-16340847 sorted_reads.bam > coverage_exon.txt	Программа считает покрытие чтениями для нуклеотидов внутри указанного региона
Скрипт (вместо depth_bases.txt на вход надо подавать покрытие экзона чтениями coverage_exon.txt)	Скрипт так же посчитает распределение покрытия нуклеотидов чтениями и посчитает среднее покрытие

Среднее покрытие экзона чтениями - 50,66

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

samtools mpileup -u -g -f chr21.fasta -o snp_mapped.bcf sorted_reads.bam	Программа принимает на вход референсную последовательность генома в формате fasta, отсортированные и откартированные на геном риды в формате bam и создает файл с полиморфизмами в несжатом формате bcf
bcftools call -v -c -o variants.vcf snp_mapped.bcf	Программа принимает на вход файл в формате bcf и создает файл в формате vcf с описанием полиморфизмов

Всего было получено 87 полиморфизмов, из них 81 - SNP (56 транзиций и 25 трансверсий) и 6 INDEL
Качество и покрытие полиморфизмов сильно варьируют. Так существуют "хорошие" полиморфизмы c качеством около 200 и плохие с качеством 3, покрытие находится в примерных границах от 1 до 75. Полиморфизмы с одновременно хорошим покрытием и качеством встречаются и довольно часто. Среднее качество 53, среднее покрытие ридами - 8. Средние посчитаны для всех вариантов и для snp и для indel.

Описание полиморфизмов

Кордината	Тип полиморфизма	Референс	Чтения	Глубина покрытия данного места	Качество чтений в данном месте
43824106	SNP (транзиция)	A	G	8	126,008
16336494	INDEl (polymerase slippage)	CTTTTTTTTTT	CTTTTTTTTT,CTTTTTTTT	94	27,4634
43846956	SNP (трансверсия)	T	G	12	150,002

Полиморфизм 1

Полиморфизм 2

Полиморфизм 3

Описание полиморфизмов

convert2annovar.pl variants.vcf -format vcf4 -outfile variants.avinput	Принимает на вход полиморфизмы, записанные в формате vcf, и записывает их в формате, пригодном для дальнейшей аннотации при помощи annovar
cat variants.avinput \| grep -v "-" > only_snp_variants.avinput	Команда удаляет индели из файла с полиморфизмами
annotate_variation.pl -filter -out /srv/databases/ngs/tinaferryman/filter_dbsnp -buildver hg19 -dbtype snp138 /srv/databases/ngs/tinaferryman/only_snp_variants.avinput /srv/databases/annovar/humandb.old/	Программа принимает файл с информацией о snp и ищет такие же snp (filter based annotation) в базе данных dbsnp, описывая их (присваивает идентификатор в базе данных dbsnp - rs) Файл с информацией о имеющихся в dbsnp полиморфизмах - filter_dbsnp.hg19_snp138_dropped. Информация о не прошедших фильтрацию snp хранится в файле filter_dbsnp.hg19_snp138_filtered.
annotate_variation.pl -geneanno -out /srv/databases/ngs/tinaferryman/genan_refgene -buildver hg19 -dbtype refGene /srv/databases/ngs/tinaferryman/only_snp_variants.avinput /srv/databases/annovar/humandb.old/	Программа ищет аналогичные snp в базе данных refGene (gene based annotation) и описывает их (локализацию: в экзоне, интроне или нетранслируемом регионе с одного из концов мРНК, синонимичность или несинонимичность для замен в экзоне)
annotate_variation.pl -regionanno -out /srv/databases/ngs/tinaferryman/reg_gwas -buildver hg19 -dbtype gwasCatalog /srv/databases/ngs/tinaferryman/only_snp_variants.avinput /srv/databases/annovar/humandb.old/	Программа аннотирует snp по базе данных gwas catalog (region based annotation) и записывает в файл reg_gwas.hg19_gwasCatalog. Аннотация заключается в наличии ассоциаций snp с каим либо фенотипами или заболеваниями
annotate_variation.pl -filter -out /srv/databases/ngs/tinaferryman/filter_clinvar -buildver hg19 -dbtype clinvar_20150629 /srv/databases/ngs/tinaferryman/only_snp_variants.avinput /srv/databases/annovar/humandb.old/	Программа аннотирует snp по базе данных clinvar (filter based annotation) и записывает в файл filter_clinvar.hg19_clinvar_20150629_dropped, но у меня этот файл пустой, что говорит об отсутствии этих snp в базе данных clinvar
annotate_variation.pl -filter -out /srv/databases/ngs/tinaferryman/filter_1000g -buildver hg19 -dbtype 1000g2014oct /srv/databases/ngs/tinaferryman/only_snp_variants.avinput /srv/databases/annovar/humandb.old/	Программа аннотирует snp по filter based annotation и базе данных 1000 genomes и записывает результат в файл filter_1000g.hg19_ALL.sites.2014_10_dropped

После gene based annotation snp по базе данных refGene кроме файла .log , содержащего информацию о запуске программы, создается два файла:genan_refgene.variant_function и genan_refgene.exonic_variant_function. Первый файл разделяет все snp по локализации: интронные варианты (68 штук), экзонные варианты (4 штуки), варианты находящиеся в нетранслируемых регионах (в моем случае в 3'- 9 штук) (intronic, exonic, 3'UTR). Во втором файле отдельно содержится информация об экзонных вариантах и они делятся на две группы: синонимичные (3 штуки) и несинонимичные (4 штуки).
Гены, в которых содержатся полиморфизмы: NRIP1, UBASH3A, AGPAT3,
Несинонимичная замена в экзоне единственная: A меняется на G и серин меняется на глицин (:c.A52G:p.S18G). Остальные замены в экзонах синонимичные. Ссылка на таблицу с заменами
В файле filter_dbsnp.hg19_snp138_dropped лежат 60 SNP имеющих rs а в файле filter_dbsnp.hg19_snp138_filtered лежат 21 SNP не имеющий rs

Сводная таблица по базам данных GWAS refGene dbsnp Создана при помощи скрипта

Частота snp полученная при filter based annotation по базе данных 1000 genomes в среднем довольно высокая (десятки процентов или даже единичная), один вариант имеет частоту 0.0131789, что является наименьшим значением. Таким образом все найденные snp довольно распространены в популяции

Клинической аннотации нет так как в clinvar не проаннотированы найденные мной snp с 21 хромосомы пациента.