Поиск и описание полиморфизмов у пациента






Следите за обновлениями и дополнениями
Если Вы заметили опечатки, или ссылка испортилась, пожалуйста, напишите мне.

Чтобы ознакомиться с командами, использовавшимися при выполнении практикума, перейдите по ссылке.


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

0. Создание рабочей директории.

Для выполнения практикума была создвна рабочая директория /nfs/srv/databases/ngs/pavel-kravchenko, в которую были скопированы файл с ридами - chr3.fastq из директории /P/y14/term3/block4/SNP/reads и c хромосомой человеческого генома (сборка версии hg19) - chr3.fasta из директории /nfs/srv/databases/ngs/Human

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

Была установлена локальная версия программы FastQC и выполнен анализ качества ридов.
Команда: fastqc chr3.fastq

Программа выдаёт несколько массивов информации о интересующих нас ридах:
Basic Statistics - общая информация о числе ридов, методе секвенирования, GC составе и средней длине прочтения.
Так, для рассмтриваемого файла была показано наличие 20932 ридов, секвенирование проводилось технологией Sanger / Illumina 1.9, длина ридов оказалась равной 30-100 нуклеотидам, а GC состав был определён значением в 38%.

Per base sequence quality - график качества ридов, который показывает, какого качества оказались те или иные прочтения.
Per sequence quality scores - характеризует среднее качество прочтений. По графику можно понять, какое качество свойственно ридам, полученным после секвенирования.
Per base sequence content - данный график демонстрирует качество прочтения по нуклеотидам и может быть использован для оценки шума прибора и контроля качества прочтений. Также можно посмотреть на частоту встречаемости того или иного нуклеотида в конкретной позиции ридов.
Per sequence GC content - показывает, какова средняя доля GC пар в ридах.
Sequence Length Distribution - среднее распределение длины ридов в файле.
Kmer Content - график, отражающий положение повторяющихся последовательностей в ридах.

Отчёт по работе программы FastQC

Рисунок 1. Качество ридов до очистки

Чтобы разобраться в графике нужно помнить, как минимум, две вещи.
Первая - как соотносятся quality scores с вероятностью ошибок в позициях прочтений. Перевод из шкалы quality scores в вероятность ошибки показан на рисунке 2.
Вторая - что такое "ящик с усами".


Рисунок 2. Sequence quality scores
Рисунок 3. "Ящик с усами"

На рисунке 1 показаны те же ящики с усами. Красным, в перделах бокса, отмечена медиана, синим - среднее качество прочтений. Жёлтые прямоугольники - это интерквартирльные размахи (IQR на рисунке 3), то есть расстояние между нижней и верхней квартилью качества для каждой позиции прочтения. Зелёная зона имеет quality scores от 28 до 40, что является оптимальным для использования в генетическом анализе и интерпретировании результатов. Попадание в жёлтую и красную зоны говорит о более высоком уровне ошибок. На рисунке 4 показан классический ящик с усами.

Рисунок 4. Ящик с усами

С помощью программы Trimmomatic, установленной на kodomo, была произведена очистка чтений с помощью программы Trimmomatic. Для этого с конца каждого чтения были отрезаны нуклеотиды с качеством ниже 20 (TRAILING:20), и оставлены только чтения длиной не меньше 50 нуклеотидов (MINLEN:50).
Команда: java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr3.fastq chr3_outfile.fastq TRAILING:20 MINLEN:50
Программа получила на вход 20932 ридов и, после очистки, выдала 20570 ридов, что составило 98,27% от исходного их количества. Было отброшено 362 ридов.

Отчёт по работе программы FastQC

Рисунок 5. Качество ридов после очистки
Рисунок 6. Качество ридов до очистки
Рисунок 7. Качество ридов после очистки

Заметно, что после очистки среднее качество ридов заметно возросло - теперь практически все значения лежат в зелёной зоне графика 5, 7. За счёт обрезания низко качественных и коротких ридов длина усов заметно сократилась. Скорее всего, шумы, присутствовавшие в однонаправленных прочтениях до этого, были связаны с накоплением пошибок полимеразы ближе к 100ому нуклеотиду. Известно, что в протоколе Illumina могут возникать ошибки встраивания полимеразой меченных нуклеотидов, приводящие к понижению качества сигнала.

