I

Для определения сайтов связывания транскрипционного фактора в заданном участке хромосомы человека использовался файл chipseq_chunk50.fastq, в котором содержатся риды, полученные при помощи Illumina в результате ChIP-Seq эксперимента. Изначально осуществлялся контроль качества чтений при помощи FastQC, результат до очистки представлен ниже.

Если верить краткой сводке программы FastQC, то файл не содержит последовательностей низкого качества. Однако видно, что наблюдаемые минимумы («усы») последовательностей часто заходят в красную область. Поэтому было решено провести очистку при помощи программы Trimmomatic. Результат работы программы FastQC после очистки представлен ниже.

Наблюдаемая теперь картина отличается от того, что было раньше. Теперь качество по каждому нуклеотиду полностью находится в благоприятной (> 28) зоне (выделена зеленым). Этому способствовало использование trailing-отсчечения: удаления концевых нуклеотидов, качество которых ниже порогового.

Далее улучшенные риды были откартированы на проиндексированный геном человека при помощи алгоритма BWA-MEM. Анализ картирования (sam-выравнивание: chunk50.sam) проводился при помощи samtools: view для перевода в бинарный формат, sort для сортировки выравнивания по координате начала рида в референсе, index для индексирования, и, наконец, idxstats для получения информации о количестве прочтений на каждой из хромосом.

Результаты картирования отображены в данной таблице. Видно, что подавляющее большинство чтений легли на вторую хромосому (18117), поэтому можно предположить, что изначально анализировалась именно она. Всего изначально было 22730 ридов, и только часть из них (18824) картировалась на геном (остальные 3906 — нет).

Процедура поиска пиков осуществлялась с использованием программы MACS2, запуск которой был проведен следующей командой:

Параметр outdir задает директорию (callpeak для файлов, в названии которых присутствует имя эксперимента, указываемое в -n (--name). Флаг --nomodel был прописан, поскольку обычный запуск программы включает в себя построение некоторой модели shifting model, что оказывается невозможным в условиях малого количества пиков (что у нас и наблюдалось).

В результате работы было получено три файла: [narrowPeak], [xls], [bed]. Содержание chunk43_peaks.xls приведено в таблице ниже.

Всего нашлось 15 пиков на второй хромосоме, ширина которых варьирует от 235 (NA_peak_14) до 746 (NA_peak_7). Стоит отметить, что колонки 7 и 9 содержат отрицательные десятичные логарифмы p-value и q-value, а значит трактовать их надо не так, как обычные значения. Здесь работает принцип «чем больше, тем лучше», поскольку при повышении этих значений p-value и q-value падают. По ним можно судить о достоверности пиков, и, согласно этому критерию, все находки оказались достаточно надежными. Например, самый надёждный пик (NA_peak_7) имеет показатели -LOG10(pvalue) и -LOG10(qvalue) 92.6 и 83.25 соответственно. Самые плохие показатели получились у NA_peak_8 (значения 10.1 и 4.2 соответственно). Но тем не менее очевидно, что и этот пик можно назвать достоверным с большой долей уверенности несмотря на то, что на фоне остальных он вышел хуже всего.

Следующая задача — визуализация полученных пиков. Для этого был использован геномный браузер UCSC Genome Browser. Для того чтобы использовать файл [narrowPeaks] в качестве трэк-файла, в него добавились следующие строки:

На рисунке ниже приведён скриншот из геномного браузера.

Для начала был рассмотрен самый достоверный пик — седьмой. Вот его характеристики:

• Координаты в геноме: 102183035..102183781;
• Ширина пика: 352;
• Расположение вершины пика относительно начала: 155;
• Достоверность: -LOG10(q-value) = 35.21821, -LOG10(p-value) = 42.06512.

Далее было решено рассмотреть на третий пик. У него наблюдаются следующие характеристики:

• Координаты в геноме: 101733691..101734034;
• Ширина пика: 832 (самый длинный);
• Расположение вершины пика относительно начала: 466;
• Достоверность: -LOG10(q-value) = 83.25476, -LOG10(p-value) = 92.61709 (самый достоверный).