Предварительная обработка данных микрочипов Affymetrix ® на уровне зонда
В этом примере показано, как использовать Toolbox™ MATLAB ® и Bioinformatics для предварительной обработки данных уровня зонда олигонуклеотидных микрочипов Affymetrix ® с двумя методами предварительной обработки, Robust Multi-array Average (RMA) и GC Robust Multi-array avid (Gray avid)
Введение
С платформами для микромассивов олигонуклеотидов Affymetrix экспрессию генов измеряют с использованием наборов зондов, состоящих из 11-20 точно соответствующих (PM) зондов (длина 25 нуклеотидов), комплементарных последовательностям мРНК-мишеней. Каждый набор зондов также имеет одинаковое количество зондов несовпадения (ММ), в которых 13-й нуклеотид был изменен на его комплемент. PM-зонды предназначены для ген-специфической гибридизации. Предполагается, что измерения контрольного MM-зонда включают большую часть фонового неспецифического связывания, такого как перекрестная гибридизация. PM-зонд и соответствующий ему MM-зонд называются парой зондов.
Измеренная интенсивность зонда и местоположения из гибридизированной микромассива сохраняются в файле CEL. Для каждой платформы микромассивов Affymetrix информация, связывающая пары зондов с идентификаторами наборов зондов и с местоположениями в массиве, хранится в файле библиотеки CDF. Информация о последовательности зондов хранится в файле последовательности (FASTA или формате, разделенном табуляцией).
В целом предварительная обработка данных экспрессии на уровне зонда Affymetrix состоит из трех шагов: корректировка фона, нормализация и суммирование на уровне набора зондов как мера уровня экспрессии соответствующей мРНК. Существует много методов для статистических процедур этих трех этапов. В этом примере используют два популярных метода, RMA (Irizarry et al., 2003) и GCRMA (Wu et al., 2004).
Примечание. В этом примере показаны процедуры предварительной обработки RMA и GCRMA для вычисления значений выражений из входных файлов CEL в пошаговой детализации с использованием нескольких функций. Можно также выполнить те же методы RMA или GCRMA в одном вызове функции с помощью Bioinformatics Toolbox affyrma или affygcrma функций, соответственно.
В этом примере используется общедоступный набор данных, содержащий измерения микромассивов 42 опухолевых тканей эмбриональной центральной нервной системы (CNS, Pomeroy et al., 2002). Вы будете импортировать и получить доступ к данным уровня зонда из нескольких массивов, а затем выполнить измерения уровня выражения с помощью методов предварительной обработки RMA и GCRMA.
Импорт данных
Эксперимент ЦНС проводили с использованием массива 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 к ложному.)
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
Способ фоновой регулировки RMA корректирует элементарную величину интенсивности PM-зонда за чипом. Интенсивность ПМ зонда моделируется как сумма нормального шумового компонента и экспоненциального сигнального компонента. Использование rmabackadj в фоновом режиме корректируйте интенсивность ТЧ в данных CNS. Можно просмотреть гистограмму распределения интенсивности и предполагаемую настройку фона определенного чипа, установив параметр входа SHOWPLOT к индексу столбца чипа.
pms_bg = rmabackadj(probeData.PMIntensities, 'showplot', 1);
Несколько нелинейных методов нормализации были успешно применены к данным микромассивов Affymetrix. Процедура RMA нормализует данные уровня зонда с помощью метода квантильной нормализации. Использование quantilenorm для нормализации фоновой скорректированной интенсивности ТЧ в данных CNS. Примечание: Если вас интересует метод нормализации инвариантного набора рангов, используйте affyinvarsetnorm вместо этого функция.
pms_bgnorm = quantilenorm(pms_bg);
Медианная процедура полирования применяется к интенсивности ТЧ при суммировании. Чтобы вычислить значения выражения, используйте rmasummary суммирование интенсивности зондирования каждого набора зондов на нескольких микросхемах. Значения выражения являются сводными данными интенсивности набора зондов по шкале журнала -2.
cns_rma_exp = rmasummary(probeData.ProbeIndices, pms_bgnorm);
Использование процедуры GCRMA
Процедура GCRMA регулирует оптический шум и неспецифическое связывание (NSB) с учетом эффекта более сильного связывания G/C пар (Naef et al., 2003, Wu et al., 2004). GCRMA использует информацию о последовательности зондов для оценки сродства зондов для вычисления неспецифического связывания. Сродство зонда моделируется как сумма позиционно-зависимых базовых эффектов. Обычно сродство зонда оценивают из интенсивности ММ эксперимента NSB. Если данные NSB недоступны, сродство зонда все еще может быть оценено из информации о последовательности и интенсивности зонда MM, нормированной средней интенсивностью зонда (Naef et al., 2003).
Для набора данных CNS используйте данные из микромассива, гибридизированного с нормальной выборкой церебеллы (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);
Примечание: Существует 496 зондов на массиве HuGeneFL, которые не имеют информации о последовательности; сродство к этим зондам было NaN.
При наличии сродства зонда количество NSB может быть оценено путем подгонки кривой LOWESS через интенсивность зонда MM по сравнению с сродством зонда MM. Функция gcrmabackadj выполняет оптические и NSB коррекции. Параметр входа SHOWPLOT показывает график интенсивности MM, скорректированной по оптическому шуму, относительно его сродства и плавности подгонки заданного чипа. Можно вычислить интенсивность фона с помощью одного из двух методов оценки, Maximum Likelihood Estimate (MLE) и Empirical-Bayes (EB), который вычисляет апостериорное среднее значения специфического связывания с учетом предшествующих наблюдаемых интенсивностей. Здесь вы будете корректировать фон четырех массивов, используя оба метода оценки. (Примечание: gcrmabackadj сообщит о прогрессе в командное окно MATLAB. Можно выключить отчет о прогрессе, установив параметр входа VERBOSE к ложному.)
Фоновая настройка первых четырех микросхем с помощью метода 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 функция, чтобы сделать все. The gcrma функция выполняет настройку фона и возвращает показатели выражения интенсивности скорректированного фона PM-зонда с помощью тех же методов нормализации и суммирования, что и RMA. Можно также передать в матрице последовательности вместо аффинностей. Функция автоматически вычисляет сродства в этом случае. (Примечание: gcrma сообщит о прогрессе в командное окно MATLAB. Можно выключить отчет о прогрессе, установив параметр входа VERBOSE к ложному.)
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.
Просмотр результатов корректировки фона
Используйте boxplots, чтобы проверить распределения интенсивности ТЧ первых четырех микросхем с тремя процедурами фоновой настройки.
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')
Используйте boxplots, чтобы просмотреть скорректированный фон и нормализованные распределения интенсивности ТЧ первых четырех микросхем с тремя процедурами настройки фона.
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., «Exploration, normalization, and сводных данных of high densition oligonucleotide массива sense level данных», 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] Naef, F. and Magnasco, M.O. «Решение загадки ярких несоответствий: маркировка и эффективное связывание в олигонуклеотидных массивах», Physical Review, E, Statistical, Nonlinear and Soft Matter Physics, 68 (1Pt1): 011906, 2003.