В этом примере показано, как использовать MATLAB ® и Bioinformatics Toolbox™ для предварительной обработки данных на уровне зондов олигонуклеотидных микроматриц Affymrix ® с помощью двух методов предварительной обработки, RMA и GC Frustive Multi-Array Average (GCRMA).
На платформах олигонуклеотидных микрочипов Affymetrix экспрессию генов измеряют с использованием наборов зондов, состоящих из 11-20 идеальных зондов совпадения (PM) (длиной 25 нуклеотидов), комплементарных последовательностям мРНК-мишени. Каждый набор зондов также имеет одинаковое количество зондов рассогласования (MM), в которых 13-й нуклеотид был изменен на его комплемент. Зонды PM предназначены для генспецифической гибридизации. Считается, что измерения контрольных ММ зондов включают большую часть фонового неспецифического связывания, такого как перекрестная гибридизация. Зонд PM и соответствующий ему зонд MM называются парой зондов.
Измеренные интенсивности зонда и местоположения из гибридного микрочипа хранятся в CEL-файле. Для каждой платформы микрочипов Affymetrix информация, относящаяся к парам зондов для идентификаторов наборов зондов и расположениям на массиве, хранится в файле библиотеки CDF. Информация о последовательности зондов хранится в файле последовательности (в формате FASTA или с разделением табуляцией).
В целом, предварительная обработка данных экспрессии на уровне зонда Affymetrix состоит из трех этапов: корректировка фона, нормализация и суммирование на уровне набора зондов как мера уровня экспрессии соответствующей мРНК. Для статистических процедур этих трех этапов существует множество методов. В этом примере используются два популярных метода: RMA (Irizarry et al., 2003) и GCRMA (Wu et al., 2004).
Примечание.В этом примере показаны процедуры предварительной обработки RMA и GCRMA для вычисления значений выражений из входных файлов CEL с пошаговой детализацией с использованием нескольких функций. Можно также выполнить те же методы RMA или GCRMA в одном вызове функции с помощью панели инструментов биоинформатики. affyrma или affygcrma соответственно.
Для этого примера используют общедоступный набор данных, содержащий измерения микрочипов Affymetrix 42 опухолевых тканей эмбриональной центральной нервной системы (CNS, Pomeroy et al., 2002). Вы импортируете и получаете доступ к данным уровня зонда нескольких массивов, а затем выполняете измерения уровня выражения с помощью методов предварительной обработки RMA и GCRMA.
Эксперимент CNS проводился с использованием массива Affymetrix HuGeneFL GeneChip ®, и данные хранились в файлах CEL. Информация, относящаяся к каждому зонду, содержится в файле библиотеки Affymetrix Hu6800 CDF.
Если у вас уже нет файла библиотеки Hu6800 CDF, загрузите zip-файл библиотеки HuGeneFL Genome Array. Извлечение файла Hu6800.CDF в каталог, например C:\Examples\affypreprocessdemo\libfiles. Примечание.Чтобы получить доступ к файлам библиотеки, необходимо зарегистрироваться, но не запускать файл setup.exe в архиве.
Набор данных CNS (файлы CEL) доступен здесь. Для завершения этого примера загрузите файлы CEL набора данных CNS в каталог, например: C:\Examples\affypreprocessdemo\data. Распакуйте архивы файлов CEL. Примечание.Этот набор данных содержит больше CEL-файлов, чем необходимо для данного примера.
CNS_DataA_Sample_CEL.txt, файл, поставляемый с Bioinformatics Toolbox, содержит список 42 имен файлов CEL, используемых для этого примера, и образцов (10 медуллобластом, 10 рабдоидов, 10 злокачественных глиом, 8 супратенториальных PNETS и 4 нормальных мозжечка человека), к которым они принадлежат. Загрузите эти данные в две переменные MATLAB.
fid = fopen('CNS_DataA_Sample_CEL.txt','r'); ftext = textscan(fid,'%q%q'); fclose(fid); samples = ftext{1}; cels = ftext{2};
Установка переменных celPath и libPath в пути к файлам CEL и каталогам библиотек.
celPath = 'C:\Examples\affypreprocessdemo\data'; libPath = 'C:\Examples\affypreprocessdemo\libfiles';
Переименуйте файлы cel так, чтобы каждое имя файла начиналось с номера MG, следующего за подчеркиванием «_» в исходном имени файла. Например, GSM1688666_MG1999060202AA.CEL переименовано в MG1999060202AA.CEL. Не нужно запускать этот код, если имена файлов уже имеют требуемый формат.
A = dir(fullfile(celPath,'*.cel')); fileNames = string({A.name}); for iFile = 1:numel(A) newName = fullfile(celPath,extractAfter(fileNames(iFile),"_")); movefile(fullfile(celPath,fileNames(iFile)),newName); end
Функция celintensityread может считывать несколько файлов CEL и получать доступ к файлу библиотеки CDF. Он возвращает структуру MATLAB, содержащую информацию зонда и интенсивности зонда. Матрицы интенсивностей PM и MM из нескольких файлов CEL хранятся в PMIntensities и MMIntensities поля. В каждой матрице интенсивности зонда индексы столбцов соответствуют порядку, в котором считывались файлы CEL, и каждая строка соответствует зонду. Создайте структуру MATLAB интенсивностей зондов PM и MM, загрузив данные из файлов CEL из каталога, где хранятся файлы CEL, и передайте путь к файлу библиотеки CDF. Примечание: celintensityread сообщает о ходе выполнения в окно команды MATLAB. Можно отключить отчет о ходе выполнения, задав входной параметр VERBOSE в значение false.)
probeData = celintensityread(cels, 'Hu6800.CDF',... 'celpath', celPath, 'cdfpath', libPath, 'pmonly', false)
Reading CDF file: Hu6800.CDF
Reading file 1 of 42: MG2000040501AA
Reading file 2 of 42: MG2000040502AA
Reading file 3 of 42: MG2000040504AA
Reading file 4 of 42: MG2000040505AA
Reading file 5 of 42: MG2000040508AA
Reading file 6 of 42: MG2000040509AA
Reading file 7 of 42: MG2000040510AA
Reading file 8 of 42: MG2000040511AA
Reading file 9 of 42: MG2000040512AA
Reading file 10 of 42: MG2000040513AA
Reading file 11 of 42: MG2000051201AA
Reading file 12 of 42: MG2000051202AA
Reading file 13 of 42: MG2000051204AA
Reading file 14 of 42: MG2000051205AA
Reading file 15 of 42: MG2000051209AA
Reading file 16 of 42: MG2000071102AA
Reading file 17 of 42: MG2000051207AA
Reading file 18 of 42: MG2000051208AA
Reading file 19 of 42: MG2000051211AA
Reading file 20 of 42: MG2000051213AA
Reading file 21 of 42: MG2000061902AA
Reading file 22 of 42: MG2000061903AA
Reading file 23 of 42: MG2000061904AA
Reading file 24 of 42: MG2000061905AA
Reading file 25 of 42: MG2000061906AA
Reading file 26 of 42: MG2000070709AA
Reading file 27 of 42: MG2000070710AA
Reading file 28 of 42: MG2000070711AA
Reading file 29 of 42: MG2000070712AA
Reading file 30 of 42: MG2000070713AA
Reading file 31 of 42: MG1999112206AA
Reading file 32 of 42: MG2000033109AA
Reading file 33 of 42: MG2000033106AA
Reading file 34 of 42: MG2000033107AA
Reading file 35 of 42: MG1999112202AA
Reading file 36 of 42: MG1999112204AA
Reading file 37 of 42: MG2000011801AA
Reading file 38 of 42: MG2000031503AA
Reading file 39 of 42: MG2000032015AA
Reading file 40 of 42: MG2000030308AA
Reading file 41 of 42: MG2000011803AA
Reading file 42 of 42: MG2000011807AA
probeData =
struct with fields:
CDFName: 'Hu6800.CDF'
CELNames: {1×42 cell}
NumChips: 42
NumProbeSets: 7129
NumProbes: 140983
ProbeSetIDs: {7129×1 cell}
ProbeIndices: [140983×1 uint8]
GroupNumbers: [140983×1 uint8]
PMIntensities: [140983×42 single]
MMIntensities: [140983×42 single]
Определите количество загруженных файлов CEL.
nSamples = probeData.NumChips
nSamples =
42
Определите количество наборов зондов в массиве HuGeneFL.
nProbeSets = probeData.NumProbeSets
nProbeSets =
7129
Определите количество зондов в массиве HuGeneFL.
nProbes = probeData.NumProbes
nProbes =
140983
Для выполнения предварительной обработки GCRMA также требуется информация о последовательности зондов массива HuGeneFL. Сайт поддержки Affymetrix предоставляет информацию о последовательности тестов для большинства доступных массивов в виде файлов с форматированием FASTA или с разделением табуляцией. В этом примере предполагается, что файл HuGeneFL_probe_tab находится в каталоге файлов библиотеки. Используйте функцию affyprobeseqread для синтаксического анализа файла последовательности и возврата последовательностей зондов в nProbes x 25-я матрица целых чисел, которая представляет базы последовательностей PM-зондов, со строками, соответствующими зондам на кристалле, и столбцами, соответствующими базовым позициям 25-мерных.
S = affyprobeseqread('HuGeneFL_probe_tab', 'Hu6800.CDF',... 'seqpath', libPath, 'cdfpath', libPath, 'seqonly', true)
S =
struct with fields:
SequenceMatrix: [140983×25 uint8]
Процедура RMA использует только интенсивности зондов PM для корректировки фона (Irizarry et al., 2003), в то время как GCRMA корректирует фон с использованием информации о последовательности зондов и интенсивности контрольных зондов MM для оценки неспецифического связывания (Wu et al., 2004). И RMA, и GCRMA предшествуют квантовая нормализация (Bolstad et al., 2003) и медианная польская суммирование (Irizarry et al., 2003) интенсивностей ТЧ.
Метод настройки фона RMA корректирует интенсивность PM зонда по кристаллу. Интенсивности зонда PM моделируются как сумма составляющей нормального шума и составляющей экспоненциального сигнала. Использовать rmabackadj для фоновой корректировки интенсивности PM в данных CNS. Вы можете проверить гистограмму распределения интенсивности и предполагаемую корректировку фона конкретной микросхемы, установив входной параметр SHOWPLOT к индексу столбца микросхемы.
pms_bg = rmabackadj(probeData.PMIntensities, 'showplot', 1);
Несколько методов нелинейной нормализации были успешно применены к данным микрочипов Affymetrix. Процедура RMA нормализует данные уровня зонда методом квантовой нормализации. Использовать quantilenorm для нормализации скорректированной в фоновом режиме интенсивности ТЧ в данных ЦНС. Примечание.Если вас интересует метод нормализации рангового инвариантного набора, используйте affyinvarsetnorm вместо этого функция.
pms_bgnorm = quantilenorm(pms_bg);
К интенсивностям ТЧ при суммировании применяется процедура медианного полирования. Чтобы вычислить значения выражения, используйте rmasummary суммировать интенсивности зондов каждого набора зондов на множестве чипов. Значения выражений представляют собой сводки интенсивности набора зондов по шкале log-2.
cns_rma_exp = rmasummary(probeData.ProbeIndices, pms_bgnorm);
Процедура GCRMA корректируется для оптического шума и неспецифического связывания (NSB) с учетом эффекта более сильного связывания пар G/C (Naef et al., 2003, Wu et al., 2004). GCRMA использует информацию о последовательности зондов для оценки сродства зондов для вычисления неспецифического связывания. Сродство зонда моделируется как сумма зависящих от положения базовых эффектов. Обычно аффинность зонда оценивают по интенсивностям ММ эксперимента NSB. Если данные NSB недоступны, аффинности зонда все еще могут быть оценены на основе информации о последовательности и интенсивностей зонда MM, нормализованных установленной медианной интенсивностью зонда (Naef et al., 2003).
Для набора данных ЦНС используйте данные из микрочипа, гибридизированного с нормальным образцом мозжечка (Brain_Ncer_1), чтобы вычислить аффинности зонда для массива HuGeneFL. Использовать affyprobeaffinities для оценки сродства зонда к микрочипу Affymetrix. Используйте SHOWPLOT входной параметр для проверки графика, показывающего влияние оснований A, C, G и T в 25 положениях.
figure idx = find(strcmpi('Brain_Ncer_1', samples)); [pmAlpha, mmAlpha] = affyprobeaffinities(S.SequenceMatrix,... probeData.MMIntensities(:, idx), 'showplot', true);
Примечание: в массиве HuGeneFL 496 зондов не имеют информации о последовательности; сродством к этим зондам был NaN.
При доступных сродствах зондов количество NSB может быть оценено путем подгонки кривой LOWESS через интенсивности зондов ММ против сродств зондов ММ. Функция gcrmabackadj выполняет коррекцию оптических и NSB. Входной параметр SHOWPLOT На показан график интенсивности ММ с регулировкой оптического шума по отношению к его сродству и гладкой посадке указанного чипа. Можно вычислить фоновые интенсивности с помощью одного из двух методов оценки, оценки максимального правдоподобия (MLE) и эмпирической оценки Байеса (EB), которая вычисляет заднее среднее специфического связывания, заданного ранее наблюдаемыми интенсивностями. Здесь выполняется фоновая корректировка четырех массивов с использованием обоих методов оценки. Примечание: gcrmabackadj сообщает о ходе выполнения в окно команды MATLAB. Можно отключить отчет о ходе выполнения, задав входной параметр VERBOSE в значение false.)
В фоновом режиме отрегулируйте первые четыре микросхемы с использованием метода GCRMA-MLE и проверьте график зависимости интенсивности от сродства для данных из третьего массива.
pms_MLE_bg = gcrmabackadj(probeData.PMIntensities(:,1:4),... probeData.MMIntensities(:, 1:4),... pmAlpha, mmAlpha, 'showplot', 3); % Adjust YLIM for better view ylim([4 16]);
Adjusting background for chip # 1 of 4 using MLE method. Adjusting background for chip # 2 of 4 using MLE method. Adjusting background for chip # 3 of 4 using MLE method. Adjusting background for chip # 4 of 4 using MLE method.
Настройка первых четырех микросхем в фоновом режиме с использованием метода GCRMA-EB. Обработка с помощью этого метода является более вычислительной интенсивностью и займет больше времени.
pms_EB_bg = gcrmabackadj(probeData.PMIntensities(:,1:4),... probeData.MMIntensities(:, 1:4),... pmAlpha, mmAlpha, 'method', 'EB');
Adjusting background for chip # 1 of 4 using EB method. Adjusting background for chip # 2 of 4 using EB method. Adjusting background for chip # 3 of 4 using EB method. Adjusting background for chip # 4 of 4 using EB method.
Можно продолжить предварительную обработку с помощью quatilenorm и rmasummary функции или используйте gcrma функция, чтобы сделать все. gcrma функция выполняет корректировку фона и возвращает измерения выражений скорректированных по фону интенсивностей зонда PM с использованием тех же методов нормализации и суммирования, что и RMA. Можно также передать матрицу последовательности вместо сродств. В этом случае функция автоматически вычисляет аффинности. Примечание: gcrma сообщает о ходе выполнения в окно команды MATLAB. Можно отключить отчет о ходе выполнения, задав входной параметр VERBOSE в значение false.)
cns_mle_exp = gcrma(probeData.PMIntensities, probeData.MMIntensities,...
probeData.ProbeIndices, pmAlpha, mmAlpha);
Adjusting background for chip # 1 of 42 using MLE method. Adjusting background for chip # 2 of 42 using MLE method. Adjusting background for chip # 3 of 42 using MLE method. Adjusting background for chip # 4 of 42 using MLE method. Adjusting background for chip # 5 of 42 using MLE method. Adjusting background for chip # 6 of 42 using MLE method. Adjusting background for chip # 7 of 42 using MLE method. Adjusting background for chip # 8 of 42 using MLE method. Adjusting background for chip # 9 of 42 using MLE method. Adjusting background for chip # 10 of 42 using MLE method. Adjusting background for chip # 11 of 42 using MLE method. Adjusting background for chip # 12 of 42 using MLE method. Adjusting background for chip # 13 of 42 using MLE method. Adjusting background for chip # 14 of 42 using MLE method. Adjusting background for chip # 15 of 42 using MLE method. Adjusting background for chip # 16 of 42 using MLE method. Adjusting background for chip # 17 of 42 using MLE method. Adjusting background for chip # 18 of 42 using MLE method. Adjusting background for chip # 19 of 42 using MLE method. Adjusting background for chip # 20 of 42 using MLE method. Adjusting background for chip # 21 of 42 using MLE method. Adjusting background for chip # 22 of 42 using MLE method. Adjusting background for chip # 23 of 42 using MLE method. Adjusting background for chip # 24 of 42 using MLE method. Adjusting background for chip # 25 of 42 using MLE method. Adjusting background for chip # 26 of 42 using MLE method. Adjusting background for chip # 27 of 42 using MLE method. Adjusting background for chip # 28 of 42 using MLE method. Adjusting background for chip # 29 of 42 using MLE method. Adjusting background for chip # 30 of 42 using MLE method. Adjusting background for chip # 31 of 42 using MLE method. Adjusting background for chip # 32 of 42 using MLE method. Adjusting background for chip # 33 of 42 using MLE method. Adjusting background for chip # 34 of 42 using MLE method. Adjusting background for chip # 35 of 42 using MLE method. Adjusting background for chip # 36 of 42 using MLE method. Adjusting background for chip # 37 of 42 using MLE method. Adjusting background for chip # 38 of 42 using MLE method. Adjusting background for chip # 39 of 42 using MLE method. Adjusting background for chip # 40 of 42 using MLE method. Adjusting background for chip # 41 of 42 using MLE method. Adjusting background for chip # 42 of 42 using MLE method. Normalizing. Calculating expression.
Для проверки распределения интенсивности PM первых четырех чипов с тремя процедурами фоновой регулировки используйте боксплоты.
figure subplot(4,1,1) maboxplot(log2(probeData.PMIntensities(:, 1:4)), samples(1:4),... 'title','Raw Intensities', 'orientation', 'horizontal') subplot(4,1,2) maboxplot(log2(pms_bg(:,1:4)), samples(1:4),... 'title','After RMA Background Adjustment','orient','horizontal') subplot(4,1,3) maboxplot(log2(pms_MLE_bg), samples(1:4),... 'title','After GCRMA-MLE Background Adjustment','orient','horizontal') subplot(4,1,4) maboxplot(log2(pms_EB_bg), samples(1:4),... 'title','After GCRMA-EB Background Adjustment','orient','horizontal')
Используйте сводные графики для проверки скорректированных и нормированных распределений интенсивности ТЧ первых четырех микросхем с помощью трех процедур корректировки фона.
pms_MLE_bgnorm = quantilenorm(pms_MLE_bg); pms_EB_bgnorm = quantilenorm(pms_EB_bg); figure subplot(3,1,1) maboxplot(log2(pms_bgnorm(:, 1:4)), samples(1:4),... 'title','Normalized RMA Background Adjusted Intensity',... 'orientation', 'horizontal') subplot(3,1,2) maboxplot(log2(pms_MLE_bgnorm), samples(1:4),... 'title','Normalized GCRMA-MLE Background Adjusted Intensity',... 'orientation', 'horizontal') subplot(3,1,3) maboxplot(log2(pms_EB_bgnorm), samples(1:4),... 'title','Normalized GCRMA-EB Background Adjusted Intensity',... 'orientation', 'horizontal')
Импорт данных из файлов CEL и всех трех этапов предварительной обработки методов RMA и GCRMA, показанных в этом примере, можно выполнить с помощью команды affyrma и affygcrma соответственно.
Дополнительные сведения об анализе экспрессии генов с помощью Bioinformatics Toolbox см. в разделе Исследование данных экспрессии генов микрочипов.
Affymetrix и GeneChip являются зарегистрированными товарными знаками Affymetrix, Inc.
[1] Pomeroy, S.L., et al., «Предсказание исхода эмбриональной опухоли центральной нервной системы на основе экспрессии генов», Nature, 415 (6870): 436-42, 2002.
[2] Irizarry, R.A., et al., «Исследование, нормализация и резюме данных уровня зонда массива олигонуклеотидов высокой плотности», Biostatistics, 4 (2): 249-64, 2003.
[3] Wu, Z., et al., «Основанная на модели корректировка фона для массивов экспрессии олигонуклеотидов», Journal of the American Statistical Association, 99 (468): 909-17, 2004.
[4] Bolstad, B.M., et al., «Сравнение методов нормализации данных массива олигонуклеотидов высокой плотности, основанных на дисперсии и смещении», Bioinformatics, 19 (2): 185-93, 2003.
[5] Наеф, Ф. и Магнаско, М.О. «Решение загадки ярких несовпадений: мечение и эффективное связывание в олигонуклеотидных массивах», Physical Review, E, Statistical, Nonlinear and Soft Matter Physics, 68 (1Pt1): 011906, 2003.