Определение сайтов связывания данного транскрипционного фактора в данном участке хромосомы человека
Файлы .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. Качество определения основания в каждой позиции рида.
Ось Y на графике отображает "quality scores" - оценки качества. Чем лучше оценка, тем более точно определено основание. Задний фон графика разделяет ось Y на очень хорошее качество определения (зеленый), разумное (среднее) качество (оранжевый) и плохое качество (красный). Качество определения на большинстве платформ обычно падает с течением процесса, поэтому даже падение качества до оранжевой области в конце является нормой. [1] Однако, как видно из Рис. 1, качество ридов оказалось очень высоким даже в конце, поэтому не было нужды в фильтрации данных программой Trimmomatic. На Рис. 2 представлена основная информация о чтениях. Как видно на изображении, всего было 6572 чтения длиной 36 нуклеотидов каждое.
Рис 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. Информация о пиках, найденных программой 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. Визуализация пиков в UCSC Genome Browser.
Как видно из Рис. 4, только 9-й пик перекрывается с функциональным элементом генома, а точнее попадает в область гена OSBPL8 (Oxysterol binding protein-like 8). Остальные же пики с какими-либо генами не перекрываются и находятся далеко от них.
Два наиболее достоверных пика (4 и 6), а также пик 9 были рассмотрены в более крупном масштабе. Они представлены на Рис. 5-7. На этих изображениях также показано расположение пиков относительно других факторов транскрипции из ранее проведенных ChIP-seq экспериментов. Из изображений ясно, что достоверность пиков подтверждается результатами других экспериментов.
Рис 5. Увеличенный пик 4.
Рис 6. Увеличенный пик 6.
Рис 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. Выравнивание последовательностей, построенное сервером MUSCLE и визуализированное в программе JalView.
Как видно из Рис. 8, выравнивание получилось очень плохим. В нем очень много гэпов и слишком мало консервативных позиций. На первый взгляд в выравнивании нельзя выявить приемлимые консенсусные последовательности длиной 4-8 нуклеотидов. Для подтверждения (или опровержения) этого факта последовательности были даны на вход программе МЕМЕ. И действительно, хоть программа и смогла найти 3 мотива, но все находки обладали огромным e-value (все были > 1), и достоверными считаться никак не могли. Вероятно, поскольку экспериментальные данные были все же неполноценны, не удалось определить последовательности мотивов достаточно точно.