Posted Tuesday, December 5, 2016 by Marina Gladkova. Renewed Saturday, December 24, 2016

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

В данном практикуме мной были использованы файлы chr8.fasta, chr8.fastq. Для работы была создана папка /nfs/srv/databases/ngs/gladmarine.


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


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


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

Основые функции следующие:
  • импорт данных из BAM, SAM или FastQ-файлов (любой вариант)
  • быстрый обзор с цельую выявления проблемных участков
  • сводные диаграммы и таблицы для быстрой оценки данных
  • экспорт результатов в HTML based permanent report
  • операции оффлайн
Модули для анализа:
  • Basic statistics содержит информацию по следующим параметрам - имя файла (filename), тип файла (file type), кодировка (encoding), общее количество последовательностей (total sequences), отфильтрованные последовательности (filtered sequences), длина последовательности с учетом наименьшей и наибольшей (sequence length), общий процент ГЦ-пар (%GC);
  • Per base sequence quality показывает разброс качества определения оснований в каждой позиции fastq-файла. Для каждой позиции создана диаграмма размаха (BoxWhisker type plot). [1] Центральная красная линия соответствует медиане. Желтый прямоугольник представляет интерквартильный разброс (25-75%). Концы усов показывают статистически значимые точки 10% и 90%. Синяя линия отражает математическое ожидание. По оси ординат отмечены качественные баллы. Цвет фона показывает качество определения: зеленый, оранжевый, красный соответственно хорошей точности, допустимой и низкой. Если хотя бы для одного основания нижний квартиль меньше 10 или медиана меньше 25, программа выдает предупреждение о возможной проблеме. Если нижний квартиль меньше 5 или медиана меньше 20, выдается сообщение об ошибке;
  • Per sequence quality scores позволяет увидеть, есть ли в файле поднабор последовательностей с общим низким качеством прочтения. Частой причиной может служить плохое качество изображения (например, край видимой области), но такие подмножества ридов составляют небольшой процент. Большое количество некачественных ридов говорит о неисправности прибора. Предупреждение появляется, когда среднее значение качества становится ниже 27 (вероятность ошибки 0.2%). Сообщение об ошибке выдается при значениях меньше 20 (вероятность ошибки 1%);
  • Per base sequence content отражает графически частоты появления каждого из оснований на каждой позиции. В случайной библиотеке можно ожидать, что различия во встречаемости от цикла к циклу будут малы, следовательно линии на графике должны быть параллельны. Относительное количество каждого из оснований должно отражать общее количество данных оснований в геноме, поэтому отклонения графиков друг от друга не должны быть очень велики. В случае, когда в какой-либо из позиций различия между встречаемостью A и T или G и С превышают 10%, выдается предупреждение, с 20% - ошибка;
  • Per sequence GC content показывает GC-содержание по всей длине каждой последовательности в файле и сравнивает его с модельным нормальным распределением (синий). В случайной библиотеке распределение содержания GC близко к нормальному, где центральный пик соответсвует общему содержанию GC в геноме. Если суммарные отклонения от нормального распределения составляют более 15% выдается предупреждение, если более 30% - сообщение об ошибке;
  • Per base N content. При невозможности определения основания секвенатором появляется символ N. График отражает процентное содержание таких оснований для каждой позиции. При содержании N >5% появляется предупреждение, при N >20% - ошибка;
  • Sequence length distribution показывает распределение длин анализируемых последовательностей. В норме у графика имеется только один пик, соответсвующий наиболее часто встречающейся длине фрагмента. Предупреждение появляется, когда не все фрагменты имеют одинаковую длину. Сообщение об ошибке появляется, если хотя бы один фрагмент имеет нулевую длину;
  • Duplicate sequences подсчитывает уровень дупликации каждой последовательности библиотеки. На графике показываны относительные количества последовательностей с разными уровнями дупликации. Синяя линия - распределение уровней дупликации по всему набору последовательностей, красная - распределение после дедупликации. В хороших библиотеках показатели для большинства последовательностей смещены в левую часть графика. Лишние дупликации по последовательности и загрязнения приводят к появлению пиков в правой части. В случае, когда неуникальные последовательности составляют более 20% программа выдает предупреждение. Ошибка появляется при их количестве более 50%;
  • Overrepresented sequences. В нормальной библиотеке с высокой пропускной способностью, содержатся разнообразные последовательности, ни одна из которых не составляет значительной части от целого. Последовательность с высокой встречаемостью либо может быть биологически значимой, либо свидетельствует о загрязнении библиотеки. В данном модуле перечисляются все последовательности, составляющие более 0.1% от целого. Каждую из таких последовательностей программа также прогоняет по базе наиболее распространенных загрязнителей и выдает наиболее близкое совпадение, если такое имеется. При наличии хотя бы одной последоваетльности, встречающейся чаще 0.1%, Вы получите предупреждение, при наличии последовательности, встречающейся чаще 1% - сообщение об ошибке;
  • Adapter content показывает суммарное процентное содержание частей библиотеки, в которых были найдены адаптеры. Как только последовательноять адаптера найдена в риде, она считается присутствующей до самого конца, поэтому по мере продвижения по последовательности процент увеличивается. Если такая последовательность присутствует более, чем в 5% ридов, выдается предупреждение, при более, чем 10% - ошибка;
  • Kmer content работает из предположения, что любой маленький фрагмент последовательности не должен иметь позиционных сдвигов при условии его присутствия в нормальной библиотеке. Биологическая причина сдвигов может быть разной, однако они должны одинаково наблюдаться во всех позициях. Данный модуль измеряет число 7-меров для каждой позиции в анализируемой библиотеке, затем использует биномиальное распределение для определения значимых отклонений от равномерного покрытия на всех позициях. Определяются все k-меры, для 6 из них представлен график распределения. Для ускорения работы программа анализирует только 2% библиотеки, далее экстраполируя результаты. Последовательности длиной больше 500 основнаний обрезаются до указанной границы. В случае, если биномиальный p-value <0.01, Вы получите сообщение об ошибке. При биномиальном p-value <105 - сообщение об ошибке. В результате работы программы Вы получите архив (.zip), который содержит отчет о программе в виде html-файла.


