Этот пример показывает, как использовать Байесов метод скрытой модели Маркова (HMM), чтобы идентифицировать изменение номера копии в основанных на массиве данных о сравнительной геномной гибридизации (CGH).
Основанный на массиве CGH является методом высокой пропускной способности, чтобы измерить изменение номера копии DNA через геном. Фрагменты DNA или "клоны" теста и ссылочных выборок гибридизируются к сопоставленным фрагментам массивов. Отношения интенсивности 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 является неизвестным параметром для каждого состояния с этим ограничением:
Уголовное прошлое для средних изменений номера копии:
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 suptitle('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 suptitle('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., и др., "Круговая бинарная сегментация для анализа основанного на массиве DNA копирует данные о номере", Биостатистика, 5 (4):557-7, 2004.
[4] Шах, S.P., и др., "Интегрируя полиморфизмы номера копии в массив анализ CGH с помощью устойчивого HMM", Биоинформатика, 22 (14): e431-e439, 2006