Определение сайтов связывания данного транскрипционного фактора в данном участке хромосомы человека
Файлы .fastq с ридами Illumina, полученные в результате сhip-seq эксперимента, были разделены на отдельные файлы, соответствующие участкам хромосом. Для выполнения данного задания был дан файл chipseq_chunk5.fastq
Контроль качества чтений
Для начала при помощи програмы FastQC (команда: fastqc chipseq_chunk5.fastq) был сделан контроль качества чтений. Результат работы программы представлен в файлах chipseq_chunk5_fastqc.html и chipseq_chunk5_fastqc.zip.
На Рис. 1. представлено качество определения основания в каждой позиции рида в FastQ-файле. Для каждой позиции нарисована диаграмма размаха, т.н. "ящик с усами" - график, компактно изображающий одномерное распределение вероятностей (показывает медиану (или, если нужно, среднее), нижний и верхний квартили, минимальное и максимальное значение выборки и выбросы). На графике присутствуют следующие элементы [1]:
- горизонтальные красные линии отображают средние значения;
- желтые "ящики" представляют межквартильный диапазон (25-75%);
- верхние и нижние "усы" представляют собой 10-й и 90-й процентили (качество 10% чтений не ниже нижней границы, 90% - не выше верхней);
- синяя линия показывает среднее качество.
![](1.png)
Рис 1. Качество определения основания в каждой позиции рида.
Ось Y на графике отображает "quality scores" - оценки качества. Чем лучше оценка, тем более точно определено основание. Задний фон графика разделяет ось Y на очень хорошее качество определения (зеленый), разумное (среднее) качество (оранжевый) и плохое качество (красный). Качество определения на большинстве платформ обычно падает с течением процесса, поэтому даже падение качества до оранжевой области в конце является нормой. [1] Однако, как видно из Рис. 1, качество ридов оказалось очень высоким даже в конце, поэтому не было нужды в фильтрации данных программой Trimmomatic. На Рис. 2 представлена основная информация о чтениях. Как видно на изображении, всего было 6572 чтения длиной 36 нуклеотидов каждое.
![](2.png)
Рис 2. Основная информация о чтениях.
Картирование на геном
Далее было проведено картирование прочтений на геном человека hg19.
Команда: /srv/databases/ngs/hg19/GRCh37.p13.genome.fa chipseq_chunk5.fastq > chipseq_chunk5.sam
Поскольку для работы был доступен уже проиндексированный геном, индексирование проводить было не нужно. В результате было построено выравнивание чтений и референсной последовательности (hg19) в формате .sam.
Анализ выравнивания
Анализ был проведен при помощи пакета samtools. Выполненные при этом программы и их описания представлены ниже:
- samtools view -bSo chipseq_chunk5.bam chipseq_chunk5.sam - перевод выравнивание чтений с референсом в бинарный формат .bam
- samtools sort chipseq_chunk5.bam -T chip_temp -o chipseq_chunk5.sorted.bam - сортировка выравнивания чтений с референсом по координате в референсе начала чтения
- samtools index chipseq_chunk5.sorted.bam - индексирование отсортированного .bam файла
- samtools idxstats chipseq_chunk5.sorted.bam > chipseq_chunk5.idxstats - информация о количестве чтений, откартированных на геном
- samtools view -c chipseq_chunk5.sorted.bam - вывод на экран общего количества откартированных чтений
Изначально было 6572 рида, и столько же было откартировано на геном. На основе данных файла chipseq_chunk5.idxstats можно предположить, что для анализа были предложены чтения с 12-й хромосомы, т.к. именно на нее откартировалось больше всего чтений - 6185.
Поиск пиков
Далее при помощи программы MACS был проведен поиск пиков (команда: macs2 callpeak -t chipseq_chunk5.sorted.bam --nomodel -n MyChunk). В результате были получены 3 файла: MyChunk_peaks.narrowPeak, MyChunk_peaks.xls и MyChunk_summits.bed, содержащие информацию о найденных пиках.
Всего было найдено 9 пиков шириной от 213 до 497 нуклеотидов, расположенных в одном регионе 12-й хромосомы. Более подробная информация них представлена на Рис.3.
![](3.png)
Рис. 3. Информация о пиках, найденных программой MACS
Ширины пиков представлена в колонке "length", а позиции вершин пиков относительно начала - в колонке "position from start" Таблицы 1. Достоверность пиков можно определить, например, по p-value. Соответственно, т.к. пик 6 обладает наибольшим –log10p-value, то является самым достоверным, а наименее достоверным является пик 8.
Визуализация в UCSC Genome Browser
Информация из файла MyChunk_peaks.narrowPeak была визуализирована при помощи геномного браузера UCSC Genome Browser. Для этого в указанный файл были дописаны дополнительные строки:
track type=narrowPeak visibility=3 db=hg19 name="my_peaks" description="Peaks from chunk 5" browser position chr12:76116000-76865000
![](4.png)
Рис 4. Визуализация пиков в UCSC Genome Browser.
Как видно из Рис. 4, только 9-й пик перекрывается с функциональным элементом генома, а точнее попадает в область гена OSBPL8 (Oxysterol binding protein-like 8). Остальные же пики с какими-либо генами не перекрываются и находятся далеко от них.
Два наиболее достоверных пика (4 и 6), а также пик 9 были рассмотрены в более крупном масштабе. Они представлены на Рис. 5-7. На этих изображениях также показано расположение пиков относительно других факторов транскрипции из ранее проведенных ChIP-seq экспериментов. Из изображений ясно, что достоверность пиков подтверждается результатами других экспериментов.
![](5.png)
Рис 5. Увеличенный пик 4.
![](6.png)
Рис 6. Увеличенный пик 6.
![](7.png)
Рис 7. Увеличенный пик 9.
Выравнивание последовательностей
Файл MyChunk_peaks.narrowPeak был конвертирован в fasta-формат.
Команда: /srv/databases/ngs/tools/seqtk/seqtk subseq /srv/databases/ngs/hg19/GRCh37.p13.genome.fa MyChunk_peaks.narrowPeak > my_peaks.fa
Результатом служил сам файл my_peaks.fa, содержащий вырезанные участки генома человека с пиками. При помощи сервера MUSCLE было построено выравнивание и затем открыто в программе JalView (изображение выравнивания приведено на Рис. 8).
![](8.png)
Рис. 8. Выравнивание последовательностей, построенное сервером MUSCLE и визуализированное в программе JalView.
Как видно из Рис. 8, выравнивание получилось очень плохим. В нем очень много гэпов и слишком мало консервативных позиций. На первый взгляд в выравнивании нельзя выявить приемлимые консенсусные последовательности длиной 4-8 нуклеотидов. Для подтверждения (или опровержения) этого факта последовательности были даны на вход программе МЕМЕ. И действительно, хоть программа и смогла найти 3 мотива, но все находки обладали огромным e-value (все были > 1), и достоверными считаться никак не могли. Вероятно, поскольку экспериментальные данные были все же неполноценны, не удалось определить последовательности мотивов достаточно точно.