Введение

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

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

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

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

Считывайте и извлекайте списки peaks из 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 Da 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 секунд. Это предотвратит точный сравнительный анализ, даже если вы повторно отобразите наборы данных на один и тот же временной вектор. В дополнение к сдвигу во времени удержания набор данных, по-видимому, ненадлежащим образом калибруется по размерности масса/заряд, потому что peaks не имеют компактной формы в смежных спектрах.

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

Перед корректировкой времени удержания можно усилить выборки с помощью подхода, аналогичного тому, который описан в примере Визуализация и предварительная обработка переносимых наборов данных масс-спектрометрии для профилирования метаболитов и белков/пептидов. Для краткости отображаем только код 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 ограничение. Путем панорамирования по обеим тепловым картам можно заметить, что общий пептидный peaks не сдвигаются более чем на 50 секунд. Используйте входной параметр SHOWCONSTRAINTS отобразить пространство ограничений для операции деформации во времени и оценить, может ли алгоритм динамического программирования (DP) обработать этот размер задачи. В этом случае у вас менее 12 000 узлов. Опуская выходные аргументы, SAMPLEALIGN отображает только ограничения без выполнения алгоритма DP. Также обратите внимание, что входные сигналы являются фильтрованными и расширенными наборами данных, но они были увеличены до 5000 МЗц значений, которые очень требовательны к вычислениям, если вы используете все. Поэтому используйте функцию MSPALIGN чтобы получить сокращенный список значений (RMZ), указывающий, где наиболее интенсивнейший peaks, используйте 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 позволяет вам задать свою собственную метрику для вычисления расстояния между спектрами, также возможно предусмотреть метрику, которая вознаграждает больше пар спектров, которые совпадают с пиками высокой интенсивности ионов, чем с шумным peaks низкой интенсивности ионов. Используйте входной параметр 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] Listgarten, J. and Emili, A., «Statistical and computational methods for comparative proteomic profiling with liquid chromatography tandem massectrometry», Molecular and Cell Proteomics, 4 (4): 419-34, 2005.

[2] Desiere, F., et al., «The Peptide Atlas Project», 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] Andrade L. and Manolakos E.S., «Автоматическая оценка коэффициентов сдвига мобильности в хроматограммах ДНК», Нейронные сети для обработки сигналов, 91-100, 2003.