В этом примере показано, как обнаружить изменения номера копии ДНК в основанных на массиве данных о сравнительной геномной гибридизации (CGH) всего генома.
Скопируйте изменения номера, или изменения форма наследственной изменчивости в геноме человека [1]. Изменения номера копии ДНК (CNAs) были соединены с разработкой и развитием рака и многих болезней.
ДНК основанная на микромассиве сравнительная геномная гибридизация (CGH) является методом, позволяет одновременный контроль количества копии тысяч генов в геноме [2,3]. В этом методе фрагменты ДНК или "клоны" от тестовой выборки и ссылочной выборки, дифференцированно помеченной красками (обычно, Cy3 и Cy5), гибридизируются к сопоставленным микромассивам ДНК и отображаются. Изменения номера копии связаны с Cy3 и отношением интенсивности флюоресценции Cy5 целей, гибридизированных к каждому зонду на микромассиве. Клоны с нормированной тестовой интенсивностью, значительно больше, чем ссылочная интенсивность, указывают на усиления номера копии в тестовой выборке в тех положениях. Точно так же значительно более низкая интенсивность в тестовой выборке является знаками потери номера копии. BAC (бактериальная искусственная хромосома) клон базирующиеся массивы CGH имеет разрешение в порядке одного миллиона пар оснований (1 МБ) [3]. Олигонуклеотид и массивы комплементарной ДНК обеспечивают более высокое разрешение 50-100kb [2].
Массив CGH находящиеся в log2 отношения интенсивности предоставляет полезную информацию о CNAs всего генома. В людях нормальный номер копии ДНК два для всех аутосом. В идеальной ситуации нормальные клоны соответствовали бы log2 отношению нуля. log2 отношения интенсивности одной потери копии были бы-1, и одно усиление копии будет 0.58. Цель состоит в том, чтобы эффективно идентифицировать местоположения прибылей или убытков номера копии ДНК.
Данные в этом примере являются массивом BAC клеточной линии Coriell данные CGH, анализируемые Snijders и др. (2001). Данные о клеточной линии Coriell широко рассматриваются как набор данных "золотого стандарта". Можно загрузить эти данные нормированных находящихся в log2 отношений интенсивности и дополнительную таблицу известных кариотипов от https://www.nature.com/articles/ng754#supplementary-information. Вы сравните эти cytogenically сопоставленные изменения с местоположениями прибылей или убытков, идентифицированных с различными функциями MATLAB и его тулбоксов.
В данном примере данные о клеточной линии Coriell обеспечиваются в файле MAT. Файл данных coriell_baccgh.mat содержит coriell_data, структура, содержащая нормированного среднего значения находящегося в log2 теста к ссылочным отношениям интенсивности 15 клеточных линий фибробласта и их геномных положений. Цели BAC упорядочены положением генома, начинающимся в 1 пункте и заканчивающимся в Xq.
load coriell_baccgh
coriell_data
coriell_data =
struct with fields:
Sample: {1x15 cell}
Chromosome: [2285x1 int8]
GenomicPosition: [2285x1 int32]
Log2Ratio: [2285x15 double]
FISHMap: {2285x1 cell}
Можно построить геном широкие находящиеся в log2 отношения интенсивности теста/ссылки клонов ДНК. В этом примере вы отобразите log2 отношения интенсивности для клеточной линии GM03576 для хромосом 1 - 23.
Найдите демонстрационный индекс для клеточной линии CM03576.
sample = find(strcmpi(coriell_data.Sample, 'GM03576'))
sample =
8
Чтобы пометить хромосомы и чертить границы хромосомы, необходимо найти количество точек данных в каждой хромосоме.
chr_nums = zeros(1, 23); chr_data_len = zeros(1,23); for c = 1:23 tmp = coriell_data.Chromosome == c; chr_nums(c) = find(tmp, 1, 'last'); chr_data_len(c) = length(find(tmp)); end % Draw a vertical bar at the end of a chromosome to indicate the border x_vbar = repmat(chr_nums, 3, 1); y_vbar = repmat([2;-2;NaN], 1, 23); % Label the autosomes with their chromosome numbers, and the sex chromosome % with X. x_label = chr_nums - ceil(chr_data_len/2); y_label = zeros(1, length(x_label)) - 1.6; chr_labels = num2str((1:1:23)'); chr_labels = cellstr(chr_labels); chr_labels{end} = 'X'; figure hold on h_ratio = plot(coriell_data.Log2Ratio(:,sample), '.'); h_vbar = line(x_vbar, y_vbar, 'color', [0.8 0.8 0.8]); h_text = text(x_label, y_label, chr_labels,... 'fontsize', 8, 'HorizontalAlignment', 'center'); h_axis = h_ratio.Parent; h_axis.XTick = []; h_axis.YGrid = 'on'; h_axis.Box = 'on'; xlim([0 chr_nums(23)]) ylim([-1.5 1.5]) title(coriell_data.Sample{sample}) xlabel({'', 'Chromosome'}) ylabel('Log2(T/R)') hold off
В графике границы между хромосомами обозначаются серыми вертикальными панелями. График показывает, что клеточная линия GM03576 является trisomic для хромосом 2 и 21 [3].
Можно также построить профиль каждой хромосомы в геноме. В этом примере вы отобразите log2 отношения интенсивности для каждой хромосомы в клеточной линии GM05296 индивидуально.
sample = find(strcmpi(coriell_data.Sample, 'GM05296')); figure; for c = 1:23 idx = coriell_data.Chromosome == c; chr_y = coriell_data.Log2Ratio(idx, sample); subplot(5,5,c); hp = plot(chr_y, '.'); line([0, chr_data_len(c)], [0,0], 'color', 'r'); h_axis = hp.Parent; h_axis.XTick = []; h_axis.Box = 'on'; xlim([0 chr_data_len(c)]) ylim([-1.5 1.5]) xlabel(['chr ' chr_labels{c}], 'FontSize', 8) end sgtitle('GM05296');
График показывает, что клеточная линия GM05296 имеет частичную трисомию в хромосоме 10 и частичная моносомия в хромосоме 11.
Заметьте, что прибыли и убытки номера копии дискретны. Эти изменения происходят в непрерывных областях хромосомы, которые покрывают несколько клонов, чтобы назвать хромосому.
Основанные на массиве данные CGH могут быть довольно шумными. Поэтому точная идентификация областей хромосомы равного номера копии, который составляет шум в данных, требует устойчивых вычислительных методов. В остальной части этого примера вы будете работать с данными хромосом 9, 10 и 11 из клеточной линии GM05296.
Инициализируйте массив структур для данных этих трех хромосом.
GM05296_Data = struct('Chromosome', {9 10 11},... 'GenomicPosition', {[], [], []},... 'Log2Ratio', {[], [], []},... 'SmoothedRatio', {[], [], []},... 'DiffRatio', {[], [], []},... 'SegIndex', {[], [], []});
Простой подход, чтобы выполнить высокоуровневое сглаживание должен использовать непараметрический фильтр. Функциональный mslowess реализует линейную подгонку к выборкам в окне перемены, этот пример, вы используете SPAN из 15 выборок.
for iloop = 1:length(GM05296_Data) idx = coriell_data.Chromosome == GM05296_Data(iloop).Chromosome; chr_x = coriell_data.GenomicPosition(idx); chr_y = coriell_data.Log2Ratio(idx, sample); % Remove NaN data points idx = ~isnan(chr_y); GM05296_Data(iloop).GenomicPosition = double(chr_x(idx)); GM05296_Data(iloop).Log2Ratio = chr_y(idx); % Smoother GM05296_Data(iloop).SmoothedRatio = ... mslowess(GM05296_Data(iloop).GenomicPosition,... GM05296_Data(iloop).Log2Ratio,... 'SPAN',15); % Find the derivative of the smoothed ratio GM05296_Data(iloop).DiffRatio = ... diff([0; GM05296_Data(iloop).SmoothedRatio]); end
Чтобы лучше визуализировать и позже подтвердить местоположения изменений номера копии, нам нужна cytoband информация. Считайте человеческие cytoband информации из hs_cytoBand.txt файл данных с помощью cytobandread функция. Это возвращает структуру человеческой cytoband информации [4].
hs_cytobands = cytobandread('hs_cytoBand.txt') % Find the centromere positions for the chromosomes. acen_idx = strcmpi(hs_cytobands.GieStains, 'acen'); acen_ends = hs_cytobands.BandEndBPs(acen_idx); % Convert the cytoband data from bp to kilo bp because the genomic % positions in Coriell Cell Line data set are in kilo base pairs. acen_pos = acen_ends(1:2:end)/1000;
hs_cytobands =
struct with fields:
ChromLabels: {862x1 cell}
BandStartBPs: [862x1 int32]
BandEndBPs: [862x1 int32]
BandLabels: {862x1 cell}
GieStains: {862x1 cell}
Можно смотреть данные путем графического вывода находящихся в log2 отношений, сглаживавших отношений и производной сглаживавших отношений вместе. Можно также отобразить положение центромеры хромосомы в графиках данных. Пурпурная вертикальная панель отмечает центромеру хромосомы.
for iloop = 1:length(GM05296_Data) chr = GM05296_Data(iloop).Chromosome; chr_x = GM05296_Data(iloop).GenomicPosition; figure hold on plot(chr_x, GM05296_Data(iloop).Log2Ratio, '.'); line(chr_x, GM05296_Data(iloop).SmoothedRatio,... 'Color', 'r', 'LineWidth', 2); line(chr_x, GM05296_Data(iloop).DiffRatio,... 'Color', 'k', 'LineWidth', 2); line([acen_pos(chr), acen_pos(chr)], [-1, 1],... 'Color', 'm', 'LineWidth', 2, 'LineStyle', '-.'); if iloop == 1 legend('Raw','Smoothed','Diff', 'Centromere'); end ylim([-1, 1]) xlabel('Genomic Position') ylabel('Log2(T/R)') title(sprintf('GM05296: Chromosome %d ', chr)) hold off end
Производные сглаживавшего отношения по определенному порогу обычно указывают на существенные изменения с большим peaks и обеспечивают оценку индексов точки перехода. Для этого примера вы выберете порог 0,1.
thrd = 0.1; for iloop = 1:length(GM05296_Data) idx = find(abs(GM05296_Data(iloop).DiffRatio) > thrd ); N = numel(GM05296_Data(iloop).SmoothedRatio); GM05296_Data(iloop).SegIndex = [1;idx;N]; % Number of possible segments found fprintf('%d segments initially found on Chromosome %d.\n',... numel(GM05296_Data(iloop).SegIndex) - 1,... GM05296_Data(iloop).Chromosome) end
1 segments initially found on Chromosome 9. 4 segments initially found on Chromosome 10. 5 segments initially found on Chromosome 11.
Гауссова смесь (GM) или кластеризация Максимизации ожидания (EM) могут предоставить точные настройки индексам точки перехода [5]. Сходимость в статистически оптимальные индексы точки перехода может быть упрощена путем окружения каждого индекса с набором равной длины смежных индексов. Таким образом каждое ребро сопоставлено с левыми и правыми распределениями. Кластеризация GM изучает параметры наибольшего правдоподобия этих двух распределений. Это затем оптимально настраивает индексы, учитывая изученные параметры.
Можно установить длину для набора смежных положений, распределенных вокруг индексов точки перехода. В данном примере вы выберете длину 5. Можно также смотреть каждую точку перехода путем графического вывода ее кластеров GM. В этом примере вы построите кластеры GM для Хромосомы 10 данных.
len = 5; for iloop = 1:length(GM05296_Data) seg_num = numel(GM05296_Data(iloop).SegIndex) - 1; if seg_num > 1 % Plot the data points in chromosome 10 data if GM05296_Data(iloop).Chromosome == 10 figure hold on; plot(GM05296_Data(iloop).GenomicPosition,... GM05296_Data(iloop).Log2Ratio, '.') ylim([-0.5, 1]) xlabel('Genomic Position') ylabel('Log2(T/R)') title(sprintf('Chromosome %d - GM05296', ... GM05296_Data(iloop).Chromosome)) end segidx = GM05296_Data(iloop).SegIndex; segidx_emadj = GM05296_Data(iloop).SegIndex; for jloop = 2:seg_num ileft = min(segidx(jloop) - len, segidx(jloop)); iright = max(segidx(jloop) + len, segidx(jloop)); gmx = GM05296_Data(iloop).GenomicPosition(ileft:iright); gmy = GM05296_Data(iloop).SmoothedRatio(ileft:iright); % Select initial guess for the cluster index for each point. gmpart = (gmy > (min(gmy) + range(gmy)/2)) + 1; % Create a Gaussian mixture model object gm = gmdistribution.fit(gmy, 2, 'start', gmpart); gmid = cluster(gm,gmy); segidx_emadj(jloop) = find(abs(diff(gmid))==1) + ileft; % Plot GM clusters for the change-points in chromosome 10 data if GM05296_Data(iloop).Chromosome == 10 plot(gmx(gmid==1),gmy(gmid==1), 'g.',... gmx(gmid==2), gmy(gmid==2), 'r.') end end % Remove repeat indices zeroidx = [diff(segidx_emadj) == 0; 0]; GM05296_Data(iloop).SegIndex = segidx_emadj(~zeroidx); end % Number of possible segments found fprintf('%d segments found on Chromosome %d after GM clustering adjustment.\n',... numel(GM05296_Data(iloop).SegIndex) - 1,... GM05296_Data(iloop).Chromosome) end hold off;
1 segments found on Chromosome 9 after GM clustering adjustment. 3 segments found on Chromosome 10 after GM clustering adjustment. 5 segments found on Chromosome 11 after GM clustering adjustment.
Если вы определяете оптимальные индексы точки перехода, также необходимо определить, представляет ли каждый сегмент статистически существенные изменения в номере копии ДНК. Вы выполните t-тесты сочетания, чтобы оценить значение идентифицированных сегментов. Сегмент включает все точки данных от одной точки перехода до следующей точки перехода или конца хромосомы. В этом примере вы выполните 10 000 сочетаний точек данных на двух последовательных сегментах вдоль хромосомы на уровне значения 0,01.
alpha = 0.01; for iloop = 1:length(GM05296_Data) seg_num = numel(GM05296_Data(iloop).SegIndex) - 1; seg_index = GM05296_Data(iloop).SegIndex; if seg_num > 1 ppvals = zeros(seg_num+1, 1); for sloop = 1:seg_num-1 seg1idx = seg_index(sloop):seg_index(sloop+1)-1; if sloop== seg_num-1 seg2idx = seg_index(sloop+1):(seg_index(sloop+2)); else seg2idx = seg_index(sloop+1):(seg_index(sloop+2)-1); end seg1 = GM05296_Data(iloop).SmoothedRatio(seg1idx); seg2 = GM05296_Data(iloop).SmoothedRatio(seg2idx); n1 = numel(seg1); n2 = numel(seg2); N = n1+n2; segs = [seg1;seg2]; % Compute observed t statistics t_obs = mean(seg1) - mean(seg2); % Permutation test iter = 10000; t_perm = zeros(iter,1); for i = 1:iter randseg = segs(randperm(N)); t_perm(i) = abs(mean(randseg(1:n1))-mean(randseg(n1+1:N))); end ppvals(sloop+1) = sum(t_perm >= abs(t_obs))/iter; end sigidx = ppvals < alpha; GM05296_Data(iloop).SegIndex = seg_index(sigidx); end % Number segments after significance tests fprintf('%d segments found on Chromosome %d after significance tests.\n',... numel(GM05296_Data(iloop).SegIndex) - 1, GM05296_Data(iloop).Chromosome) end
1 segments found on Chromosome 9 after significance tests. 3 segments found on Chromosome 10 after significance tests. 4 segments found on Chromosome 11 after significance tests.
Цитогенетическое исследование указывает на клеточную линию, GM05296 имеет трисомию в 10q21-10q24 и моносомию в 11p12-11p13 [3]. Постройте средние значения сегмента этих трех хромосом по исходным данным с полужирными красными линиями и добавьте идеограммы хромосомы в графики с помощью chromosomeplot функция. Обратите внимание на то, что геномные положения в наборе данных клеточной линии Coriell находятся в парах оснований килограмма. Поэтому необходимо будет преобразовать cytoband данные от BP до килограмма BP при добавлении идеограмм в график.
for iloop = 1:length(GM05296_Data) figure; seg_num = numel(GM05296_Data(iloop).SegIndex) - 1; seg_mean = ones(seg_num,1); chr_num = GM05296_Data(iloop).Chromosome; for jloop = 2:seg_num+1 idx = GM05296_Data(iloop).SegIndex(jloop-1):GM05296_Data(iloop).SegIndex(jloop); seg_mean(idx) = mean(GM05296_Data(iloop).Log2Ratio(idx)); line(GM05296_Data(iloop).GenomicPosition(idx), seg_mean(idx),... 'color', 'r', 'linewidth', 3); end line(GM05296_Data(iloop).GenomicPosition, GM05296_Data(iloop).Log2Ratio,... 'linestyle', 'none', 'Marker', '.'); line([acen_pos(chr_num), acen_pos(chr_num)], [-1, 1],... 'linewidth', 2,... 'color', 'm',... 'linestyle', '-.'); ylabel('Log2(T/R)') ax = gca; ax.Box = 'on'; ylim([-1, 1]) title(sprintf('Chromosome %d - GM05296', chr_num)); chromosomeplot(hs_cytobands, chr_num, 'addtoplot', gca, 'unit', 2) end
Как показано в графиках, никакие изменения номера копии не были найдены на хромосоме 9, существует промежуток усиления номера копии от 10q21 до 10q24, и область номера копии потерь от 11p12 до 11p13. CNAs, найденный соответствием известные результаты в клеточной линии GM05296, определяется цитогенетическим анализом.
Можно также отобразиться, CNAs клеточной линии GM05296 выравниваются к представлению сводных данных идеограммы хромосомы с помощью chromosomeplot функция. Определите геномные положения для CNAs на хромосомах 10 и 11.
chr10_idx = GM05296_Data(2).SegIndex(2):GM05296_Data(2).SegIndex(3)-1; chr10_cna_start = GM05296_Data(2).GenomicPosition(chr10_idx(1))*1000; chr10_cna_end = GM05296_Data(2).GenomicPosition(chr10_idx(end))*1000; chr11_idx = GM05296_Data(3).SegIndex(2):GM05296_Data(3).SegIndex(3)-1; chr11_cna_start = GM05296_Data(3).GenomicPosition(chr11_idx(1))*1000; chr11_cna_end = GM05296_Data(3).GenomicPosition(chr11_idx(end))*1000;
Создайте структуру, содержащую данные об изменении номера копии из данных о клеточной линии GM05296 согласно входным требованиям chromosomeplot функция.
cna_struct = struct('Chromosome', [10 11],... 'CNVType', [2 1],... 'Start', [chr10_cna_start, chr11_cna_start],... 'End', [chr10_cna_end, chr11_cna_end])
cna_struct =
struct with fields:
Chromosome: [10 11]
CNVType: [2 1]
Start: [69209000 34420000]
End: [105905000 35914000]
chromosomeplot(hs_cytobands, 'cnv', cna_struct, 'unit', 2) title('Human Karyogram with Copy Number Alterations of GM05296')
В этом примере показано, как MATLAB и его тулбоксы обеспечивают инструменты для анализа и визуализации изменений номера копии в основанных на массиве данных CGH.
[1] Редон, R., и др., "Глобальное изменение номера копии в геноме человека", Природа, 444 (7118):444-54, 2006.
[2] Pinkel, D., и др., "Анализ высокого разрешения ДНК копирует изменения номера с помощью сравнительной геномной гибридизации для микромассивов", Генетика Природы, 20 (2):207-11, 1998.
[3] Snijders, Утра, и др., "Блок микромассивов для измерения всего генома ДНК копирует номер", Генетика Природы, 29 (3):263-4, 2001.
[4] Геном человека сборка NCBI 36.
[5] Майерс, C.L., и др., "Точное обнаружение анеуплоидий в массиве CGH и микроданные массива экспрессии гена", Биоинформатика, 20 (18):3533-43, 2004.