Введение

Основанный на массивах 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 является неизвестным параметром для каждого состояния с этим ограничением:

Априорные значения для изменения среднего числа копий:

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

Классификация профилей CGH массива

Для каждого розыгрыша 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