Введение

В типичном эксперименте с дефенированной масс-спектрометрией белки и метаболиты собирают из клеток, тканей или жидкостей организма, растворяют и денатурируют в растворе и ферментативно расщепляют в смеси. Эти смеси затем разделяют либо высокоэффективной жидкостной хроматографией (ВЭЖХ), либо капиллярным электрофорезом (СЕ), либо газовой хроматографией (ГХ) и соединяют с методом идентификации методом масс-спектрометрии, таким как электроспрей-ионизационная масс-спектрометрия (ESI-MS), матричная ионизация (MALDI или

Открытые репозитории данных и формат файлов mzXML

Для этого примера мы используем тестовый образец данных LC-ESI-MS с смесью семи белков. Данные в этом примере (7MIX_STD_110802_1) взяты из репозитария данных Sashimi. Набор данных не распределяется с MATLAB ®. Для завершения этого примера необходимо загрузить набор данных в локальный каталог или собственный репозиторий. В качестве альтернативы вы можете попробовать другие наборы данных, доступные в других общедоступных базах данных для данных экспрессии белка, таких как Пептидный атлас в Институте системной биологии [1].

Большинство текущих масс-спектрометров могут преобразовывать или сохранять данные сбора данных с помощью схемы mzXML. Этот стандарт представляет собой общий формат файлов на основе XML (eXtensible Markup Language), разработанный в рамках проекта Sashimi для решения проблем, связанных с представлением наборов данных от различных производителей и из различных экспериментальных установок в общую и расширяемую схему. Файлы mzXML, используемые в масс-спектрометрии с переносами, обычно очень велики. MZXMLINFO функция позволяет получать основную информацию о наборе данных, не считывая его в память. Например, можно получить количество сканирований, диапазон времени хранения и количество тандемных приборов MS (уровней).

info = mzxmlinfo('7MIX_STD_110802_1.mzXML','NUMOFLEVELS',true)

info = 

  struct with fields:

                            Filename: '7MIX_STD_110802_1.mzXML'
                         FileModDate: '01-Feb-2013 11:54:30'
                            FileSize: 26789612
                       NumberOfScans: 7161
                           StartTime: 'PT0.00683333S'
                             EndTime: 'PT200.036S'
       DataProcessingIntensityCutoff: 'N/A'
            DataProcessingCentroided: 'true'
            DataProcessingDeisotoped: 'N/A'
    DataProcessingChargeDeconvoluted: 'N/A'
       DataProcessingSpotIntegration: 'N/A'
                    NumberOfMSLevels: 2

MZXMLREAD функция считывает XML-документ в структуру MATLAB. Области scan и index размещают на первом уровне выходной структуры для улучшения доступа к спектральным данным. Оставшаяся часть дерева документов mzXML анализируется в соответствии со спецификациями схемы. Этот набор данных LC/MS содержит 7161 сканирование с двумя уровнями MS. В этом примере будут использоваться только сканирования первого уровня. Спектры второго уровня обычно используются для идентификации пептида/белка и приходят на более поздней стадии в некоторых типах анализа рабочего процесса. MZXMLREAD можно фильтровать необходимые сканирования без загрузки всего набора данных в память:

mzXML_struct = mzxmlread('7MIX_STD_110802_1.mzXML','LEVEL',1)

mzXML_struct = 

  struct with fields:

     scan: [2387x1 struct]
    mzXML: [1x1 struct]
    index: [1x1 struct]

Если при загрузке возникают ошибки, связанные с памятью или пространством кучи Java, попробуйте увеличить пространство кучи Java, как описано здесь.

Более подробная информация, относящаяся к масс-спектрометру и экспериментальным условиям, находится в полевых условиях. msRun.

mzXML_struct.mzXML.msRun

ans = 

  struct with fields:

         scanCount: 7161
         startTime: "PT0.00683333S"
           endTime: "PT200.036S"
        parentFile: [1x1 struct]
      msInstrument: [1x1 struct]
    dataProcessing: [1x1 struct]

Для облегчения обработки данных, MZXML2PEAKS функция извлекает список пиков из каждого сканирования в массив ячеек (peaks]) и их соответствующее время удержания в векторе столбца (time). Можно извлечь спектры определенного уровня, установив LEVEL входной параметр.

[peaks,time] = mzxml2peaks(mzXML_struct);
numScans = numel(peaks)

numScans =

        2387

MSDOTPLOT создает обзор наиболее интенсивных пиков во всем наборе данных. В этом случае мы визуализируем только наиболее интенсивные 5% пики интенсивности ионов, задав входной параметр QUANTILE до 0,95.

h = msdotplot(peaks,time,'quantile',.95);
title('5 Percent Overall Most Intense Peaks')