Я пользовалась FastQC, установленной на kodomo, но для удобства установила также версию с графическим интерфейсом (Java-application).

Использованная командаРезультат
fastqc chr8.fastq
chr8_fastqc.zip
fastqc_report.html

По результатам работы программы можно сказать, что качество последовательности неплохое, за исключением 3 критериев, которые нужно проверить посредством очистки: Per base sequence content, Per sequence GC content, Sequence length distribution.

2. Очистка чтений.


Trimmomatic позволяет обрезать последовательности адаптеров Illumina (спаренных и неспаренных концевых участков). Выбор режима тримминга и ассоциированных параметров осуществляется в командной строке.

Доступные команды:
  • ILLUMINACLIP: вырезает адаптер и другие Illumina-специфичные последовательности рида
  • SLIDINGWINDOW: тримминг за счет сдвигающегося окна, вырезает нуклеотиды как только среднее качество области меньше определенной границы
  • LEADING: отрезает основания от начала рида, если качество меньше граничного балла
  • TRAILING: отрезает основания от конца рида, если качество меньше граничного балла
  • CROP: обрезает рид до определенной длины
  • HEADCROP: обрезает специфическое число оснований от начала до конца рида
  • MINLEN: отбрасывает рид, если его длина меньше заранее установленной
  • TOPHRED33: конвертирует баллы качества в Phred-33
  • TOPHRED64: конвертирует баллы качества в Phred-64

