Введение

Использование комплексных пептидных и метаболитных смесей в ЖХ/МС требует процедур выравнивания без меток. Анализ данных этого типа требует поиска статистически значимых различий между биологически связанными наборами данных без необходимости полной идентификации всех соединений в образце (либо пептидов/белков, либо метаболитов). Сравнение соединений требует выравнивания в двух измерениях: измерении массового заряда и измерении времени удержания [1]. В примерах Предварительная обработка необработанных данных масс-спектрометрии и визуализация и предварительная обработка наборов данных масс-спектрометрии для профилирования метаболитов и белков/пептидов можно узнать, как использовать MSALIGN, MSPALIGN, и SAMPLEALIGN функции для деформации или калибровки различных типов аномалий в размерности масса/заряд. В этом примере вы узнаете, как использовать SAMPLEALIGN функция также позволяет корректировать нелинейные и непредсказуемые изменения измерения времени удержания.

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

Описание набора данных

В этом примере используются два образца в PAe000153 и PAe000155, доступных из Пептидного атласа [2]. Образцы представляют собой LC-ESI-MS-сканирование четырех фракций солевого белка из saccharomyces cerevisiae, каждая из которых содержит более 1000 пептидов. Образцы дрожжей обрабатывали различными химическими веществами (глицином и серином), чтобы получить два биологически разнообразных образца. Согласование по времени этих двух наборов данных является одним из наиболее сложных случаев, описанных в [3]. Наборы данных не распределяются с MATLAB ®. Необходимо загрузить наборы данных в локальный каталог или собственный репозиторий. Кроме того, можно попробовать другие наборы данных, доступные в общедоступных базах данных для белковых данных, такие как Sashimi Data Repository. При получении ошибок, связанных с памятью или пространством кучи Java, попробуйте увеличить пространство кучи Java, как описано здесь. Анализ данных LC/MS требует больших объемов памяти от операционной системы; если вы получаете"Out of memory" ошибки при выполнении этого примера, попытка увеличения объема виртуальной памяти (или пространства подкачки) операционной системы или попытка установки переключателя 3GB (только для 32-разрядной Windows ® XP). Эти методы описаны в этом документе.

Прочтите и извлеките списки пиков из XML-файлов, содержащих данные интенсивности для образца, обработанного серином, и образца, обработанного глицином.

ser = mzxmlread('005_1.mzXML')
[ps,ts] = mzxml2peaks(ser,'level',1);
gly = mzxmlread('005a.mzXML')
[pg,tg] = mzxml2peaks(gly,'level',1);

ser = 

  struct with fields:

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


gly = 

  struct with fields:

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

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

[MZs,Ys] = msppresample(ps,5000);
[MZg,Yg] = msppresample(pg,5000);

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

fh1 = msheatmap(MZs,ts,log(Ys),'resolution',0.15);
title('Serine Treatment')
fh2 = msheatmap(MZg,tg,log(Yg),'resolution',0.15);
title('Glycine Treatment')

Детальный осмотр проблем с несоответствием

Обратите внимание, что наборы данных можно визуализировать отдельно; однако векторы времени имеют различный размер, и поэтому карты теплоты имеют различное количество строк (или Ys и Yg имеют различное количество столбцов). Кроме того, частота дискретизации не является постоянной, и сдвиг между векторами времени не является линейным.

whos('Ys','Yg','ts','tg')
figure
plot(1:numel(ts),ts,1:numel(tg),tg)
legend('Serine','Glycine','Location','NorthWest')
title('Time Vectors')
xlabel('Spectrum Index')
ylabel('Retention Time (seconds)')

  Name         Size                Bytes  Class     Attributes

  Yg        5000x921            18420000  single              
  Ys        5000x937            18740000  single              
  tg         921x1                  7368  double              
  ts         937x1                  7496  double              

Чтобы увидеть одну и ту же интересующую область в обоих наборах данных, необходимо вычислить соответствующие индексы строк в каждой матрице. Например, чтобы проверить пики пептида в диапазоне 480-520 Да MZ и диапазоне времени удержания 1550-1900 секунд, необходимо найти ближайшие совпадения для этого диапазона в векторах времени и затем увеличить каждый рисунок:

ind_ser = samplealign(ts,[1550;1900]);
figure(fh1);
axis([480 520 ind_ser'])

ind_gly = samplealign(tg,[1550;1900]);
figure(fh2);
axis([480 520 ind_gly'])

Несмотря на то, что вы увеличили масштаб в том же диапазоне, вы все еще можете наблюдать, что верхний правый пептид в осях смещен в измерении времени удерживания. В образце, обработанном серином, центр этого пика, по-видимому, находится примерно через 1595 секунд, в то время как в образце, обработанном глицином, предполагаемый тот же самый пептид находится примерно через 1630 секунд. Это предотвратит точный сравнительный анализ, даже если выполнить повторную выборку наборов данных в один и тот же временной вектор. В дополнение к сдвигу времени удерживания набор данных, по-видимому, неправильно откалиброван по размеру масса/заряд, поскольку пики не имеют компактной формы в смежных спектрах.

Калибровка массы/заряда и расширение матриц

Прежде чем корректировать время удерживания, можно улучшить образцы, используя подход, аналогичный описанному в примере Визуализация и предварительная обработка наборов данных масс-спектрометрии для профилирования метаболитов и белков/пептидов. Для краткости мы отображаем только код MATLAB без каких-либо дополнительных объяснений:

SF = @(x) 1-exp(-x./5e7); % scaling function
DF = @(R,S) sqrt((SF(R(:,2))-SF(S(:,2))).^2 + (R(:,1)-S(:,1)).^2);
CMZ = (315:.10:680)'; % Common Mass/Charge Vector with a finer grid

% Align peaks of the serine sample in the MZ direction
LAI = zeros(size(CMZ));
for i = 1:numel(ps)
    if ~rem(i,250), fprintf(' %d...',i); end
    [k,j] = samplealign([CMZ,LAI],double(ps{i}),'band',1.5,'gap',[0 2],'dist',DF);
    LAI = LAI*.25;
    LAI(k) = LAI(k) + ps{i}(j,2);
    psa{i,1} = [CMZ(k) ps{i}(j,2)];
end

% Align peaks of the glycine sample in the MZ direction
LAI = zeros(size(CMZ));
for i = 1:numel(pg)
    if ~rem(i,250), fprintf(' %d...',i); end
    [k,j] = samplealign([CMZ,LAI],double(pg{i}),'band',1.5,'gap',[0 2],'dist',DF);
    LAI = LAI*.25;
    LAI(k) = LAI(k) + pg{i}(j,2);
    pga{i,1} = [CMZ(k) pg{i}(j,2)];
end

% Peak-preserving resample
[MZs,Ys] = msppresample(psa,5000);
[MZg,Yg] = msppresample(pga,5000);

% Gaussian Filtering
Gpulse = exp(-.5*(-10:10).^2)./sum(exp(-.05*(-10:10).^2));
Ysf = convn(Ys,Gpulse,'same');
Ygf = convn(Yg,Gpulse,'same');

% Visualization
fh3 = msheatmap(MZs,ts,log(Ysf),'resolution',0.15);
title('Serine Treatment (Enhanced)')
axis([480 520 ind_ser'])
fh4 = msheatmap(MZg,tg,log(Ygf),'resolution',0.15);
title('Glycine Treatment (Enhanced)')
axis([480 520 ind_gly'])

 250... 500... 750... 250... 500... 750...

Хроматографическое выравнивание

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

Вы можете использовать SAMPLEALIGN для коррекции дрейфа в хроматографическом домене. Во-первых, вы должны проверить данные и искать худший сдвиг, это помогает оценить BAND ограничение. При панорамировании по обеим тепловым картам можно наблюдать, что общие пептидные пики не сдвигаются более чем на 50 секунд. Использовать входной аргумент SHOWCONSTRAINTS отображение пространства ограничений для операции искажения времени и оценка того, может ли алгоритм динамического программирования (DP) обрабатывать этот размер проблемы. В этом случае имеется менее 12 000 узлов. Опустив выходные аргументы, SAMPLEALIGN отображает только ограничения без выполнения алгоритма точки данных. Также следует отметить, что входными сигналами являются отфильтрованные и улучшенные наборы данных, но они были увеличены до 5000 MZ значений, которые очень требовательны к вычислениям, если вы используете все. Поэтому используйте функцию MSPALIGN для получения уменьшенного списка значений массы/заряда (RMZ), указывающего, где наиболее интенсивные пики, затем используйте SAMPLEALIGN функция также для поиска индексов MZ (или MZg), которые наилучшим образом соответствуют вектору уменьшенной массы/заряда:

RMZ = mspalign([ps;pg])';
idx = samplealign(MZs,RMZ,'width',1); % with these input parameters this
                                      % operation is equivalent to find the
                                      % nearest neighbor for each RMZ in
                                      % MZs.

samplealign([ts Ysf(idx,:)'],[tg Ygf(idx,:)'],'band',50,'showconstraints',true)

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

cc = @(Xu,Yu) (mean(Xu.*Yu,2).^2)./mean(Xu.*Xu,2)./mean(Yu.*Yu,2);
ub = @(X) bsxfun(@minus,X,mean(X,2));
df = @(x,y) 1-cc(ub(x),ub(y));

[i,j] = samplealign([ts Ysf(idx,:)'],[tg Ygf(idx,:)'],'band',50,...
                    'distance',df ,'weight',[false true(1,numel(idx))]);

fh5 = figure;
plot(ts(i),tg(j)); grid
title('Warp Function')
xlabel('Retention Time in Data Set with Serine')
ylabel('Retention Time in Data Set with Glycine')

fh6 = figure; hold on
plot((ts(i)+tg(j))/2,ts(i)-tg(j))
title('Shift Function')
xlabel('Retention Time')
ylabel('Retention Time Shift between Data Sets')

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

[p,p_struct,mu] = polyfit((ts(i)+tg(j))/2,ts(i)-tg(j),20);
sf = @(z) polyval(p,(z-mu(1))./mu(2));
figure(fh6)
plot(tg,sf(tg),'r')
legend('Observed shifts','Estimated shift curve','Location','SouthEast')

Сравнительный анализ

Для проведения сравнительного анализа повторная выборка наборов данных LC/MS на общий временной вектор. При повторной выборке для корректировки измерения времени удержания используется расчетная функция сдвига. В этом примере каждый спектр в обоих наборах данных сдвинут на половину расстояния функции сдвига. В случае множественного совмещения наборов данных можно вычислить парные функции сдвига между всеми наборами данных, а затем применить поправки к общей ссылке таким образом, что общие сдвиги минимизируются [4].

t = (max(min(tg),min(ts)):mean(diff(tg)):min(max(tg),max(ts)))';

% Visualization
fh7 = msheatmap(MZs,t,log(interp1q(ts,Ysf',t+sf(t)/2)'),'resolution',0.15);
title('Serine Treatment (Enhanced & Aligned)')
fh8 = msheatmap(MZg,t,log(interp1q(tg,Ygf',t-sf(t)/2)'),'resolution',0.15);
title('Glycine Treatment (Enhanced & Aligned)')

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

linkaxes(findobj([fh7 fh8],'Tag','MSHeatMap'))
axis([480 520 1550 1900])

Ссылки

[1] Листгартен, Дж. и Эмили, А., «Статистические и вычислительные методы сравнительного протеомного профилирования с использованием жидкостной хроматографии тандемной масс-спектрометрии», Molecular and Cell Proteomics, 4 (4): 419-34, 2005.

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

[3] Prince, J.T. and Marcotte, E.M., «Хроматографическое выравнивание наборов данных протеомики ESI-LC-MS по упорядоченному биективному интерполированному деформированию», Analytical Chemistry, 78 (17): 6140-52, 2006.

[4] Андраде Л. и Манолакос Е. С., «Автоматическая оценка коэффициентов сдвига подвижности в ДНК-хроматограммах», Neural Networks for Signal Proceedings, 91-100, 2003.