Можно также фильтровать пики отдельно для каждого сканирования, используя процентиль интенсивности базового пика. Базовый пик является наиболее интенсивным пиком, обнаруженным в каждом сканировании [2]. Этот параметр задается автоматически большинством спектрометров. Эта операция требует запроса в структуру mxXML для получения базовой пиковой информации. Обратите внимание, что можно также найти интенсивность базового пика, выполнив итерацию MAX функция над списком пиков.

basePeakInt = [mzXML_struct.scan.basePeakIntensity]';
peaks_fil = cell(numScans,1);
for i = 1:numScans
    h = peaks{i}(:,2) > (basePeakInt(i).*0.75);
    peaks_fil{i} = peaks{i}(h,:);
end

whos('basePeakInt','level_1','peaks','peaks_fil')
msdotplot(peaks_fil,time)
title('Peaks Above (0.75 x Base Peak Intensity) for Each Scan')

  Name                Size               Bytes  Class     Attributes

  basePeakInt      2387x1                19096  double              
  peaks            2387x1             14031800  cell                
  peaks_fil        2387x1               289568  cell                

Можно настроить 3-D обзор отфильтрованных пиков с помощью STEM3 функция. STEM3 функция требует, чтобы поместить данные в три вектора, элементы которых образуют триплеты (время удержания, масса/заряд и значение интенсивности), которые представляют каждый стебель.

peaks_3D = cell(numScans,1);
for i = 1:numScans
    peaks_3D{i}(:,[2 3]) = peaks_fil{i};
    peaks_3D{i}(:,1) = time(i);
end
peaks_3D = cell2mat(peaks_3D);

figure
stem3(peaks_3D(:,1),peaks_3D(:,2),peaks_3D(:,3),'marker','none')
axis([0 12000 400 1500 0 1e9])
view(60,60)
xlabel('Retention Time (seconds)')
ylabel('Mass/Charge (M/Z)')
zlabel('Relative Ion Intensity')
title('Peaks Above (0.75 x Base Peak Intensity) for Each Scan')

Цветные стебли можно распечатать с помощью PATCH функция. Для каждой тройни в peaks_3Dперемежать новый триплет со значением интенсивности, равным нулю. Затем создайте цветовой вектор в зависимости от интенсивности стебля. Логарифмическое преобразование увеличивает динамический диапазон карты цветов. Для перемежающихся триплетов назначьте NaN, так что PATCH функция не рисует линии, соединяющие смежные стебли.

peaks_patch = sortrows(repmat(peaks_3D,2,1));
peaks_patch(2:2:end,3) = 0;

col_vec = log(peaks_patch(:,3));
col_vec(2:2:end) = NaN;

figure
patch(peaks_patch(:,1),peaks_patch(:,2),peaks_patch(:,3),col_vec,...
      'edgeColor','flat','markeredgecolor','flat','Marker','x','MarkerSize',3);
axis([0 12000 400 1500 0 1e9])
view(90,90)
xlabel('Retention Time (seconds)')
ylabel('Mass/Charge (M/Z)')
zlabel('Relative Ion Intensity')
title('Peaks Above (0.75 x Base Peak Intensity) for Each Scan')

view(40,40)

Создание тепловых карт наборов данных LC/MS

Общие методы в отрасли работают с пиковой информацией (a.k.a. центроидальные данные) вместо необработанных сигналов. Это может сэкономить память, но некоторые важные детали не видны, особенно когда необходимо проверить образцы со сложными смесями. Для дальнейшего анализа этого набора данных можно создать общую сетку в измерении массы/заряда. Поскольку не все сканирования имеют достаточно информации для восстановления исходного сигнала, мы используем метод повторной выборки с сохранением пика. Путем выбора соответствующих параметров для MSPPRESAMPLE функция позволяет гарантировать, что разрешение спектров не потеряно и что максимальные значения пиков коррелируют с исходной пиковой информацией.

[MZ,Y] = msppresample(peaks,5000);
whos('MZ','Y')

  Name         Size                 Bytes  Class     Attributes

  MZ        5000x1                  20000  single              
  Y         5000x2387            47740000  single              

С этой матрицей интенсивностей ионов, Yможно создать цветную тепловую карту. MSHEATMAP функция автоматически настраивает панель цветов, используемую для отображения статистически значимых пиков с горячими цветами и шумных пиков с холодными цветами. Алгоритм основан на кластеризации значимых пиков и шумных пиков путём оценки смеси гауссов с подходом ожидания-максимизации. Кроме того, можно использовать MIDPOINT входной параметр для выбора произвольного порога, чтобы отделить шумные пики от значительных пиков, или можно в интерактивном режиме сдвинуть карту цветов, чтобы скрыть или отобразить пики. При работе с тепловыми картами обычно отображается логарифм интенсивностей ионов, что увеличивает динамический диапазон карты цветов.

