Exponenta Banner
exponenta event banner

Анализ данных CGH на основе массива с использованием байесовского моделирования скрытого Маркова

Открыть сценарий

В этом примере показано, как использовать метод байесовской скрытой модели Маркова (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 значительно отличаются от теоретических значений. Отношения обычно уменьшаются до нуля. Одним из основных факторов является загрязнение образцов опухоли нормальными клетками. Существует также зависимость между соотношениями интенсивности соседних клонов. Эти проблемы требуют использования эффективных статистических алгоритмов, характеризующих геномные профили.

Байесовский HMM

Guha et al., (2006) предложил байесовский статистический подход в зависимости от скрытой модели Маркова (HMM) для анализа данных массива CGH. На скрытую марковскую модель приходится зависимость между соседними клонами. Отношения интенсивности генерируются скрытыми состояниями числа копий. Байесовское обучение используется для идентификации общегеномных изменений в количестве копий из данных. Задние выводы сделаны относительно прироста и потерь числа копий.

В этом байесовском алгоритме HMM существуют четыре состояния, определяемые как состояние потери числа копий (1), состояние нейтрали числа копий (2), состояние усиления одной копии (3) и состояние усиления (множественное усиление) (4). С каждым клоном связано состояние номера копии. Нормализованные отношения log2 считаются распределенными как

$$Y_k \sim N(\mu_{sk}, \sigma_{sk}^2)$$

mu является неизвестным параметром для каждого состояния с этим ограничением:

$$\mu_1<\mu_2<\mu_3<\mu.4$$

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

$$\mu_1 \sim N \mathrm(-1,\tau_1^2) \cdot I\mathrm(\mu_1< -\epsilon)$$

$$\mu_2 \sim N(0,\tau_2^2) \cdot I(-\epsilon<\mu_2<\epsilon)$$

$$\mu_3 \sim N(0.58,\tau_3^2) \cdot I(\epsilon<\mu_3<0.58)$$

$$[\mu_4|\mu_3, \sigma_3] \sim N(1,\tau_4^2) \cdot I(\mu_4&#62;\mu_3 +&#xA;3\sigma_3)$$

Guha et al., (2006) также описал алгоритм Метрополиса-внутри-Гиббса для генерации задних образцов. Алгоритм MCMC независимо выполняется для каждой хромосомы для генерации образца MCMC для параметров хромосомы. Начальные значения параметров генерируются на основе предварительных значений. Сгенерированные состояния числа копий представляют собой розыгрыши из краевого заднего конца, представляющего интерес. Для каждого розыгрыша MCMC генерируемые состояния проверяются и классифицируются как фокальные абляции, точки перехода, амплификации, отклонения и целые хромосомные изменения.

В этом примере для анализа матричных профилей CGH некоторых образцов рака поджелудочной железы применяется байесовский алгоритм HMM [2].

Загрузка данных

Данные в этом примере являются матричными CGH-профилями 24 клеточных линий аденокарциномы поджелудочной железы и 13 первичных образцов опухолей от Alguirre et al., (2004). Меченые фрагменты ДНК гибридизовали с микрочипами кДНК человека Agilent ®, содержащими 14 160 клонов кДНК. Около 9420 клонов имеют уникальные положения карты со средним интервалом между отображенными элементами 100,1 кб. Подробнее с данными и экспериментом можно ознакомиться в [2]. Для удобства данные уже нормализованы, и отношения интенсивности на основе log2 обеспечиваются файлом MAT.pancrea_oligocgh.mat.

Вы примените байесовский алгоритм 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 использует алгоритм Метрополиса-внутри-Гиббса для генерации задних выборок параметров [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];

Задайте параметр 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 = 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 путем вычисления муса и блок обновления 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

Классификация профилей 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')

Байесовский алгоритм HMM обнаружил на графике 3 точки перехода, обозначенные разорванными вертикальными линиями. Байесовский алгоритм HMM идентифицировал на графике две высокоуровневые усиленные области, отмеченные синими восходящими треугольниками. Две усиленные области высокого уровня соответствуют двум минимальным общим областям (MCR) [2] на хромосоме 12, связанным с увеличением числа копий, как объясняется Aguirre et al., (2004). Байесовский ПММ объявил первый набор высокоинтенсивных пайков единой областью высокоуровневого усиления. Для сравнения, процедура CBS не смогла обнаружить второй MCR и сегментировала первый MCR на две области. В этом примере отклонение не обнаружено.

Ссылки

[1] Гуха, С., Ли, Я. и Нойберг, Д., «Байесовское скрытое марковское моделирование данных массива CGH», Journal of the American Statistical Association, 103 (482): 485-497, 2008.

[2] Агирре, A.J., et al., «Характеристика высокого разрешения генома аденокарциномы поджелудочной железы», PNAS, 101 (24): 9067-72, 2004.

[3] Olshen, A.B., et al., «Круговая бинарная сегментация для анализа данных числа копий ДНК на основе массива», Biostatistics, 5 (4): 557-7, 2004.

[4] Shah, S.P., et al., «Интеграция полиморфизмов числа копий в анализ CGH массива с использованием надежной HMM», Bioinformatics, 22 (14): e431-e439, 2006

Документация по инструментарию биоинформатики

Поддержка