Рисунок 8. Распределение нуклеотидов по позициям ридов. Остаётся неизменным до и после очистки.
Рисунок 9. GC контент до чистки
Рисунок 10. GC контент после чистки приблизился к теоретически рассчитанному.


Рисунок 11. Длина прочтений до коррекции
Рисунок 12. Длина прочтений после коррекции уменьшилась, что и ожидалось.


Рисунок 13. Quality scores до коррекции
Рисунок 14. Quality scores после коррекции потеряли узкую модальность.


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

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

Чтобы картировать прочтения, нужно произвести ряд манипуляций с командной строкой UNIX:

  1. Сначала необходимо проиндексировать референсную последовательность.
  2. Затем построить выравнивание прочтений и референса в формате .sam.
  3. Сохранить вывод программы hisat2 в отдельный файл.



    Таблица 1. Использовавшиеся команды.
    Команда Результат
    export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5 Прописывает временный путь до исполняемой программы
    hisat2-build chr3.fasta chr3 Индексирует референсную последовательность, записывает результат в *.ht2 файлы
    hisat2 -x chr3 -U chr3_outfile.fastq --no-spliced-alignment --no-softclip > alignment.sam Строит выранивание прочтения и референса, записывает результат в .sam файл

    Программа hisat2 выдаласледующую информацию о выравниваниях. 

    20570 reads; of these:
  20570 (100.00%) were unpaired; of these:
    92 (0.45%) aligned 0 times
    20467 (99.50%) aligned exactly 1 time
    11 (0.05%) aligned >1 times
99.55% overall alignment rate

    Было поручено 20570 ридов, из которых 92 были выровнены 0 раз, а 20467 лишь один раз; 11 были выровнены более 1 раза.
    Был получен файл с выравниванием

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

    Для анализа выравнивания использовалась программа samtools, работающая с бинарными файлами.



    Таблица 2. Использовавшиеся команды.
    Команда Результат
    samtools view alignment.sam -b -o alignment.bam Переводит выравнивание в бинарный формат
    samtools sort alignment.bam -T out_sort.txt -o alignment_sorted.bam Индексирует референсную последовательность, записывает результат в *.ht2 файлы
    samtools index alignment_sorted.bam Индексирует отсортированный .bam файл
    samtools idxstats alignment_sorted.bam > out.txt Записывает откартировавшиеся прочтения

    Файл out.txt содержит информацию о картировании ридов: 

    chr3	198022430	20489	0
*	0	0	92


    Данная информация говорит о том, что на хромосому было картировано 20489 ридов, тогда как 92 ридов не было картировано. Длина хромосомы составила 198022430 нуклеотидов.

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

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

    Для поиска SNP и инделей были использованы команды, представленные в таблице 3.

    Таблица 3. Использовавшиеся команды.
    Команда Результат
    samtools mpileup -uf chr3.fasta alignment_sorted.bam > snp.bcf Создаёт файл с полиморфизмами в формате .bcf
    bcftools call -cv snp.bcf -o snp.vcf Создаёт файл со списком отличий между референсом и чтениями в формате .vcf

    В выдаче обнаружено более ста SNP, а также 11 инделей. Был сформирован файл, содержащий информацию о глубине прочтения каждого SNP.

    Команда: samtools depth alignment_sorted.bam >depth_sorted.tsv

    Cамыми интересными, на мой взгляд, оказались три, приведённые в таблице 4.

    Первые два были выбраны из-за высокой глубины покрытия и одинаковым качеством прочтений, а второе из SNP с худшими параметрами наугад.


    Таблица 4. Найденные SNP
    Кордината Tип полиморфизма: замена, вставка или делеция В референсной последовательности Что найдено в чтениях Глубина покрытия Качество прочтений
    41937051 Замена T C 256 225.009
    41925423 Замена T C 110 225.009
    41607701 Замена C G 36 137.008

    Визуализация SNPs...


    Рисунок 15. SNP в визуализаторе IGV выглядит так


    Рисунок 15. SNP 41937051
    Рисунок 16. SNP 41925423
    		 /
    Рисунок 17. SNP 41607701


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

    С помощью программы annovar были проаннотированы полученные snp. Использовались базы данных: refgene, dbsnp, 1000 genomes, GWAS, Clinvar. Для подготовки входных файлов из .vsf файла были удалены все индели, и для очищенного файла использовался скрипт convert2annovar.pl

    Команда: perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 /nfs/srv/databases/ngs/pavel-kravchenko/snp.vcf > /nfs/srv/databases/ngs/pavel-kravchenko/snp.avinput

    Результат:

    Finished reading 250 lines from VCF file
