Чтения экзома, кртирующиеся на участок хромосомы человека - chr2.fastq
Хромосомы человеческого генома (сборка версии hg19) - chr2.fasta
С помощью программы FastQC был проведён контроль качества чтений. Команда:
fastqc chr2.fastq |
Получив на вход файл в формате .fastq (данный формат характерен для файлов с данными после NGS), программа выдала архив (.zip) с отчетом по работе в виде текстового файла и рисунков, в котором помимо прочей информации находился файл chr2_fastqc.html, содержащий визуализацию качества чтений.
Используя программу Trimmomatic, производилась очистка чтений. Команда:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr2.fastq chr2_clean.fastq TRAILING:20 MINLEN:50 |
Здесь, параметр TRAILING отвечает за удаление нуклеотидов с конца каждого чтения, если их quality score меньше, чем 20, а параметр MINLEN оставляет риды длиной 50 и более, а риды с меньшей длиной удаляются. В результате число прочтений уменьшилось с 10410 до 10191, то есть на 2,1%. Далее, с помощью FastQC был вновь проведен анализ качества чтений и получен архив с файлом - chr2_clean_fastqc.html. Команда:
fastqc chr2_clean.fastq |
В Таблице 1 приведены различия по чтениям до и после очистки.
Таблица 1. Характеристика ридов до и после очистки. | ||
Категория |
До очистки |
После очистки |
Число ридов |
10410 |
10191 |
Длина ридов |
40-100 |
50-100 |
%GC |
42 |
42 |
На Рис. 1-2 представлена информация, полученная с помощью программы FastQC, а именно - характеристика качества каждого нуклеотида в ридах. По абсциссе расположены позиции в ридах, а по ординате - Quality score, причем данные представлены в виде диаграмм boxplot, где границами желтых столбцов служат первый и третий квартили, красная линия - медиана, синяя линия - среднее значение, а отходящие от столбцов "усы" демонстрируют размах между 10-м и 90-м процентилями. Помимо этого само поле графика покрашено в три цвета: нахождение диаграмм в зеленой плотности говорит о том, что качество данной позиции хорошее и она достоверна (выше 28), в оранжевой - среднее (от 20 до 28), в красной - плохое, таким последовательностям верить нельзя (от 0 до 20).
На Рис. 1 можно учидеть, что в 14 позициях "усы" диаграмм попадают в оранжевую область, а в трех уходят в красную, такой результат говорит о падении качества чтения последних нуклеотидов. Однако, в целом, большая часть позиций обладает высоким качеством прочтения. | |
Рис. 1. Характеристика качества последовательностей ридов до очистки. | |
Как и следовало ожидать, после очистки исчезли концевые нуклеотиды, обладающие quality score < 20, а следовательно, на графике это отразилось "укорачиванием усов" нуклеотидов, особенно значительно это для последних позиций. Теперь лишь в последней позиции диаграмма заходит в оранжевую область, а такую совокупность ридов можно достоверно использовать дальше. | |
Рис. 2. Характеристика качества последовательностей ридов после очистки. |
Помимо этого, как уже отмечалось выше, программа FastQC выдает различные графики, ниже на Рис. 3-4 представлена характеристика распределения quality score среди чтений.
Рис. 3. Распределение среднего качества до очистки. |
Рис. 4. Распределение среднего качества после очистки. |
Как видно из Рис. 3, до очистки среди всех ридов присутствуют чтения с качеством ниже 20, однако большая часть обладает качеством выше 30, что является очень хорошим показателем. |
После очистки можно увидеть, что среди всех ридов отсутствуют рида качеством ниже 20, а превалирующая часть чтений всё так же обладает средним качеством выше 30. |
С помощью программы Hisat2 было осуществлено картирование чтений. Сначала была проиндексирована референсная последовательность с помощью команды:
hisat2-build chr2.fasta chr2 |
В результате чего были получены 8 файлов в формате .ht2. Затем построено выравнивание очищенных чтений и референса в формате .sam с помощью команды, где параметр --no-softclip запрещает подрезать чтения, а --no-spliced-alignment требует проводить картирование без разрывов:
hisat2 -x chr2 -U chr2_clean.fastq --no-spliced-alignment --no-softclip > align_clean.sam |
В результате выдача сохранена в файл align_clean.sam.
С помощью команды, приведенной ниже, выравнивание из формата .sam было переведено в бинарный формат .bam:
samtools view align_clean.sam -bo align_clean.bam |
Для этого использовалась команда view из пакета samtools. Затем выравнивание чтений с референсом было отсортировано по координате в референсе начала чтения с помощью следующей команды:
samtools sort align_clean.bam -T 1.txt -o chr2_sort.bam |
И далее отсортированный файл chr2_sort.bam был проиндексирован командой:
samtools index chr2_sort.bam |
После сего получен файл chr2_sort.bam.bai. Затем выяснялось, сколько чтений откартировано на геном, эта информация выдавалась сразу после работы программы Hisat2:
10191 reads; of these: 10191 (100.00%) were unpaired; of these: 48 (0.47%) aligned 0 times 10139 (99.49%) aligned exactly 1 time 4 (0.04%) aligned >1 times 99.53% overall alignment rate
Из которой видно, что 48 чтений не было совсем откартировано, 10139 откартировалось ровно 1 раз, а 4 - больше, чем один раз. В то же время эту информацию можно получить с помощью команды:
samtools idxstats chr2_sort.bam > out.txt |
В результате в файле out.txt была записана следующая информация: длина референсной последовательности - 243199373 и число откартированных чтений, однако никаких пояснений не даётся.
Для выполнения данного задания с помощью команды, приведенной ниже, было вычислено покрытие для каждого из нуклеотидов:
samtools depth chr2_sort.bam > depth.txt |
Затем в полученном файле depth.txt был выбран нуклеотид с покрытием 206 в позиции 234189668. С помощью GenomeBrowser (версия генома hg19) определены координаты экзона, содержащего выбранный нуклеотид: 234160217-234204320. После этого была перезапущена команда samtools depth с указанием координат экзона:
samtools depth -r chr2:234160217-234204320 chr2_sort.bam > depth_ex.txt |
И получен файл depth_ex.txt, в котором обозначалось вычисленное покрытие для нуклеотидов из экзона, а с помощью программы Excel вычислялось среднее покрытие чтениями, оказавшееся равным 31,1568.
samtools mpileup -uf chr2.fasta chr2_sort.bam > snp.bcf |
Был создан файл с полиморфизмами в формате .bcf. Затем командой:
bcftools call -cv snp.bcf > snp.vcf |
Получен файл со списком отличий между референсом и чтениями в формате .vcf. В файле snp.vcf по результатам работы программы bcftools call было найдено 79 полиморфизмов, из которых 7 инделей и 72 замены. В Таблице 2 собрана информация о 3 полиморфизмах.
Таблица 2. Описание полиморфизмов. | |||||
Координата |
Тип полиморфизма |
Референс |
Чтения |
Глубина покрытия |
Качество чтений |
234181593 |
Замена |
G |
C |
28 |
134.008 |
234202274 |
Вставка |
TCC |
TCCC |
1 |
3.80767 |
238426328 |
Делеция |
gtggatggatggatggatggatggatggatggat |
gtggatggatggatggatggatggatggat |
1 |
22.4955 |
С помощью программы IGV - Integrative Genomics Viewer, были визуализированы полиморфизмы, приведённые в Таблице 2. На вход программа получает отсортированный файл с выравниванием в формате .bam. Рис. 5-7 демонстрируют результат работы этой программы, для каждого полиморфизма программа выдает краткую информацию. Верхней линией на рисунках показана картированная хромосома, что позволяет довольно легко найти нужный участок.
Рис. 5. Замена 234181593. |
Рис. 6. Вставка 234202274. |
Рис. 7. Делеция 238426328. |
Программа IGV выдает информацию о том, в каком числе ридов встречается тот или иной нуклеотид, например, в данном случае видно, что цитозин встречается в 9 из 28 чтений, а гуанин в 19 из 28, в референсной последовательности гуанин. |
В данном примере произошедшая вставка описана в пояснениях и показана синей точкой. |
Здесь показана делеция в пояснениях и в виде пропуска нуклеотидов в покрытии. |
Программа annovar позволяет аннотировать snp, однако сначала необходимо перевести файл в формат, с которым способна работать данная программа. Сперва были удалены индели, следовательно в файле осталось только 72 snp, затем с помощью приведенной ниже команды файл snp.vcf был переведён в файл snp.avinput.
perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 snp.vcf -outfile snp.avinput |
Затем проводилось аннотирование полиморфизмов по следующим базам данных:
refgene - gene-based annotation
dbsnp - filter-based annotation
1000 genomes - filter-based annotation
Gwas - region-based annotation
Clinvar - filter-based annotation
Команда:
perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs_find.snp -build hg19 -dbtype snp138 snp.avinput /nfs/srv/databases/annovar/humandb/ |
По результатам работы программы annovar были получены три файла:
rs_find.snp.log - описание работы программы
rs_find.snp.hg19_snp138_dropped - список аннотированных в базе данных snp (имеющих rs)
rs_find.snp.hg19_snp138_filtered - список не аннотированных в базе данных snp (не имеющих rs)
снено, что 69 полиморфизмов имеют rs, а 3 не имеют. По этой ссылке Вы можете ознакомиться с примером аннотации полиморфизма в базе данных. Здесь описывается, из каких источников была получена информация о наличии полиморфизмов, характерность встречаемости данных полиморфизмов в популяциях и другая полезная информация. Неаннотированные полиморфизмы характеризуются низким качеством чтения (< 10).
Команда:
perl /nfs/srv/databases/annovar/annotate_variation.pl -out rs.refgene -build hg19 snp.avinput /nfs/srv/databases/annovar/humandb/ |
По результатам работы программы annovar были получены три файла:
rs.refgene.log - описание работы программы
rs.refgene.exonic_variant_function - описание полиморфизмов, попавших в экзоны (часть информации из файла rs.refgene.variant_function)
rs.refgene.variant_function - информация о расположении полиморфизмов относительно генов (в нетранслируемых областях, в интронах, в экзонах и др)
В Таблице 3 приведена характеристика полиморфизмов по базе данных refgene.
Таблица 3. Расположение snp. | ||
Зона |
Описание |
Число snp |
UTR3 |
snp в 3'нетранслируемой области |
5 |
intronic |
snp в интронах |
57 |
UTR5 |
snp в 5'нетранслируемой области |
2 |
exonic |
snp в экзонах |
6 |
ncRNA_exonic |
snp в транскрибируемой РНК, не имеющей аннотированного кодирующего участка |
1 |
Из данной информации следует, что наибольшее количество полиморфизмов приходится на интроны. Примерно одинаковое число полиморфизмов приходится на экзоны и 3'нетранслируемую область. Помимо этого, присутствует информация, является ли данный полиморфизм гетерозиготным (het - 30) или гомозиготным (hom - 42). Информация о полиморфизмах, попавших в экзоны представлена в Таблице 4.
Таблица 4. Характеристика snp, попавших в экзоны. | ||||||||
Ген |
Координата |
Референс |
Чтение |
Синонимичность замены |
Тип замены |
Глубина |
Качество |
Аминокислотная замена |
ATG16L1 |
234183368 |
A |
G |
Несинонимичная |
Гетерозиготная |
41 |
225.009 |
T184A |
MLPH |
238427194 |
T |
C |
Несинонимичная |
Гетерозиготная |
22 |
192.009 |
L153P |
MLPH |
238427251 |
G |
A |
Несинонимичная |
Гетерозиготная |
24 |
131.008 |
G172D |
MLPH |
238443226 |
A |
G |
Несинонимичная |
Гомозиготная |
86 |
221.999 |
H347R |
MLPH |
238449007 |
T |
C |
Несинонимичная |
Гомозиготная |
17 |
224.127 |
V374A |
MLPH |
238449107 |
A |
G |
Синонимичная |
Гомозиготная |
13 |
219.003 |
E407E |
Обозначенные в Таблице 4 гены соответствуют: MLPH (NM_024101) - ген, кодирующий белок, участвующий в связывании пигментных органелл, называемых меланосомами, с актиновым цитоскелетом в мелоноцитах, необходимый для видимой пигментации в волосах и коже [1]; ATG16L1 (NM_017974) - ген, кодирующий белок, являющийся частью белкового комплекса, необходимого для аутофагии - основного процесса, благодаря которому внутриклеточные компоненты деградируют с помощью лизосом [2].
perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs.1000g -buildver hg19 -dbtype 1000g2014oct_all snp.avinput /nfs/srv/databases/annovar/humandb/ |
По результатам работы программы annovar были получены три файла:
rs.1000g.log - описание работы программы
rs.1000g.hg19_ALL.sites.2014_10_dropped - информация о частоте встречаемости конкретных snp
rs.1000g.hg19_ALL.sites.2014_10_filtered - информация о snp, для которых не известна частота встречаемости
В результате работы программы для трех полиморфизмов не было найдено информации о частоте встречаемости, для остальных 69 она указана, причем наибольшая частота встречаемости равна 0,994209, а минимальная - 0,00738818. Информация по частоте встречаемости полиморфизмов, попавших в экзоны представлена в Таблице 5 (так как нас интересуют именно попавшие в экзоны полиморфизмы).
Таблица 5. Характеристика snp, попавших в экзоны. | ||||||||
Ген |
Координата |
Референс |
Чтение |
Тип замены |
Аминокислотная замена |
Частота встречаемости |
ATG16L1 |
234183368 |
A |
G |
Гетерозиготная |
T184A |
0.395966 |
MLPH |
238427194 |
T |
C |
Гетерозиготная |
L153P |
0.178714 |
MLPH |
238427251 |
G |
A |
Гетерозиготная |
G172D |
0.178714 |
MLPH |
238443226 |
A |
G |
Гомозиготная |
H347R |
0.371406 |
MLPH |
238449007 |
T |
C |
Гомозиготная |
V374A |
0.505591 |
MLPH |
238449107 |
A |
G |
Гомозиготная |
E407E |
0.505591 |
Как видно из Таблицы 5, частота встречаемости данных полиморфизмов довольно высока. Также с помощью программы Excel были подсчитаны среднее значение и медиана частот встречаемости полиморфизмов, они оказались равными 0,372 и 0,375 соответственно. Можно заметить, что в данном случае медиана и среднее значение оказались почти одинаковыми, из чего следует, что распределение частот встречаемости полиморфизмоф близко к нормальному.
Команда:
perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -out rs.gwas -build hg19 -dbtype gwasCatalog snp.avinput /nfs/srv/databases/annovar/humandb/ |
По результатам работы программы annovar были получены два файла:
rs.gwas.log - описание работы программы
rs.gwas.hg19_gwasCatalog - список snp, имеющих клиническое значение
В полученном файле содержалась информация по 4 из 72 полиморфизмам, значит, что лишь 4 могут иметь своё фенотипическое проявление. В Таблице 6 собрана информация по этим snp.
Таблица 6. Характеристика snp. | |||||||
Фенотипическое проявление | Координата | Референс | Чтение | Тип замены | Глубина | Качество | Расположение snp |
Рост | 55516323 | T | G | Гетерозиготная | 112 | 225.009 | UTR3 |
Болезнь Крона | 234173503 | G | A | Гетерозиготная | 23 | 113.008 | Интрон |
Болезнь Крона | 234183368 | A | G | Гетерозиготная | 41 | 225.009 | Экзон |
Рак простаты | 238443226 | A | G | Гомозиготная | 86 | 221.999 | Экзон |
Замечание: болезнь Крона (гранулематозный энтерит) - тяжёлое хроническое иммуноопосредованное гранулематозное воспалительное заболевание желудочно-кишечного тракта, которое может поражать все его отделы, начиная от полости рта и заканчивая прямой кишкой, с преимущественным поражением терминального отрезка подвздошной кишки и илеоколитом в 50 % случаев [3].
С помощью программы IGV приведенные в Таблице 6 полиморфизмы были визуализованы, так как, возможно, являются клинически значимыми и интересно узнать их расположение. На Рис. 8-11 показаны результаты такой работы. Довольно интересно, что в программе можно узнать, в какую зону входит полиморфизм, так, для первого полиморфизма видим, что под аминокислотной последовательностью нет аминокислотной последовательности, что согласуется с информацией, которая раньше была получена в базе данных Gwas: из неё следует, что данный snp входит в 3'нетранслируемую область, что отражено толстой синий полосой под нуклеотидной последовательностью, однако, если бы полиморфизм входил в состав интрона, полоса была бы тонкой, как это представлено на Рис. 9. Если полиморфизм входит в экзон, то под нуклеотидной последовательностью идёт аминокислотная, как на Рис. 10-11. Помимо этого, в отличие от первых трех полиморфизмов, у которых встречаемость разных нуклеотидов в чтениях составляет примерно 50 на 50, в последнем полиморфизме замена аденина на гуанин наблюдается в 81 из 86 чтениях, что является довольно значим результатом и согласуется с тем, что в базе данных Clinvar было найдено заболевание именно по этому полиморфизму, которое ассоциировано с данной мутацией.
Рис. 8. Фенотипическое проявление - рост. |
Рис. 9. Фенотипическое проявление - болезнь Крона (234173503). |
Рис. 10. Фенотипическое проявление - болезнь Крона (234183368). |
Рис. 11. Фенотипическое проявление - рак простаты. |
perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs.clinvar -dbtype clinvar_20150629 -buildver hg19 snp.avinput /nfs/srv/databases/annovar/humandb/ |
По результатам работы программы annovar были получены три файла:
rs.clinvar.log - описание работы программы
rs.clinvar.hg19_clinvar_20150629_dropped - список snp, аннотированных в Clinvar
s.clinvar.hg19_clinvar_20150629_filtered - список snp, не аннотированных в Clinvar
Лишь один полиморфизм является аннотированным в базе данных Clinvar: в позиции 234183368 замена T184A асооциирована с воспалительным заболеванием кишечника, этот результат согласуется с возможным фенотипическим проявлением болезни Крона, которое было аннотировано в базе данных Gwas. Остальные 71 замены не аннотированы в Clinvar. Таким образом, база данных Gwas дала более полный список возможных проблем, вызываемых присутствующими полиморфизмами. Была составлена сводная таблица, в которую вошли все snp и их характеристики по использованным для аннотации базам данных. Цветами выделены наиболее интересные и вероятные заболевания и проблемы, обусловленные данными полиморфизмами. Таблица доступна по этой ссылке snp.xlsx.
С помощью программы IGV было так же проведен анализ распределения чтений относительно экзонов, на Рис. 12 можно познакомиться с результатом этого анализа.
Рис. 12. Распределение чтений. |
Как видно на Рис. 12 чтения распределяются внутри экзонов и немного вокруг них, на далеком расстоянии от экзонов их нет (экзоны обозначены толстыми синими полосами, при увеличении изображения показывается аминокислотная последовательность).
Таблица 7. Команды использованные для выполнения данного практикума. | ||||||||
Команда |
Пояснение |
|||||||
fastqc chr2.fastq |
Проводит контроль качества чтений |
|||||||
java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr2.fastq chr2_clean.fastq TRAILING:20 MINLEN:50 |
Очищает чтения |
|||||||
fastqc chr2_clean.fastq |
Проводит контроль качества очищенных чтений |
|||||||
export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5 |
Доступ к файлам для запуска программы Hisat2 |
|||||||
hisat2-build chr2.fasta chr2 |
Индексирование референсной последовательности, выдает файлы в формате .ht2 |
|||||||
hisat2 -x chr2 -U chr2_clean.fastq --no-spliced-alignment --no-softclip > align_clean.sam |
строит выравнивание очищенных чтений и референса в формате .sam. Параметр --no-softclip запрещает подрезать чтения, а --no-spliced-alignment требует проводить картирование без разрывов. |
|||||||
samtools view align_clean.sam -bo align_clean.bam |
Переводит выравнивание в бинарный формат |
|||||||
samtools sort align_clean.bam -T 1.txt -o chr2_sort.bam |
Сортирует выравнивание чтений с референсом в формате .bam по координате в референсе начала чтения, в файл 1.txt записана выдача stdout после работы программы |
|||||||
samtools index chr2_sort.bam |
Индексирует сортированный файл |
|||||||
samtools idxstats chr2_sort.bam > out.txt |
Выводит информацию об откартированных чтениях на геном |
|||||||
samtools depth chr2_sort.bam > depth.txt |
Вычисляет покрытие каждого из нуклеотидов, и выводит эту информацию в файл depth.txt |
|||||||
samtools depth -r chr2:238395053-238463961 chr2_sort.bam > depth_ex.txt |
Вычисляет покрытие для каждого нуклеотида из указанной области |
|||||||
perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 snp.vcf -outfile snp.avinput |
Переводит файл snp.vcf в формат, с которым может работать программа annovar |
|||||||
perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs_find.snp -build hg19 -dbtype snp138 snp.avinput /nfs/srv/databases/annovar/humandb/ |
Позволяет узнать, какие snp имеют rs |
|||||||
perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -out rs.gwas -build hg19 -dbtype gwasCatalog snp.avinput /nfs/srv/databases/annovar/humandb/ |
Поиск по Gwas |
|||||||
erl /nfs/srv/databases/annovar/annotate_variation.pl -out rs.refgene -build hg19 snp.avinput /nfs/srv/databases/annovar/humandb/ |
Поиск по RefGene |
|||||||
perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs.clinvar -dbtype clinvar_20150629 -buildver hg19 snp.avinput /nfs/srv/databases/annovar/humandb/ |
Поиск по Clinvar |
|||||||
perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs.1000g -buildver hg19 -dbtype 1000g2014oct_all snp.avinput /nfs/srv/databases/annovar/humandb/ |
Поиск по 1000 genome |
[1]Melanophilin [2]ATG16L1 [3]Болезнь Крона |