Адаптеры в исследуемой последовательности были удалены заранее. Чтобы обрезать с конца чтений нуклеотиды с качеством менее 20, был указан параметр TRAILING:20. С помощью MINLEN:50 были исключены чтения длиной меньше 50 нуклеотидов. В результате работы программы были удалены 140 ридов из всего 8367, таким образом число чтений после чистки составило 8227.

Использованная командаРезультат
java -jar /usr/share/java/trimmomatic.jar SE -phred33 
chr8.fastq chr8_trimmomatic.fastq TRAILING:20 MINLEN:5
chr8_trimmomatic.fastq

Аналогично полученный файл был проанилизирован программой FastQC.

Использованная командаРезультат
fastqc chr8_trimmomatic.fastq
chr8_trimmomatic_fastqc.html


Сравнение результатов до и после работы с Trimmomatic
Per base sequence quality
До очистки чтений
После очистки чтений
Можно сказать, что Per base quality было достаточно хорошим и до чистки. После работы Trimmomatic не осталось усов, выходящих за пределы оранжевой области. Почти не изменилось значение медианы, и немного увеличилось математическое ожидание. Trimmomatic отбросил чтения низкого качества на конце.
Per sequence GC content
До очистки чтений
После очистки чтений
Диаграммы данного модуля были мне интересны тем, что после очистки чтений изменился его "статус". В первом случае возникает предупреждение (!), то есть отклонения составляют более 15%, что хорошо заметно на графике. График не имеет гладкого пика и даже немного сдвинут вправо относительно графика нормального распределения. После работы Trimmomatic Warning не появляется (✓), мы видим, что красный график хорошо накладывается на синий. Правая часть полученного графика становится более гладкой. Наблюдаемые расширенные пики могут свидетельствовать о присутсвии различных загрязнителей.



Часть II. Картирование чтений.


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


Burrows-Wheeler aligner (BWA) - пакет программного обеспечения для картирования низкодивергирующих последовательностей при сравнении с большими референсными геномами (например, человеческим). Состоит из 3 алгоритмов: BWA-backtrack, BWA-SW и BWA-MEM. Первый разработан для анализа чтений Illumina длиной до 100 п.о., в то время как остальные два оптимизированы для работы с более длинными последовательностями от 70 п.о. до миллиона. BWA-MEM и BWA-SW обладают сходными характеристиками (long-read support и split alignment), но BWA-MEM, самый новый, соответственно наиболее быстрый и точный, рекомендуется для анализа высококачественных последовательностей. BWA-MEM более хорош в использовании при исследовании чтений Illumina 70-100 п.о. (нежели BWA-backtrack). Программа BWA использует преобразование Барроуза — Уилера. [2]

Используя программу BWA, установленную на kodomo, я откартировала чтения при помощи нижеуказанных команд:

Использованная командаОписаниеРезультат
bwa index chr8.fasta
Индексирование референсной последовательностиИндексированный chr8.fasta
bwa mem chr8.fasta
chr8_trimmomatic.fastq

> chr8.sam
Выравнивание очищенных прочтений и референса с использованием алгоритма BWA-MEMchr8.sam


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

Для выполнения задания использовался пакет Samtools. [3]

Использованная командаОписаниеРезультат
samtools view chr8.sam
-b -o chr8.bam
Перевод выравнивания чтений с референсом в бинарный формат .bam. Параметр -b меняет начальный формат на .bam, -o направляет output в файлchr8.bam
samtools sort chr8.bam -T
temporary.txt -o chr8_sort.bam
Сортировка выравнивания чтений с референсом (полученный bam-файл) по координате в референсе начала чтенияchr8_sort.bam
samtools index chr8_sort.bam
Индексирование отсортированного bam-файлаИндексированный chr8_sort.bam
samtools idxstats chr8_sort.bam
> chr8_result.out
Подсчет числа чтений, откартированных на геном chr8_result.out



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


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

При выполнении использовались также программы пакета BCFtools. [4]

