В предыдущих статьях (1,2) мы познакомились с тем, как могут выглядеть данные в зависимости от проведенного биологического эксперимента. На основании этих визуализированных данных были сделаны предположения о том, что же происходит внутри клетки. Теперь остановимся на том, как математически и алгоритмически проанализировать данные для того, чтобы машины за нас могли выполнить рутинную работу. К сожалению, после прочтения множества статей по анализу данных у меня сложилось впечатление, что однозначного или наиболее универсального решения не существует. Есть алгоритмы, которые хорошо себя показывают на некотором наборе данных, а в других случаях уже не отвечают поставленным задачам.
Для начала хочу напомнить и немного пояснить основные используемые в статье понятия. Каждый раз, когда рассказываю о ридах, и каким образом они получаются, наталкиваюсь на удивление и недопонимание. И тут случайно, закидывая пачку бумаги в шредер, понял, как можно попробовать объяснить, что такое риды. Допустим, у вас есть книга, сама книга — это клетка, а каждая глава книги — это ДНК. У вас в наличии аж несколько миллионов экземпляров этой книги. Вы скормили книги шредеру, который разрезал их на полоски. Вот вы берете контейнер с накромсанными бумажными полосками, хорошо перемешиваете их и достаете то, что первое попалось. Далее аккуратно раскладываете полоски, с двух сторон отрезаете по маленькому кусочку и складываете отрезанные части в отдельную коробочку. Эти маленькие кусочки — и есть риды, полоски — это фрагменты. Далее с помощью некой программы мы определяем, в какой части книги находился отрезанный кусочек. Иногда нам везёт, и мы можем точно определить место в книге, откуда был вырезан кусочек. Иногда программа определяет несколько возможных вариантов расположения этого кусочка, например, из-за того, что некоторые главы начинались с повторения предыдущей или в результате всех операций часть букв стерлась, и они не нашлись в книге. Количество вариаций на тему, что могло произойти с кусочками полосок, очень велико, если все их обсуждать, можно и до алгоритмов не дойти. Надеюсь, теперь стало понятнее, что есть что.
Описывая алгоритм, я буду рассказывать о нем, используя те данные, на которых его тестировали разработчики. Но это не означает, что его нельзя использовать на данных другого вида. Также постараюсь выявить недостатки алгоритмов, с которыми я столкнулся. Систематизировать алгоритмы я не берусь, буду описывать их группами, в том порядке, в котором они применялись при общем анализе секвенированных данных. В этой статье я расскажу о простом анализе “загрязнение адаптером” и о предварительной подготовке ландшафта и уже, наверное, в следующей статье или через одну — об оценке качества полученных данных.
Загрязнение адаптером
Компании, производящие оборудование для секвенирования, также производят и сопутствующие наборы: слайды, пробирки, реагенты. Один из видов реагентов — это полинуклеотидная цепочка, необходимая для секвенирования, которая называется адаптером.
Встречаются различные виды загрязнения подготавливаемого материала. Один из них — загрязнение адаптером. Анализ такого загрязнения прост. Необходимо подсчитать частоту нуклеотидов на первой, второй и т.д. позициях полученного FASTA/FASTQ файла. На основе полученных данных составляется таблица. Рассмотрим процесс анализа более подробно, используя следующий пример, ниже представлены первые 9 столбцов таблицы (см. Таблица 1):
Таблица 1
Pos/Nucl
0
1
2
3
4
5
6
7
A
446500
2869138
599400
941816
1145756
1133404
3599581
3291736
C
1922795
1367021
1761361
3723679
1030494
1173573
963827
1284534
T
485441
1337742
3273475
918001
1016981
1006168
1494250
1291715
G
4130093
1425977
1365764
1416504
3806769
3686855
942342
1132015
N
15171
122
Используя данные таблицы, можно построить график (Рис. 1). Пики графика соответствуют наиболее часто встречающейся на текущей позиции букве (нуклеотиду). В частности, буква «G» на первой, буква «A» на второй, «T» на третьей позициях.
Рис. 1
Если составить всю цепочку из соответствующих пикам букв, то получится последовательность, которая наиболее часто встречается в исходных данных. Если эта последовательность совпала с последовательностью из списка адаптеров, (которые поставляют производители оборудования для секвенирования), как например, в рассматриваемом примере получилась одна из последовательностей адаптера компании Illumina, то значит произошло загрязнение адаптером.
Чтобы вычислить процентное соотношение загрязнения адаптером и исходного материала, выполняем следующие простые действия. 1. Подсчитываем общее количество букв в одной из позиций в таблице, в приведенном примере их ровно 7000000. 2. Если загрязнения не было, то количество букв на каждой позиции должно быть примерно одинаковым. Исключив из рассмотрения букву, соответствующую пику, считаем оставшиеся. Как следует из графика, таких букв в каждой позиции, примерно, 1200000. 3. Так как на каждой позиции должно быть одинаковое количество букв, то получаем: 1200000 * 4 = 4800000. 4. Таким образом, загрязнение адаптером составляет: 4800000 / 7000000 * 100% = 32%. Т.е. при подготовке материала в следующий раз надо убавить количество добавляемых адаптеров на 32%.
Подготовка профайла для последующего анализа
Большинство известных мне алгоритмов выполняют первоначальную обработку данных — подготавливают профайл. Т.е группируют данные, а в дальнейшем оценивают их уровень значимости. Сгруппированные данные будут определять место вероятного расположения белка на ДНК, которое мы пытаемся найти с помощью алгоритмов.
Рассмотрим данные, получающиеся в результате использования иммунопреципитации хроматина (подробнее “Секвенирование ДНК”). Напомню об особенностях этого метода: фрагменты выбираются с помощью антител, соответствующих интересующим белкам, в результате, большинство фрагментов содержит последовательность, к которой прикреплялся белок. Также при нарезке получаются фрагменты разной величины, дополнительно их фильтруют по длине (не более, чем 150-300 нуклеотидов). Большинство ридов в конечном результате расположены с двух сторон относительно белка. Этой особенностью мы и будем пользоваться в дальнейшем для определения вероятного местонахождения белка.
Рис. 2
На рисунке 2 желтыми стрелками показаны возможные места разреза ДНК, именно здесь могут начинаться или заканчиваться фрагменты. Красными линиями изображены риды, полученные с положительной спирали (на рисунке обозначена знаком “+”), а синими линиями — с отрицательной (обозначена знаком “-”). В некоторых случаях, если проводится более дорогой эксперимент (“pair-end sequencing”), аппаратура секвенирования способна определить, какие два рида были получены с одного фрагмента. Также утверждается, что использование такого метода в несколько раз увеличивает качество определения предполагаемой позиции белка [1]. Метод “pair-end” пока еще используется не часто, поэтому в большинстве случаев мы имеем дело с несвязанным набором ридов.
В идеальном случае ландшафт расположения ридов вокруг белка должен быть, как на рисунке 3. Учитывая, что возможное количество случайных составляющих эксперимента велико, такое встречается редко.
Рис. 3
Ниже приведены примеры ландшафтов из экспериментальных данных (Рис. 4). Ширина столбца (окна) соответствует 20 нуклеотидам, максимальная высота столбцов составляет всего 5-8 ридов. Если бы ширина окна была равна одному нуклеотиду, то изображение стало бы ещё более дискретным и шумным. Но, если данные обработать предлагаемыми ниже алгоритмами по подготовке профайла, то ландшафт будет более похожим на идеальный (Рис. 3).
Рис. 4
Ручная подготовка профайла
Как же привести данные к более читаемому для машины виду? В некоторых случаях мы заранее знаем, какой длины ожидаются фрагменты. Например, в случае H3K4Me3 (MNase метод) они составляют около 148bp [2]. Сдвигаем риды на половину длины фрагмента: 148 / 2 = 74. Риды, найденные на положительной спирали, сдвигаем вправо на 74, найденные на отрицательной спирали, сдвигаем влево на 74. Складываем абсолютные величины и отображаем на графике (Рис. 5) [3, 4].
Рис. 5
На рисунке (Рис. 5) фиолетовым цветом отображен пик, получившийся в результате сложения количества ридов с положительной и отрицательной сторон спирали. Вершина пика должна с некоторой точностью соответствовать центру расположения белка на ДНК. Таким простым преобразованием мы усилили сигнал, и теперь он должен стать более заметным.
Существуют различные способы построения графиков, приведённых выше. Один из способов: по оси “Y” откладываем количество ридов, берущих начало в точке “Х”. Другой способ: построение графика с помощью сегментов, по оси “Y” откладываем количество пересечений ридов в заданной точке “Х”. Иногда график строят, откладывая по оси “Y” пересечение не ридов, а фрагментов. Смысл таких графиков — плотность покрытия ДНК ридами/фрагментами.
Это хорошо, когда можно заранее оценить размер фрагмента и передать его значение в программу. Но, что делать, если длина фрагмента не известна, как мы можем вычислить величину, на которую необходимо сдвигать данные?
Кросс-корреляционная подготовка профайла
Один из способов вычисления величины сдвига — это кросс-корреляция данных о ридах с положительной и отрицательной спирали. Для каждой хромосомы “c” определим функцию , которая принимает в каждой точке “x” значение, равное количеству начал ридов на стороне спирали “s”. Определим для сдвига функцию , где:
— линейный коэффициент Пирсона между векторами “a” и “b”;
C — набор всех хромосом;
Nc — количество ридов на хромосоме “c”;
N — общее количество ридов.
Выберем диапазон, в котором может находиться величина сдвига, например, от -100 до 300 и построим график (Рис. 6) [4].
Рис. 6
Пик этого графика по оси “x” будет соответствовать сдвигу с наибольшей корреляцией. В случае, представленном на рисунке 6, сдвиг равен 100 / 2 = 50 нуклеотидам. Коррелирование данных с учетом всех хромосом может оказаться менее точным относительно коррелирования интересующего участка. Таким образом получают свою сдвига для каждого обогащенного сегмента.
Kernel density профайл
Рассмотрим алгоритм с использованием kernel density функции [4,5,6]. Переопределим плотность покрытия для каждой позиции i в соответствии с формулой , где:h — kernel density bandwidth, в статье [5] использовалось 30bp, а в статье [4] взяли 0.45 от ширины пика в середине, изображенного на рисунке 6;
— Гауссова kernel density функция;
— возвращает количество ридов, окончившихся в точке j на положительной и отрицательной спирали, соответственно.
Таким образом, пересчитав плотность для каждой спирали и усилив сигнал, для хорошо обогащенных участков мы можем определить пики на каждой стороне спирали и взять половину разности между ними за величину сдвига (Рис. 7, “Peak shift”).
Рис. 7
Геометрический метод построения профайла
Этот алгоритм использует нехитрые математические операции. Для окна в ширину w и для позиции i определим: p1(i), как количество тегов в окне на положительной спирали слева от i; p2(i), как количество тегов справа от позиции i. Определим n1(i) и n2(i) соответственно определению p, только для отрицательной спирали (Рис 8.1). Подсчитаем значение в точке i по следующей формуле: в пределах заданного окна находим максимум. Найденный максимум — это предполагаемая позиция закрепления белка. Этот метод назвали WTD от Windows Tag Density.
В статье [4] его также предложили модифицировать , где P[a,b] — это линейный коэффициент корреляции Пирсона между векторами a и b, такими что, вектор “a” — это количество ридов на положительной спирали от позиции i до i+k, а вектор “b” — это количество ридов на отрицательной спирали от позиции i до i-k, (Рис 8.2). Назвали MTC от Mirrow Tag Correlation.
Рис. 8
Есть и другие алгоритмы, используемые в построении профайла [7]. Я постарался осветить наиболее часто применяемые. У всех алгоритмов уже есть реализация в каком-либо программном пакете, ссылки приведены в статьях, вот некоторые из них:sissrs.rajajothi.com/http://home.gwu.edu/~wpeng/Software.htmmendel.stanford.edu/sidowlab/downloads/quest/compbio.med.harvard.edu/software.html.
Правильное построение профайла с помощью алгоритмов — наиболее важная часть всех программ, так как, чем точнее он построен, тем с наибольшей долей вероятности можно определить, где находится белок. Другая часть алгоритмов направлена на то, чтобы математически оценить, насколько значимы найденные участки. Но, к сожалению, все алгоритмы не способны заменить человека, а могут только помочь в рутинной обработке данных, намекая на место, где следует вести поиски. Так что пока вся система взаимодействия биологии, химии, математики и программирования далека от совершенства, и существует огромное поле для поиска, реализации и совершенствования идей.
1. Wang, C., et al., An effective approach for identification of in vivo protein-DNA binding sites from paired-end ChIP-Seq data. BMC Bioinformatics, 2010. 11: p. 81. www.ncbi.nlm.nih.gov/pmc/articles/PMC2831849/
2. Barski, A., et al., High-resolution profiling of histone methylations in the human genome. Cell, 2007. 129(4): p. 823-37. www.ncbi.nlm.nih.gov/pubmed/17512414
3. Pepke, S., B. Wold, and A. Mortazavi, Computation for ChIP-seq and RNA-seq studies. Nat Methods, 2009. 6(11 Suppl): p. S22-32. www.ncbi.nlm.nih.gov/pubmed/19844228
4. Kharchenko, P.V., M.Y. Tolstorukov, and P.J. Park, Design and analysis of ChIP-seq experiments for DNA-binding proteins. Nat Biotechnol, 2008. 26(12): p. 1351-9. www.ncbi.nlm.nih.gov/pubmed/19029915
5. Valouev, A., et al., Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data. Nat Methods, 2008. 5(9): p. 829-34. www.ncbi.nlm.nih.gov/pubmed/19160518
6. Boyle, A.P., et al., F-Seq: a feature density estimator for high-throughput sequence tags. Bioinformatics, 2008. 24(21): p. 2537-8. www.ncbi.nlm.nih.gov/pmc/articles/PMC2732284/
7. Zang, C., et al., A clustering approach for identification of enriched domains from histone modification ChIP-Seq data. Bioinformatics, 2009. 25(15): p. 1952-8.
Часть рисунков была взята из журнала nature.com и напечатана с их разрешения:
License Date: Feb 16, 2012
License Numbers: 2850830942235, 2850831421638
Review is prepared by Andrey Kartashov, Covington, KY, porter@porter.st
