Введение

Основанный на массиве 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 и др., (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

Классификация массива профили 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 идентифицировал, что два высоких уровня усилили области, отмеченные синими-треугольниками в графике. Усиленные области двух высоких уровней соответствуют двум минимальным общим областям (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