fh1 = msheatmap(MZ,time,log(Y),'resolution',.1,'range',[500 900]);
title('Original LC/MS Data Set')

Можно зумировать до небольших областей, представляющих интерес, чтобы наблюдать данные, в интерактивном или программном режиме с помощью AXIS функция. Мы наблюдаем некоторые регионы с высокой относительной интенсивностью ионов. Они представляют собой пептиды в биологическом образце.

axis([527 539 385 496])

Можно наложить исходную пиковую информацию набора данных LC/MS. Это позволяет оценить производительность операции повторной выборки с сохранением пиковых значений. Для скрытия/отображения точек можно использовать возвращенный маркер.

dp1 = msdotplot(peaks,time);

Калибровка местоположения масс/зарядов пиков по общей сетке

Двумерные пики кажутся шумными или они не имеют компактной формы в смежных спектрах. Это обычная проблема для многих масс-спектрометров. Случайные колебания значения массы/заряда, извлеченного из пиков повторных профилей, находятся в диапазоне от 0,1% до 0,3% [3]. Такая изменчивость может быть вызвана несколькими факторами, например плохой калибровкой детектора, низким отношением сигнал/шум или проблемами в алгоритмах выделения пиков. MSPALIGN функция реализует расширенные алгоритмы объединения данных, которые синхронизируют все спектры в наборе данных с общей сеткой масса/заряд (CMZ). CMZ может быть выбран произвольно или может быть оценен после анализа данных [2,4,5]. Процедура согласования пиков может использовать либо правило ближайшего соседа, либо динамическое программирование.

[CMZ, peaks_CMZ] = mspalign(peaks);

Повторите процесс визуализации с выровненными пиками: выполните повторную выборку с сохранением пиков, создайте тепловую карту, наложите информацию о выровненных пиках и увеличьте масштаб в той же области, что и раньше. При повторной калибровке спектра можно различать изотопные паттерны некоторых пептидов.

[MZ_A,Y_A] = msppresample(peaks_CMZ,5000);
fh2 = msheatmap(MZ_A,time,log(Y_A),'resolution',.10,'range',[500 900]);
title('LC/MS Data Set with the Mass/Charge Calibrated to a CMZ')
dp2 = msdotplot(peaks_CMZ,time);
axis([527 539 385 496])

Локальная калибровка местоположения масс/зарядов пиков

MSPALIGN вычисляет один CMZ для всего набора данных LC/MS. Это может быть не идеальным случаем для образцов с более сложными смесями пептидов и/или метаболитов, чем набор данных, использованный в этом примере. В случае сложных смесей можно выровнять каждый спектр по локальному набору спектров, которые содержат только информативные пики (высокая интенсивность) с аналогичными временами удерживания, в противном случае процесс калибровки в областях с небольшими пиками (низкая интенсивность) может быть смещен другими пиками, которые имеют одинаковые значения массы/заряда, но находятся в разных временах удерживания. Для выполнения более точной калибровки можно использовать SAMPLEALIGN функция. Эта функция является обобщением алгоритмов ограниченного динамического искажения времени (CDTW), обычно используемых при обработке речи [6]. В следующем for цикл, мы сохраняем буфер с интенсивностями предыдущих выровненных спектров (LAI). Интенсивности ионов спектров масштабируются с помощью анонимной функции SF (внутри SAMPLEALIGN) для уменьшения расстояния между большими пиками. SAMPLEALIGN уменьшает общее расстояние между всеми согласованными точками и при необходимости вводит промежутки. В этом случае мы используем более тонкий вектор MZ (FMZ), чтобы максимально сохранить правильное значение массы/заряда пиков. Примечание: это может занять некоторое время, так как алгоритм CDTW выполняется 2387 раз.

SF = @(x) 1-exp(-x./5e7); % scaling function
DF = @(R,S) sqrt((SF(R(:,2))-SF(S(:,2))).^2 + (R(:,1)-S(:,1)).^2);

FMZ = (500:0.15:900)';  % setup a finer MZ vector
LAI = zeros(size(FMZ)); % init buffer for the last alignment intensities

peaks_FMZ = cell(numScans,1);
for i = 1:numScans
    % show progress
    if ~rem(i,250)
        fprintf(' %d...',i);
    end
    % align peaks in current scan to LAI
    [k,j] = samplealign([FMZ,LAI],double(peaks{i}),'band',1.5,'gap',[0,2],'dist',DF);
    % updating the LAI buffer
    LAI = LAI*.25;
    LAI(k) = LAI(k) + peaks{i}(j,2);
    % save the alignment
    peaks_FMZ{i} = [FMZ(k) peaks{i}(j,2)];