A total of 218 locus in VCF file passed QC threshold, representing 218 SNPs 
(121 transitions and 97 transversions) and 1 indels/substitutions
Finished writing 217 SNP genotypes (121 transitions and 96 transversions) 
and 1 indels/substitutions for 1 sample

    В файл было записано 217 SNP. Для аннотации SNP с помощью баз данных была использована программа annotate_variation.pl



    Аннотация по dbSNP

    Был проведён поиск rc SNP по базе данных dbSNP.

    Команда: /nfs/srv/databases/annovar/annotate_variation.pl -filter -out snp.rs -build hg19 -dbtype snp138 snp.avinput /nfs/srv/databases/annovar/humandb/

    Были получены следующие файлы, содержащие информацию о работе и результат работы программы. 

    snp.rs.hg19_snp_dropped
snp.rs.hg19_snp_filtered
snp.rs.hg19_snp138_dropped
snp.rs.hg19_snp138_filtered
snp.rs.invalid_input
snp.rs.log

    В файле snp.rs.hg19_snp138_dropped были записаны аннотированные SNP. 178 SNP имеют rs.

    С содержимым файла можно ознакомиться здесь

    Аннотация по RefGene

    Была выполнена аннотация по базе данных RefGene. По аннотации можно судить, в какой из функционально или топологически нагруженных участков попал найденный SNP.

    Команда: perl /nfs/srv/databases/annovar/annotate_variation.pl -out snp.refgene -build hg19 snp.avinput /nfs/srv/databases/annovar/humandb/

    Результаты работы программы можно найти в директории, а также посмотреть распределение SNP по группам в таблице 5.

    snp.refgene.exonic_variant_function
snp.refgene.log
snp.refgene.variant_function

    С содержимым файла можно ознакомиться здесь
    Для анализа файла может пригодиться таблица 5.

    Таблица 5. Расположение SNP по гуппам
    Зона Описание Число SNP
    intronic SNP в интронах 213
    exonic SNP в экзонах 13
    UTR3 SNP в 3'-некодирующей области 4
    UTR5 SNP в 5'-некодирующей области 2
    ncRNA_exonic SNP в транскрибируемой РНК, не имеющей аннотированного кодирующего участка 0

    Интересно, что первые два из выбранных мной SNP попали в белок кодирующие области. Это может помочь при дальнейшей аннотации.


    line147	synonymous SNV	ULK4:NM_017886:exon16:c.A1536G:p.Q512Q	chr3	41937051 41937051	T	C	het	225.009	256
line144	synonymous SNV	ULK4:NM_017886:exon17:c.A1599G:p.V533V	chr3	41925423 41925423	T	C	het	225.009	110

    Аннотация по 1000 genomes

    Была выполнена аннотация по базе данных 1000 genomes. Результат будет содержать частоту встречаемости найденных SNP.

    Команда: perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out snp.1000g -buildver hg19 -dbtype 1000g2014oct_all snp.avinput /nfs/srv/databases/annovar/humandb/

    Результаты работы программы можно найти в директории.

    snp.1000g.log
