Введение

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

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

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

Этот пример использует две выборки в PAe000153 и PAe000155 доступный из Атласа Пептида [2]. Выборки являются сканами LC-ESI-MS четырех соленых частей белка от saccharomyces cerevisiae каждый содержащий больше чем 1 000 пептидов. Выборки дрожжей были обработаны с различными химикатами (глицин и серин) для того, чтобы получить две биологически разнообразных выборки. Выравнивание времени этих двух наборов данных является одним из самых сложных случаев, сообщил в [3]. Наборы данных не распределяются с MATLAB®. Необходимо загрузить наборы данных на локальную директорию или собственный репозиторий. В качестве альтернативы можно попробовать другие наборы данных, доступные в общедоступных базах данных для данных о белке, таких как Репозиторий данных сашими. Если вы получаете какие-либо ошибки, связанные с памятью или пространством "кучи" Java, попытайтесь увеличить свое пространство "кучи" Java, аналогичное описанному здесь. Анализ данных LC/MS требует расширенных объемов памяти от операционной системы; если вы получаете "Out of memory" ошибки при выполнении этого примера, попытайтесь увеличить виртуальную память (или область подкачки) операционной системы или попытайтесь установить переключатель 3GB (только 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 функция, чтобы передискретизировать наборы данных при сохранении пиковых высот и местоположений в размерности массы/заряда. Наборы данных передискретизируются, чтобы иметь обоих общая сетка с 5 000 значений массы/заряда. Общая сетка желательна для сравнительной визуализации, и для дифференциального анализа.

[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              

Чтобы наблюдать ту же необходимую область в обоих наборах данных, необходимо вычислить соответствующие индексы строки в каждой матрице. Например, смотреть пептид достигает максимума в области значений Da MZ 480-520 и области значений времени задержания 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'])

Даже при том, что вы увеличили масштаб та же область значений, можно все еще заметить, что верхний правый пептид в осях смещен в измерении времени задержания. В выборке отнесся с серином, центр этого пика, кажется, происходит приблизительно в 1 595 секунд, в то время как в выборке отнесся с глицином предполагаемое, тот же пептид происходит приблизительно в 1 630 секунд. Это предотвратит вас от точного сравнительного анализа, даже если вы будете передискретизировать наборы данных к тому же временному вектору. В дополнение к сдвигу во время задержания набор данных, кажется, неправильно калибруется в размерности массы/заряда, потому что 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. Также обратите внимание, что входные сигналы являются отфильтрованными и расширенными наборами данных, но они были сверхдискретизированы к 5,000 значений MZ, которые очень в вычислительном отношении требовательны, если вы используете все. Поэтому используйте функциональный MSPALIGN чтобы получить уменьшаемый список значений массы/заряда (RMZ), указание, где самый интенсивный peaks, затем использует SAMPLEALIGN функционируйте также, чтобы найти индексы MZs (или 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 интенсивности, а не низкой ионной интенсивностью шумный 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. и Emili, A., "Статистические и вычислительные методы для сравнительного протеомного профилирования с помощью тандемной масс-спектрометрии жидкостной хроматографии", Молекулярный и Протеомика Ячейки, 4 (4):419-34, 2005.

[2] Desiere, F., и др., "Проект Атласа Пептида", Исследование Нуклеиновых кислот, 34:D655-8, 2006.

[3] Принц, Дж.Т. и Маркотт, E.M., "Хроматографическое выравнивание наборов данных протеомики ESI-LC-MS упорядоченным bijective, интерполированным, деформируясь", Аналитическая Химия, 78 (17):6140-52, 2006.

[4] Андрэйд Л. и Мэнолэкос Э.С., "Автоматическая оценка коэффициентов сдвига мобильности в хроматограммах ДНК", нейронные сети для продолжений обработки сигналов, 91-100, 2003.