end

 250... 500... 750... 1000... 1250... 1500... 1750... 2000... 2250...

Повторите процесс визуализации и увеличьте изображение интересующей области.

[MZ_B,Y_B] = msppresample(peaks_FMZ,4000);
fh3 = msheatmap(MZ_B,time,log(Y_B),'resolution',.10,'range',[500 900]);
title('LC/MS Data Set with the Mass/Charge Calibrated Locally')
dp3 = msdotplot(peaks_FMZ,time);
axis([527 539 385 496])

В качестве последнего шага для улучшения изображения можно применить фильтр Гаусса в хроматографическом направлении для сглаживания всего набора данных.

Gpulse = exp(-.1*(-10:10).^2)./sum(exp(-.1*(-10:10).^2));
YF = convn(Y_B,Gpulse,'same');
fh4 = msheatmap(MZ_B,time,log(YF),'resolution',.10,'limits',[500 900]);
title('Final Enhanced LC/MS Data Set')
dp4 = msdotplot(peaks_FMZ,time);
axis([527 539 385 496])

Можно связать оси двух тепловых карт для интерактивного или программного сравнения областей между двумя наборами данных. В этом случае мы сравниваем исходную и окончательную расширенные матрицы LC/MS.

linkaxes(findobj([fh1 fh4],'Tag','MSHeatMap'))
axis([521 538 266 617])

Извлечение спектров с помощью общей ионной хроматограммы

Как только набор данных LC/MS сглаживается и повторно дискретизируется в обычную сетку, можно извлечь наиболее информативные спектры, рассматривая локальные максимумы общей ионной хроматограммы (TIC). TIC вычисляется путем суммирования строк YF. Затем используйте MSPEAKS поиск значений времени удержания для извлечения выбранных подмножеств спектров.

TIC = mean(YF);
pt = mspeaks(time,TIC','multiplier',10,'overseg',100,'showplot',true);
title('Extracted peaks from the Total Ion Chromatogram (TIC)')
pt(pt(:,1)>4000,:) = []; % remove spectra above 4000 seconds
numPeaks = size(pt,1)

numPeaks =

    22

Создайте 3-D график выбранных спектров.

xRows = samplealign(time,pt(:,1),'width',1); % finds the time index for every peak
xSpec = YF(:,xRows);                         % gets the signals to plot

figure;
hold on
box on
plot3(repmat(MZ_B,1,numPeaks),repmat(1:numPeaks,numel(MZ_B),1),xSpec)
view(20,85)

ax = gca;
ax.YTick = 1:numPeaks;
ax.YTickLabel = num2str(time(xRows));
axis([500 900 0 numPeaks 0 1e8])

xlabel('Mass/Charge (M/Z)')
ylabel('Time')
zlabel('Relative Ion Intensity')
title('Extracted Spectra Subset')

Маркеры наложения для извлеченных спектров по улучшенной тепловой карте.

linkaxes(findobj(fh4,'Tag','MSHeatMap'),'off')
figure(fh4)
hold on
for i = 1:numPeaks
    plot([400 1500],xRows([i i]),'m')
end
axis([500 900 100 925])
dp4.Visible = 'off';
title('Final Enhanced LC/MS Data Set with Extracted Spectra Marked')

Ссылки

[1] Desiere, F. et al., «Проект пептидного атласа», Nucleic Acids Research, 34: D655-8, 2006.

[2] Purvine, S., Kolker, N. и Kolker, E., «Спектральная оценка качества для высокопрочной тандемной масс-спектрометрии протеомики», OMICS: A Journal of Integrative Biology, 8 (3): 255-65, 2004.

[3] Kazmi, A.S., et al., «Согласование масс-спектров высокого разрешения: разработка эвристического подхода для метаболомики», Metabolomics, 2 (2): 75-83, 2006.

[4] Джеффрис, Н., «Алгоритмы выравнивания протеомных данных масс-спектрометрии», Биоинформатика, 21 (14): 3066-3073, 2005.

[5] Yu, W., et al., «Множественное выравнивание пиков при последовательном анализе данных: подход на основе масштабного пространства», IEEE ®/ACM Trans. Computational Biology and Bioinformatics, 3 (3): 208-219, 2006.

[6] Сакое, Х. и Тиба с., «Оптимизация алгоритма динамического программирования для распознавания разговорных слов», IEEE Trans. Acoustics, Speech and Signal Processing, ASSP-26 (1): 43-9, 1978.