В этом примере показано, как использовать Байесов метод скрытой модели Маркова (HMM), чтобы идентифицировать изменение номера копии в основанных на массиве данных о сравнительной геномной гибридизации (CGH).
Основанный на массиве CGH является высокопроизводительным методом, чтобы измерить изменение номера копии ДНК через геном. Фрагменты ДНК или "клоны" теста и ссылочных выборок гибридизируются к сопоставленным фрагментам массивов. Отношения интенсивности Log2 теста к ссылке предоставляют полезную информацию о профилях всего генома в номере копии. В идеальной ситуации log2 отношение нормальных (нейтральных в отношении копии) клонов является log2 (2/2) = 0, одна потери копии log2 (1/2) =-1, и одно усиления копии являются log2 (3/2) = 0.58. Несколько усилений копии или усилений имели бы значения log2 (4/2), log2 (5/2).... Потеря обеих копий или удаление соответствовала бы значению-inf. В действительных приложениях, даже после составления погрешности измерения, log2 отношения значительно отличаются от теоретических значений. Отношения обычно уменьшаются по направлению к нулю. Одним основным фактором является загрязнение выборок опухоли с нормальными ячейками. Существует также зависимость между отношениями интенсивности соседних клонов. Эти проблемы требуют использования эффективных статистических алгоритмов, характеризующих геномные профили.
Guha и др., (2006) предложил Байесов статистический подход в зависимости от скрытой модели Маркова (HMM) для анализа массива данные CGH. Скрытая модель Маркова составляет зависимость между соседними клонами. Отношения интенсивности сгенерированы скрытыми состояниями номера копии. Байесово изучение используется, чтобы идентифицировать изменения всего генома в номере копии из данных. Следующие выводы сделаны о прибылях и убытках номера копии.
В этом Байесовом алгоритме HMM существует четыре состояния, заданные как состояние номера копии потерь (1), копируют номер нейтральное состояние (2), одно состояние усиления копии (3), и усиление (несколько получают), состояние (4). Состояние номера копии сопоставлено с каждым клоном. Нормированные log2 отношения приняты, чтобы быть распределенными как
mu является неизвестным параметром для каждого состояния с этим ограничением:
Уголовное прошлое для средних изменений номера копии:
![$$[\mu_4|\mu_3, \sigma_3] \sim N(1,\tau_4^2) \cdot I(\mu_4>\mu_3 +
3\sigma_3)$$](../../examples/bioinfo/win64/acghhmmdemo_eq02275008633568004220.png)
Guha и др., (2006) также описал алгоритм Столицы в Гиббсе, чтобы сгенерировать следующие выборки. Алгоритм MCMC независимо запущен для каждой хромосомы, чтобы сгенерировать выборку MCMC для параметров хромосомы. Начальные значения параметров сгенерированы от уголовного прошлого. Сгенерированные состояния номера копии представляют, чертит от крайнего следующего из интереса, Поскольку каждый MCMC чертит, сгенерированные состояния смотрятся и классифицируются как фокальные абляции, точки перехода, усиления, выбросы и целые хромосомные изменения.
В этом примере вы примените Байесов алгоритм HMM, чтобы анализировать массив профили CGH некоторых выборок рака поджелудочной железы [2].
Данные в этом примере являются массивом профили CGH 24 клеточных линий аденокарциномы поджелудочной железы и 13 первичных экземпляров опухоли от Alguirre и др., (2004). Пометил фрагменты DNA, были гибридизированы к микромассивам комплементарной ДНК человека Agilent®, содержащим 14 160 клонов комплементарной ДНК. Приблизительно 9 420 клонов имеют уникальные положения карты со средним интервалом между сопоставленными элементами 100,1 Кбит. Больше деталей данных и эксперимента может быть найдено в [2]. Для удобства были уже нормированы данные, и log2 базировался, отношения интенсивности обеспечиваются файлом MAT pancrea_oligocgh.mat.
Вы примените Байесов алгоритм HMM, чтобы анализировать хромосому 12 из демонстрационных 6 из данных об аденокарциноме поджелудочной железы и сравнить результаты с сегментами, найденными алгоритмом круговой бинарной сегментации (CBS) Olshen и др., (2004).
Загрузите файл MAT, содержащий log2 отношения интенсивности, и сопоставил геномные положения этих 37 выборок.
load pancrea_oligocgh
pancrea_data
pancrea_data =
struct with fields:
Sample: {37x1 cell}
Chromosome: [13446x1 int8]
GenomicPosition: [13446x1 int32]
Log2Ratio: [13446x37 double]
Log2RatioMed: [13446x37 double]
Log2RatioSeg: [13446x37 double]
CloneIDs: [13446x1 int32]
Задайте номер хромосомы и выборку, чтобы анализировать.
sampleIndex = 6;
chromID = 12;
sample = pancrea_data.Sample{sampleIndex}
sample =
'PA.C.Dan.G'
Загрузите и отобразите log2 данные об отношении на графике хромосомы 12 от демонстрационного PA.C.Dan.G.
idx = pancrea_data.Chromosome == chromID; X = double(pancrea_data.GenomicPosition(idx)); Y = pancrea_data.Log2Ratio(idx, sampleIndex); % Remove NaN data points idx = ~isnan(Y); X = X(idx); Y = Y(idx); % Plot the data figure; plot(X, Y, '.', 'color', [0.6 0.6 1]) ylims = [-1.5, 3.5]; ylim(gca, ylims) title(sprintf('%s - Chromosome %d', sample, chromID)) xlabel('Genomic Position'); ylabel('Log2(Ratio)')
Количество клонов на хромосоме 12, чтобы анализироваться
N = numel(Y)
N = 437
Можно запустить анализ путем выполнения хромосомной сегментации с помощью алгоритма CBS [3], который реализован в cghcbs функция. Процесс займет несколько секунд. Можно просмотреть график средних значений сегмента по исходным данным путем определения SHOWPLOT параметр. Примечание: можно ввести doc cghcbs для получения дополнительной информации об этой функции.
PS = cghcbs(pancrea_data, 'SampleInd', sampleIndex, ... 'Chromosome', chromID, 'ShowPlot', chromID); ylim(gca, ylims)
Analyzing: PA.C.Dan.G. Current chromosome 12
Как показано на рисунке, процедура CBS объявила набор отношений высокой интенсивности как два отдельных сегмента. Процедура CBS также нашла область с потерями номера копии.
Байесов подход HMM использует алгоритм Столицы в Гиббсе, чтобы сгенерировать следующие выборки параметров [1]. Параметры модели сгруппированы в четыре блока. Алгоритм итеративно генерирует каждое четыре условных выражения блоков на оставшихся блоках и данных.
Чтобы анализировать данные с Байесовым алгоритмом HMM, необходимо инициализировать параметры. Больше деталей о предшествующих параметрах может быть найдено в ссылках [1] и [4].
Инициализируйте состояние генератора случайных чисел, чтобы гарантировать, что фигуры, сгенерированные ими, управляют, соответствуют рисункам в версии HTML этого примера.
rng('default');
Задайте количество состояний
NS = 4;
Задайте количество итераций MCMC
NMC = 100;
Определите гиперпараметры предшествующих распределений для четырех состояний.
mus_hyper = [-1, 0, 0.58, 1]; taus_hyper = [1, 1, 1, 2];
Установите эпсилон параметра, который определяет ограничивание средних значений.
eps = 0.1;
Установите границы предшествующих средних значений каждого состояния.
mu_low_bounds = [-Inf, -eps, eps, 0.58]; mu_up_bounds = [-eps, eps, 0.58, Inf];
Guha и др., (2006) принимает инверсию предшествующих ошибочных отклонений (sigma^2) как гамма распределения с нижними границами 0,41 для состояний 1, 2 и 3. Установите масштабные коэффициенты для гамма распределений для каждого состояния.
sg_alpha = [1 1 1 1]; sg_beta = [1, 1, 1, 1]; sg_bounds = [0.41 0.41 0.41 1];
Задайте переменную states чтобы сохранить номер копии утверждают последовательности клонов для каждой итерации MCMC.
states = zeros(N, NMC);
Задайте переменную st_counts содержать изменение состояния значит каждое состояние номера копии.
st_counts = zeros(NS, NS);
Итерация MCMC запускается в
iloop = 1;
Определите сигмы для четырех состояний путем выборки от гамма распределения с предшествующей альфой масштабного коэффициента и бетой.
sigmas = zeros(NS, NMC); for i = 1:NS sigmas(i, iloop) = acghhmmsample('gamma', sg_alpha(i), sg_beta(i), sg_bounds(i)); end
Определите средние значения для четырех состояний путем выборки от усеченного нормального распределения между нижними и верхними границами средних значений. Примечание: четвертая нижняя граница состояния будет определена третьим состоянием.
mus = zeros(NS, NMC); for i = 1:NS if i == 4 mu_low_bounds(4) = mus(3,iloop) + 3*sigmas(3,iloop); end mus(i, iloop) = acghhmmsample('normal', mus_hyper(i), taus_hyper(i),... mu_low_bounds(i), mu_up_bounds(i)); end
Примите независимое уголовное прошлое Дирихле для строк стохастического 4x4 матрица [1] вероятности перехода. Сгенерируйте стохастическую предшествующую матрицу А перехода от распределений Дирихле.
a = ones(NS, NS);
A = acghhmmsample('dirichlet', a, NS);
Матрица перехода имеет уникальное стационарное распределение. Стационарный PI распределения является собственным вектором матрицы перехода, сопоставленной с собственным значением 1.
PI =@(x, n) (ones(1,n)/(eye(n) -x + ones(n)))';
Сгенерируйте предшествующий стационарный PI распределения.
Pi = PI(A, NS);
Сгенерируйте начальную матрицу B эмиссии
B = zeros(NS, N); for i = 1:NS B(i,:) = normpdf(Y, mus(i,iloop), sigmas(i,iloop)); end
Декодируйте начальную букву скрытые состояния клонов с помощью стохастического прямого обратного алгоритма [4].
states(:, iloop) = acghhmmfb(Pi, A, B);
Для каждой итерации MCMC четыре блока параметров сгенерированы можно следующим образом [1]: Обновитесь блок B1 с помощью Гастингса Столицы продвигаются, чтобы сгенерировать матрицу перехода, блок B2 обновления, состояния номера копии с помощью стохастического форварда распространяют назад производящий алгоритм, блок B3 обновления путем вычисления mus и блока B4 обновления, чтобы сгенерировать сигмы.
for iloop = 2:NMC % Compute the number of transitions from state i to state j for i =1:NS for j = 1:NS st_counts(i, j) = sum((states(1:N-1, iloop-1) == i) .* (states(2:N, iloop-1) == j)); end end % Updating block B1 % Generate the transition matrix from the Dirichlet distributions C = acghhmmsample('dirichlet', st_counts + 1, NS); % Compute the state probabilities under stationary distribution of a % given transition matrix C. PiC = PI(C, NS); % Compute the accepting probability using a Metropolis-Hastings step beta = min([1, exp(log(PiC(states(1, iloop-1))) - log(Pi(states(1, iloop-1))))]); if rand < beta A = C; Pi = PiC; end % Updating block B2 % Generate copy number states using Forward propagate, backward sampling [4]. states(:, iloop) = acghhmmfb(Pi, A, B); % Updating blocks B3 and B4 for i = 1:NS idx_s = states(:, iloop) == i; num_states = sum(idx_s); % If state i is not observed, then draw from its prior parameters if num_states == 0 mus(i, iloop) = acghhmmsample('normal', mus_hyper(i),... taus_hyper(i), mu_low_bounds(i), mu_up_bounds(i)); sigmas(i, iloop)= acghhmmsample('gamma', sg_alpha(i),... sg_beta(i), sg_bounds(i)); else Y_avg = mean(Y(idx_s)); theta_prec = 1/taus_hyper(i)^2 + num_states/sigmas(i,iloop-1)^2; weight_means = (mus_hyper(i)/(taus_hyper(i)^2) +... Y_avg * num_states/(sigmas(i, iloop-1)^2))/theta_prec; % Compute mus - B3 mus(i, iloop) = acghhmmsample('normal', weight_means, ... 1/sqrt(theta_prec), mu_low_bounds(i), mu_up_bounds(i)); % Compute sigmas - B4 Y_v = sum((Y(idx_s) - mus(i, iloop)).^2); sigmas(i, iloop) = acghhmmsample('gamma', sg_alpha(i)+num_states/2,... sg_beta(i)+Y_v/2, sg_bounds(i)); end % Update the emission matrix with new mus and sigmas. B(i,:) = normpdf(Y, mus(i,iloop),sigmas(i,iloop)); end end
Постройте следующее среднее значение mu распределения четырех состояний.
figure; for j = 1:NS subplot(2,2,j) ksdensity(mus(j,:)); title(sprintf('State %d', j)) end sgtitle('Distribution of Mu of States'); hold off;
Постройте следующие распределения сигмы четырех состояний.
figure; for j = 1:NS subplot(2,2,j) ksdensity(sigmas(j,:)); title(sprintf('State %d', j)) end sgtitle('Distribution of Sigma of States'); hold off;
Чертите метку состояния для каждого клона от выборки MCMC и вычислите апостериорные вероятности каждого состояния.
clone_states = zeros(1, N); state_prob = zeros(NS, N); state_count = zeros(NS, N); for i = 1:N % for each clone state = states(i, :); for j=1:NS state_count(j, i) = sum(state == j); end selState = find(state_count(:,i) == max(state_count(:,i))); if length(selState) > 1 if i ~= 1 clone_states(i) = clone_states(i-1); else clone_states(i) = min(selState); end else clone_states(i) = selState; end state_prob(:, i) = state_count(:,i)/NMC; end clone_states = clone_states';
Постройте метку состояния для каждого клона на хромосоме 12 из демонстрационного PA.C.Dan.G.
figure; leg = zeros(1,4); for i = 1:N if clone_states(i) == 1 leg(1) = plot(i,Y(i),'v', 'MarkerFaceColor', [1 0.2 0.2],... 'MarkerEdgeColor', 'none'); elseif clone_states(i) == 2 leg(2) = plot(i,Y(i),'o', 'Color', [0.4 0.4 0.4]); elseif clone_states(i) == 3 leg(3) = plot(i,Y(i),'^', 'MarkerFaceColor', [0.2 1 0.2],... 'MarkerEdgeColor', 'none'); elseif clone_states(i) == 4 leg(4) = plot(i, Y(i), '^', 'MarkerFaceColor', [0.2 0.2 1],... 'MarkerEdgeColor', 'none'); end hold on; end ylim(gca, ylims) legend(leg, 'State 1', 'State 2','State 3','State 4') xlabel('Index') ylabel('Log2(ratio)') title('State Label') hold off
Поскольку каждый MCMC чертит, сгенерированные состояния могут быть классифицированы как фокальные аберрации, точки перехода, усиления, выбросы и целые хромосомные изменения [1]. В этом примере вы найдете высокоуровневые усиления, точки перехода и выбросы на хромосоме 12 из демонстрационного PA.C.Dan.G.
Клон с состоянием = 4 рассматривается высокоуровневым усилением [1]. Найдите высокоуровневые усиления.
high_lvl_amp_idx = find(clone_states == 4);
Точка перехода сопоставлена с крупномасштабными областями прибылей и убытков и объявляется, когда ширина измененной области превышает 5 мега пар оснований [1]. Найдите точки перехода.
region_lim = 5e6; focalabr_idx=[1;find(diff(clone_states)~=0);N]; istranspoint = false(length(focalabr_idx), 1); for i = 1:length(focalabr_idx)-1 region_x = X(focalabr_idx(i+1)) - X(focalabr_idx(i)); istranspoint(i+1) = region_x > region_lim; end trans_idx = focalabr_idx(istranspoint); % Remove adjacent trans_idx that have the same states. hasadjacentstate = [diff(clone_states(trans_idx))==0; true]; trans_idx = trans_idx(~hasadjacentstate) focalabr_idx = focalabr_idx(~istranspoint); focalabr_idx = focalabr_idx(2:end-1);
trans_idx = 107 135 323
Выброс для усилений является фокальной аберрацией, удовлетворяющей ее z-счету, больше, чем 2, в то время как выброс за потери имеет z-счет меньше чем-2 [1].
Найдите выбросы за потери
outlier_loss_idx = focalabr_idx(clone_states(focalabr_idx) == 1) if ~isempty(outlier_loss_idx) [F,Xi] = ksdensity(mus(1,:)); [dummy, idx] = max(F); mu_1 = Xi(idx); [F,Xi] = ksdensity(sigmas(1,:)); [dummy, idx] = max(F); sigma_1 = Xi(idx); outlier_loss_idx = outlier_loss_idx((Y(outlier_loss_idx) - mu_1)/sigma_1 < -2) end
outlier_loss_idx = 0x1 empty double column vector
Найдите выбросы для усилений
outlier_gain_idx = focalabr_idx(clone_states(focalabr_idx) == 3); if ~isempty(outlier_gain_idx) [F,Xi] = ksdensity(mus(3,:)); [dummy, idx] = max(F); mu_1 = Xi(idx); [F,Xi] = ksdensity(sigmas(3,:)); [dummy, idx] = max(F); sigma_1 = Xi(idx); outlier_gain_idx = outlier_gain_idx((Y( outlier_gain_idx) - mu_1)/sigma_1 > 2) end
outlier_gain_idx = 0x1 empty double column vector
Добавьте классифицированные метки в график отношения интенсивности хромосомы 12 из демонстрационного PA.C.Dan.G. Постройте средние значения сегмента из процедуры CBS для сравнения.
figure; hl1 = plot(X, Y, '.', 'color', [0.4 0.4 0.4]); hold on; if ~isempty(high_lvl_amp_idx) hl2 = line(X(high_lvl_amp_idx), Y(high_lvl_amp_idx),... 'LineStyle', 'none',... 'Marker', '^',... 'MarkerFaceColor', [0.2 0.2 1],... 'MarkerEdgeColor', 'none'); end if ~isempty(trans_idx) for i = 1:numel(trans_idx) hl3 = line(ones(1,2)*X(trans_idx(i)), [-3.5, 3.5],... 'LineStyle', '--',... 'Color', [1 0.6 0.2]); end end if ~isempty(outlier_gain_idx) line(X(outlier_gain_idx), Y(outlier_gain_idx),... 'LineStyle', 'none',... 'Marker', 'v',... 'MarkerFaceColor', [1 0 0],... 'MarkerEdgeColor', 'none'); end if ~isempty(outlier_loss_idx) hl4 = line(X(outlier_loss_idx), Y(outlier_loss_idx),... 'LineStyle', 'none',... 'Marker', 'v',... 'MarkerFaceColor', [1 0 0],... 'MarkerEdgeColor', 'none'); end % Plot segment means from the CBS procedure. for i = 1:numel(PS.SegmentData.Start) hl5 = line([PS.SegmentData.Start(i) PS.SegmentData.End(i)],... [PS.SegmentData.Mean(i) PS.SegmentData.Mean(i)],... 'Color', [1 0 0],... 'LineWidth', 1.5); end ylim(gca, ylims) ylabel('Log2(Ratio)') title(sprintf('%s - Chromosome %d', sample, chromID)) % Adding chromosome 12 ideogram and legends to the plot. chromosomeplot('hs_cytoBand.txt', chromID, 'addtoplot', gca) legend([hl1, hl2, hl3,hl5], 'IntensityRatio', 'Amplification',... 'TransitionPoint', 'CBS SegmentMean')
Байесов алгоритм HMM нашел 3 точки перехода обозначенными поврежденными вертикальными линиями в графике. Байесов алгоритм HMM идентифицировал, что два высоких уровня усилили области, отмеченные синими-треугольниками в графике. Усиленные области двух высоких уровней соответствуют двум минимальным общим областям (MCRs) [2] на хромосоме 12, сопоставленный с усилениями номера копии, как объяснил Агирре и др., (2004). Байесов HMM объявил первый набор порций высокой интенсивности как одна область высокоуровневого усиления. В сравнении процедура CBS не удалась обнаружить второй MCR и сегментировала первый MCR на две области. Никакой выброс не был обнаружен в этом примере.
[1] Guha, S., Литий, Y. и Neuberg, D., "Байесово скрытое моделирование Маркова массива данные CGH", Журнал американской Статистической Ассоциации, 103 (482):485-497, 2008.
[2] Агирре, A.J., и др., "Характеристика с высоким разрешением генома аденокарциномы поджелудочной железы", PNAS, 101 (24):9067-72, 2004.
[3] Olshen, A.B., и др., "Круговая бинарная сегментация для анализа основанной на массиве ДНК копирует данные о номере", Биостатистика, 5 (4):557-7, 2004.
[4] Шах, S.P., и др., "Интегрируя полиморфизмы номера копии в массив анализ CGH с помощью устойчивого HMM", Биоинформатика, 22 (14): e431-e439, 2006