Использованная командаОписаниеРезультат
samtools mpileup -uf chr8.fasta
chr8_sort.bam -o chr8_snp.bcf
Создание файла с полиморфизмами в формате .bcfchr8_snp.bcf
bcftools call -cv chr8_snp.bcf
-o chr8_snp.vcf
Создание файла со списком отличий между референсом и чтениями в формате .vcfchr8_snp.vcf

В итоговом файле chr8_snp.vcf описано 100 полиморфизмов. Встречаются как индели, так и замены. Примерно 2/5 от общего числа полиморфизмов имеют хорошее качество прочтения. Глубина покрытия превышает показатель в 50 достаточно редко.

Примеры полиморфизмов
КоординатаТип полиморфизмаРеференсРидКачество прочтения на участке Глубина покрытия на участке
116654469ИнсерцияCCATA97.43183
27454785ДелецияTAATGAATAA113.4677
116423424ЗаменаTC221.99949


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

Работая с программой ANNOVAR, нужно было аннотировать полученные ранее SNP. [5] Для этого из файла была удалена информация об инделях, к новому файлу (chr8_snp_noindl.vcf) я применила скрипт convert2annovar.pl.

Использованная командаОписаниеРезультат
perl /nfs/srv/databases/annovar/convert2annovar.pl

-format vcf4 chr8_snp_noindl.vcf

-outfile chr8_snp.avinput
Подготовка файла (перевод в формат ANNOVAR)chr8_snp.avinput

Были произведены аннотации SNP по следующим базам данных - refgene, dbsnp, 1000 genomes, GWAS и Clinvar посредством скрипта annotate_variation.pl. В ANNOVAR существуют 3 типа аннотаций, основанных на:
1) gene-based annotation - генной разметке;
2) region-based annotation - разметке других регионов генома;
3) filter-based annotation - фильтрации.
*Для каждой базы указан тип использованной аннотации.


1. Аннотация по базе refgene (gene-based).

Использованная командаОписаниеРезультат
perl /nfs/srv/databases/annovar/
annotate_variation.pl
  
-geneanno -dbtype refGene

-buildver hg19

chr8_snp.avinput -outfile

chr8_refgene /nfs/srv/databases/
annovar/humandb/
Аннотация файлов по генной разметкеchr8_refgene.exonic_variant_function - описание всех полиморфизмов;
chr8_refgene.variant_function - описание полиморфизмов в экзонах;
chr8_refgene.log - отчет о работе программы

В первой колонке файла *.variant_annotation отображены категории, указывающие местоположение полиморфизма (в скобках указано число снипов данной категории, встретившихся в моем файле):
  • exonic - полиморфизм, попадающий в экзон полностью или частично (5);
  • splicing - полиморфизм в пределах n п.о. от границы сплайсинга (n по умолчанию равно 2);
  • ncRNA - полиморфизм входит в транскрипт некодирующих РНК;
  • UTR5 - полиморфизм входит в 5′-НТО;
  • UTR3 - полиморфизм входит в 3′-НТО (13);
  • intronic - интронный полиморфизм (58);
  • downstream - полиморфизм на расстоянии 1-kb downstream от сайта окончания транскрипции (граница может быть изменена);
  • upstream - полиморфизм на расстоянии 1-kb upstream от сайта начала транскрипции (граница может быть изменена);
  • intergenic - полиморфизм межгенной области (18).
Таким образом, всего 94 полиморфизма, 23 из них - гетерозиготные, 71 - гетерозиготные. Если полиморфизм в одной из категорий: exonic, intronic или ncRNA, то во второй колонке указывается имя соответсвующего гена (несколько генов разделяются запятой). Иначе указаны два соседних гена и расстояние между ними. Больше всего замен произошло в интронных зонах, так как они не влияют на итоговые продукты.

Таблица экзонных полиморфизмов