snp.1000g.hg19_ALL.sites.2014_10_dropped
snp.1000g.hg19_ALL.sites.2014_10_filtered

    С содержимым файла можно ознакомиться здесь
    Были посчитаны частоты встречаемости SNP. Максимальная составила 0.996805, а минимальная - 0.000599042.


    Аннотация по Gwas

    Была выполнена аннотация по базе данных Gwas. Она позволяет понять, с чем, из ранее описанного в литературе, ассоциированы SNP.

    Команда: perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -out snp.gwas -build hg19 -dbtype gwasCatalog snp.avinput /nfs/srv/databases/annovar/humandb/

    Результаты работы программы можно найти в директории. Было найдено четыре SNP, ассоциированных с уровнем кровяного давления, уровнем адипонектина и ростом человека. Использую поиск по GWAS Catalog, можно подробнее изучить все исследования, ассоциации и испытания, связанные с конкретным SNP. Поиск по базе осуществлялся с помощью rs идентификатора, полученного из базы данных 

    snp.gwas.hg19_gwasCatalog
snp.gwas.log

    С содержимым файла можно ознакомиться здесь


    Аннотация по Clinvar

    Была выполнена аннотация по базе данных Clinvar.

    Команда: perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out snp.clin -dbtype clinvar_20150629 -buildver hg19 snp.avinput /nfs/srv/databases/annovar/humandb/

    Результаты работы программы можно найти в директории.

    snp.clin.hg19_clinvar_20150629_dropped
snp.clin.hg19_clinvar_20150629_filtered
snp.clin.log

    Оказалось, что в Clinvar не аннотровано ни одного SNP. Как говорится, отрицательный результат - тоже результат.


    Ознакомиться с полной информацией о найденных SNP вы можете в таблице.

    Положение ридов на выравнивании вы можете увидеть на рисунке 18.

    Рисунок 18. Распределение ридов.



    Крайняя таблица. Использовавшиеся команды.
    Команда Результат
    fastqc chr3.fastq Анализ качества ридов
    java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr3.fastq chr3_outfile.fastq TRAILING:20 MINLEN:50 Очистка чтений
    export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5 Прописывает временный путь до исполняемой программы
    hisat2-build chr3.fasta chr3 Индексирует референсную последовательность, записывает результат в *.ht2 файлы
    hisat2 -x chr3 -U chr3_outfile.fastq --no-spliced-alignment --no-softclip > alignment.sam Строит выранивание прочтения и референса, записывает результат в .sam файл
    samtools view alignment.sam -b -o alignment.bam Переводит выравнивание в бинарный формат
    samtools sort alignment.bam -T out_sort.txt -o alignment_sorted.bam Индексирует референсную последовательность, записывает результат в *.ht2 файлы
    samtools index alignment_sorted.bam Индексирует отсортированный .bam файл
    samtools depth alignment_sorted.bam >depth_sorted.tsv Создаёт файл, содержащий информацию о глубине прочтения каждого SNP
    perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 /nfs/srv/databases/ngs/pavel-kravchenko/snp.vcf > /nfs/srv/databases/ngs/pavel-kravchenko/snp.avinput Подготовка файла для аннотирования
    /nfs/srv/databases/annovar/annotate_variation.pl -filter -out snp.rs -build hg19 -dbtype snp138 snp.avinput /nfs/srv/databases/annovar/humandb/ Поиск rc SNP по базе данных dbSNP
    perl /nfs/srv/databases/annovar/annotate_variation.pl -out snp.refgene -build hg19 snp.avinput /nfs/srv/databases/annovar/humandb/ Аннотация по базе данных RefGene
    perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out snp.1000g -buildver hg19 -dbtype 1000g2014oct_all snp.avinput /nfs/srv/databases/annovar/humandb/ Аннотация по базе данных 1000 genomes
    perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -out snp.gwas -build hg19 -dbtype gwasCatalog snp.avinput /nfs/srv/databases/annovar/humandb/ Аннотация по базе данных Gwas
    perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out snp.clin -dbtype clinvar_20150629 -buildver hg19 snp.avinput /nfs/srv/databases/annovar/humandb/ Аннотация по базе данных Clinvar

    Ссылки

    1. Распределение файлов с ридами
    2. FastQC
    3. Trimmomatic
    4. Руководство Trimmomatic
    5. samtools
    6. bcftools
    7. annotate_variation.pl
    8. convert2annovar.pl.
    9. Adiponectin