В этом примере показано, как использовать байесовский метод скрытой марковской модели (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 et al., (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_eq02646086227142414759.png)
Guha et al., (2006) также описал алгоритм Metropolis-with-Gibbs для генерации апостериорных выборок. Алгоритм MCMC независимо запускается для каждой хромосомы, чтобы сгенерировать выборку MCMC для параметров хромосомы. Начальные значения параметров генерируются из априоров. Сгенерированные состояния числа копий представляют ничьи из маргинального апостериора интереса, Для каждого розыгрыша MCMC сгенерированные состояния проверяются и классифицируются как фокальные абляции, точки перехода, усиление, выбросы и целые хромосомные изменения.
В этом примере вы примените алгоритм Bayesian HMM для анализа профилей CGH массива некоторых выборок рака поджелудочной железы [2].
Данные в этом примере являются профилями массива CGH 24 клеточных линий аденокарциномы поджелудочной железы и 13 первичных образцов опухоли от Alguirre et al., (2004). Меченные фрагменты ДНК гибридизуют с микромассивами кДНК человека Agilent ®, содержащими 14160 кДНК-клонов. Около 9420 клонов имеют уникальные положения карты с медианным интервалом между отображенными элементами 100,1 кб. Более подробную информацию о данных и эксперименте можно найти в [2]. Для удобства данные уже нормированы, и основанные на логах коэффициенты интенсивности 2 предоставлены файлом MAT pancrea_oligocgh.mat.
Вы примените алгоритм Bayesian HMM, чтобы проанализировать хромосому 12 выборки 6 данных аденокарциномы поджелудочной железы и сравнить результаты с сегментами, найденными алгоритмом круговой двоичной сегментации (CBS) Olshen et al., (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 использует алгоритм Metropolis-with-Gibbs, чтобы сгенерировать апостериорные выборки параметров [1]. Параметры модели сгруппированы в четыре блока. Алгоритм итерационно генерирует каждый из четырех блоков, обусловленных оставшимися блоками и данными.
Чтобы проанализировать данные с помощью алгоритма Bayesian HMM, необходимо инициализировать параметры. Более подробную информацию о предыдущих параметрах можно найти в ссылках [1] и [4].
Инициализируйте состояние генератора случайных чисел, чтобы убедиться, что рисунки, сгенерированная этими командами, совпадают с рисунками в HTML версии этого примера.
rng('default');
Определите количество состояний
NS = 4;
Определите количество итераций MCMC
NMC = 100;
Определите гиперпараметры предыдущих распределений для четырех состояний.
mus_hyper = [-1, 0, 0.58, 1]; taus_hyper = [1, 1, 1, 2];
Установите параметр epsilon, который определяет ограничения средств.
eps = 0.1;
Установите границы предшествующих средств каждого состояния.
mu_low_bounds = [-Inf, -eps, eps, 0.58]; mu_up_bounds = [-eps, eps, 0.58, Inf];
Guha et al., (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 из распределений Дирихле.
a = ones(NS, NS);
A = acghhmmsample('dirichlet', a, NS);
Переходная матрица имеет уникальное стационарное распределение. Стационарное распределение PI является собственным вектором переходной матрицы, сопоставленной с собственным значением 1.
PI =@(x, n) (ones(1,n)/(eye(n) -x + ones(n)))';
Сгенерируйте предыдущий ПИ стационарного распределения.
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 с помощью шага Metropolis-Hastings, чтобы сгенерировать матрицу перехода, обновить блок B2 состояниями числа копий с помощью стохастического алгоритма прямого распространения дискретизации назад, обновить B3 блоков путем вычисления mus и обновить B4 блоков, чтобы сгенерировать sigmas
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')
Алгоритм Bayesian HMM обнаружил 3 точки перехода, обозначенные ломаными вертикальными линиями на графике. Алгоритм Bayesian HMM определил две высокоуровневые усиленные области, отмеченные синими вверх-треугольниками на графике. Две высокоуровневые усиленные области соответствуют двум минимальным общим областям (MCR) [2] на хромосоме 12, связанным с увеличением числа копий, как объяснено Aguirre et al., (2004). Байесовская HMM объявила первый набор рационов высокой интенсивности как одну область высокоуровневой амплификации. Для сравнения, процедура CBS не обнаружила второй MCR и сегментировала первый MCR на две области. В этом примере не было обнаружено выбросов.
[1] Guha, S., Li, Y. and Neuberg, D., «Bayesian hidden Markov моделирования of массива CGH данных», Journal of the American Statistical Association, 103 (482): 485-497, 2008.
[2] Aguirre, A.J., et al., «High-resolution charcerization of the pancreatic adenocarcinoma genome», PNAS, 101 (24): 9067-72, 2004.
[3] Olshen, A.B., et al., «Circular binary segmentation for the analysis of array-based ДНК copy number данных», Biostatistics, 5 (4): 557-7, 2004.
[4] Shah, S.P., et al., «Интеграция полиморфизмов числа копий в анализ CGH массива с помощью устойчивой HMM», Bioinformatics, 22 (14): e431-e439, 2006