Этот пример показывает, как обнаружить изменения числа копий ДНК в данных сравнительной геномной гибридизации (CGH) на основе массива генома.
Изменения или изменения числа копий являются формой генетических изменений в геноме человека [1]. Изменения числа копий ДНК (CNAs) были связаны с развитием и прогрессированием рака и многих заболеваний.
ДНК микромассива основе сравнительная геномная гибридизация (CGH) является методом, позволяющей одновременно контролировать количество копий тысяч генов по всему геному [2,3]. В этом методе фрагменты ДНК или «клоны» из тестовой выборки и ссылки выборки дифференциально меченные красителями (обычно Cy3 и Cy5) гибридизуются для отображения микромассивов ДНК и визуализируются. Изменения числа копий связаны с Cy3 и Cy5 интенсивностью флуоресценции мишеней, гибридизованных с каждым зондом на микромассиве. Клоны с нормированной интенсивностью теста, значительно большей, чем эталонная интенсивность, указывают на увеличение числа копий в тестовой выборке в этих положениях. Точно так же значительно более низкие интенсивности в тестовой выборке являются признаками потери числа копий. Массивы CGH на основе клонов BAC (бактериальная искусственная хромосома) имеют разрешение порядка одного миллиона пар оснований (1Mb) [3]. Массивы олигонуклеотидов и кДНК обеспечивают более высокое разрешение 50-100kb [2].
Основанные на массиве коэффициенты интенсивности CGH 2 предоставляют полезную информацию о общегеномных CNA. У людей нормальное число копий ДНК составляет два для всех аутосомов. В идеальной ситуации нормальные клоны будут соответствовать соотношению log2 нуля. Коэффициенты интенсивности потери единственного экземпляра log2 будут -1, а коэффициент усиления единственного экземпляра составит 0,58. Цель состоит в том, чтобы эффективно идентифицировать местоположения добавок или потерь номера копии ДНК.
Данные в этом примере являются данными CGH массива BAC камеры Кориелла, проанализированными Snijders et al. (2001). Данные клеточной линии Кориелла широко рассматриваются как набор данных «золотого стандарта». Вы можете загрузить эти данные нормированных основанных на логах 2 коэффициентов интенсивности и дополнительной таблицы известных кариотипов из https://www.nature.com/articles/ng754#supplementary-information. Вы сравните эти цитогенно сопоставленные изменения с местоположениями добавок или потерь, идентифицированных различными функциями MATLAB и его тулбоксов.
В данном примере данные линии камер Кориелла предоставляются в файле MAT. Файл данных coriell_baccgh.mat содержит coriell_data, структуру, содержащую нормированное среднее значение основанного на логарифмике теста для ссылки на коэффициенты интенсивности 15 линий камеры фибробластов и их геномных положений. Цели BAC упорядочены по положению генома, начиная с 1p и заканчивая 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}
Можно построить график коэффициентов интенсивности теста/эталона ДНК клонов на основе genome wide 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 клеточная линия является трисомной для хромосом 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
Чтобы лучше визуализировать и позже подтвердить местоположение изменений номера копии, нам нужна информация о цитобанде. Считайте информацию о человеческом цитобанде из hs_cytoBand.txt файл данных с использованием cytobandread функция. Он возвращает структуру информации о цитобанде человека [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}
Можно просмотреть данные, построив график основанных на логах коэффициентов, сглаженных коэффициентов и производной сглаженных коэффициентов вместе. Можно также отобразить положение центромеры хромосомы на графиках данных. Пурпурный вертикальный стержень помечает центромеру хромосомы.
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.
Если вы определяете оптимальные индексы точек изменения, вам также нужно определить, представляет ли каждый сегмент статистически значимые изменения в количестве копий ДНК. Вы выполните сочетания, чтобы оценить значимость идентифицированных сегментов. Сегмент включает все точки данных от одной точки изменения до следующей точки изменения или конца хромосомы. В этом примере вы выполните 10000 сочетания точек данных на двух последовательных сегментах вдоль хромосомы на уровне значимости 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 указаны в килограммах основы парах. Поэтому при добавлении идеограмм к графику вам нужно будет преобразовать данные цитобэнда из bp в kilo 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. Найденные CNA совпадают с известными результатами в клеточной линии GM05296 определяемыми цитогенетическим анализом.
Можно также отобразить CNA клеточной линии GM05296, выровненной с суммарным представлением хромосомной идеограммы, используя chromosomeplot функция. Определите геномные положения для CNA на хромосомах 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] Redon, R. et al., «Global изменения in copy number in the human genome», Nature, 444 (7118): 444-54, 2006.
[2] Pinkel, D., et al., «High разрешения analysis of ДНК copy number изменений with comparative genomic hybridization to микромассивов», Nature Genetics, 20 (2): 207-11, 1998.
[3] Snijders, A.M., et al., «Assembly of микромассивов for genome-wide измерения of ДНК copy number», Nature Genetics, 29 (3): 263-4, 2001.
[4] Геном человека NCBI Build 36.
[5] Myers, C.L., et al., «Accurate detection of aneuploidies in массива CGH and gene expression микромассива данных», Bioinformatics, 20 (18): 3533-43, 2004.