Этот пример показывает, как обнаружить изменения числа копий ДНК в данных сравнительной геномной гибридизации (CGH) на основе массива по всему геному.
Изменения или изменения числа копий - это форма генетической вариации в геноме человека [1]. Изменения числа копий ДНК (CNA) были связаны с развитием и прогрессированием рака и многих заболеваний.
Метод сравнительной геномной гибридизации (CGH) на основе микрочипов ДНК позволяет осуществлять одновременный мониторинг количества копий тысяч генов по всему геному [2,3]. В этом методе фрагменты ДНК или «клоны» из тестируемого образца и эталонного образца, дифференциально меченные красителями (обычно Cy3 и Cy5), гибридизуются с картированными микрочипами ДНК и изображаются. Изменения числа копий связаны с отношением интенсивности флуоресценции Cy3 и Cy5 мишеней, гибридизированных с каждым зондом на микрочипе. Клоны с нормализованной интенсивностью тестирования, значительно превышающей эталонную интенсивность, указывают увеличение числа копий в тестируемом образце в этих положениях. Аналогично, значительно меньшие интенсивности в тестируемом образце являются признаками потери числа копий. Матрицы CGH на основе клонов BAC (бактериальной искусственной хромосомы) имеют разрешение порядка одного миллиона пар оснований (1Mb) [3]. Олигонуклеотидные и кДНК массивы обеспечивают более высокое разрешение 50-100 кб [2].
Соотношения интенсивности на основе массива CGH log2 предоставляют полезную информацию об общегеномных CNA. У людей нормальное число копий ДНК составляет два для всех аутосом. В идеальной ситуации нормальные клоны соответствуют нулевому отношению log2. Отношения интенсивности log2 потери единичной копии будут равны -1, а коэффициент усиления единичной копии будет равен 0,58. Цель состоит в том, чтобы эффективно идентифицировать места получения или потери количества копий ДНК.
Данные в этом примере представляют собой данные CGH матрицы BAC клеточной линии Coriell, проанализированные Snijders et al. (2001). Данные клеточной линии Кориелла широко рассматриваются как набор данных «золотого стандарта». Вы можете загрузить эти данные нормализованных отношений интенсивности на основе log2 и дополнительную таблицу известных кариотипов из https://www.nature.com/articles/ng754#supplementary-information. Вы сравните эти цитогенно отображенные изменения с местоположениями прибылей или убытков, идентифицированными с различными функциями MATLAB и его панелей инструментов.
Для этого примера данные линии клеток Coriell представлены в файле MAT. Файл данных coriell_baccgh.mat содержит coriell_dataструктура, содержащая нормализованное среднее отношение теста на основе log2 к эталонным соотношениям интенсивности 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}
Вы можете построить график отношений тест/эталонная интенсивность ДНК-клонов на основе генома в широком диапазоне 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}
Можно проверить данные, выведя на график коэффициенты на основе 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
Производные сглаженного отношения по определенному порогу обычно указывают на существенные изменения при больших пиках и обеспечивают оценку индексов точки изменения. В этом примере выбирается пороговое значение 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]. Сходимость к статистически оптимальным индексам точки изменения может быть облегчена окружением каждого индекса набором смежных индексов равной длины. Таким образом, каждое ребро связано с левым и правым распределением. Кластеризация ГМ изучает параметры максимального правдоподобия двух распределений. Затем он оптимально корректирует индексы с учетом полученных параметров.
Можно задать длину для набора смежных позиций, распределенных вокруг индексов точек изменения. В этом примере выбирается длина 5. Можно также проверить каждую точку изменения путем печати ее кластеров GM. В этом примере мы построим графики GM-кластеров для данных Chromosome 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-тесты перестановок. Сегмент включает все точки данных от одной точки изменения до следующей точки изменения или конца хромосомы. В этом примере выполняется 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 функция. Отметим, что геномные положения в наборе данных клеточной линии Кориелла находятся в парах килобазов. Поэтому при добавлении идеограмм к сюжету нужно будет конвертировать данные цитобанда из 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 на основе массива.
Redon, R., et al., «Глобальная вариация числа копий в геноме человека», Nature, 444 (7118): 444-54, 2006.
Pinkel, D., et al., «Анализ вариаций числа копий ДНК с высоким разрешением с использованием сравнительной геномной гибридизации с микрочипами», Nature Genetics, 20 (2): 207-11, 1998.
Snijders, A.M., et al., «Сборка микрочипов для измерения числа копий ДНК по всему геному», Nature Genetics, 29 (3): 263-4, 2001.
[4] Геном человека NCBI Build 36.
[5] Myers, C.L., et al., «Точное обнаружение анеуплоидий в массиве CGH и данных микрочипов экспрессии генов», Bioinformatics, 20 (18): 3533-43, 2004.