№КоординатыГенНуклеотидная заменаТип заменыКачество чтений Глубина покрытияАK замена
127462481CLU [6]A > G (синонимичная)het80.0075 16H263H
276452313HNF4G [7]G > A (несинонимичная)hom221.999 30S29N
376468228HNF4GG > A (синонимичная)het225.00999L209L
476468282HNF4GG > A (несинонимичная)hom221.999104 M227I
5116631902TRPS1 [8]C > A (синонимичная)hom221.99935 P141P, P134P, P132P


2. Аннотация по базе dbsnp (filter-based).

Использованная командаОписаниеРезультат
perl /nfs/srv/databases/annovar/
annotate_variation.pl

-filter -dbtype snp138 -buildver

hg19 chr8_snp.avinput -outfile

chr8_snp /nfs/srv/databases/
annovar/humandb/
Аннотация файлов по фильтрации chr8_snp.hg19_snp138_dropped - полиморфизмы, имеющие идентификатор rs в базе данных dbsnp;
chr8_snp.hg19_snp138_filtered - оставшиеся полиморфизмы (не имеют описания rs);
chr8_snp.log - отчет о работе программы

76 полиморфизмов имеют rs-идентификатор, соответственно 18 - не имеют. Интересно отметить, что максимальное покрытие в случае snp без rs составило 5, у большинства этот показатель вообще составляет 1 или 2, в отличие от snp из файла *_dropped.


3. Аннотация по базе 1000 genomes (filter-based).

Использованная командаОписаниеРезультат
perl /nfs/srv/
databases/annovar/
annotate_variation.pl

-filter -dbtype
1000g2014oct_all

-buildver hg19
chr8_snp.avinput

-outfile
chr8_1000gnm

/nfs/srv/databases/
annovar/humandb/
Аннотация файлов по фильтрации chr8_1000gnm.hg19_ALL.sites.2014_10_dropped - полиморфизмы, имеющие идентификатор rs в базе данных 1000 genomes, кроме этого указана частота;
chr8_1000gnm,hg19_ALL.sites.2014_10_filtered - оставшиеся полиморфизмы (не имеют описания rs);
chr8_1000gnm.log - отчет о работе программы

По-прежнему, 76 полиморфизмов имеют указанный идентификатор, частота варьирует от примерно 0,005 до 1.


4. Аннотация по базе Gwas (region-based).

Использованная командаОписаниеРезультат
perl /nfs/srv/databases/
annovar/annotate_variation.pl

-regionanno -dbtype

gwasCatalog -buildver hg19

chr8_snp.avinput -outfile

chr8_gwas /nfs/srv/databases/
annovar/humandb/
Аннотация файлов по разметке регионов генома chr8_gwas.hg19_gwasCatalog - полиморфизмы, имеющие описанное клиническое значение;
chr8_1000gnm.log - отчет о работе программы


По данным файла можно видеть, что клиническое значение аннотировано в четырёх snp (это болезнь Альцгеймера, регуляция уровня мочевой кислоты, HDL ("хороший") cholesterol), причем среди них нет ни одного экзонного. Следовательно, можно сделать вывод о сложности процессов регуляции и экспрессии генов с данным клиническим значением.


5. Аннотация по базе Clinvar (filter-based).

Использованная командаОписаниеРезультат
perl /nfs/srv/databases/
annovar/
annotate_variation.pl

-filter -dbtype
clinvar_20150629

-buildver hg19
chr8_snp.avinput

-outfile chr8_clinvar
/nfs/srv/databases/
annovar/humandb/
Аннотация файлов по фильтрации chr8_clinvar.hg19_clinvar_20150629_dropped - полиморфизмы, имеющие аннотацию в базе Clinvar;
chr8_clinvar.hg19_clinvar_20150629_filtered - полиморфизмы, не имеющие аннотацию в базе Clinvar;
chr8_clinvar.log - отчет о работе программы


По результатам работы последней программы, не было найдено аннотаций в базе Clinvar (файл *_dropped оказался пустым).



Источники