Введение

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

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

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

Большинство текущих масс-спектрометров могут перемещать или сохранять данные сбора с помощью схемы mzXML. Этот стандарт представляет собой основанный на XML (eXtensible Разметки Language) общий формат файла, разработанный в рамках проекта Sashimi для решения проблем, связанных с представлением наборов данных от различных производителей и от различных экспериментальных настроек до общей и расширяемой схемы. Файлы mzXML, используемые в расставленной масс-спектрометрии, обычно являются очень большими. The 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

The MZXMLREAD функция считывает XML-документ в структуру MATLAB. Поля scan и index размещены на первом уровне структуры output для улучшения доступа к спектральным данным. Остальная часть дерева документов 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]

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

Более подробная информация, касающаяся масс-спектрометра и экспериментальных условий, найдена в поле 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 из каждого скана в массив ячеек (peaks]) и их соответствующее время удержания в вектор-столбец (time). Можно извлечь спектры определенного уровня, установив LEVEL входной параметр.

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

numScans =

        2387

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

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

Можно также фильтровать 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 обзор отфильтрованного peaks с помощью STEM3 функция. The 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 функция, можно гарантировать, что разрешение спектров не потеряно, и что максимальные значения peaks коррелируют с исходной пиковой информацией.

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

  Name         Size                 Bytes  Class     Attributes

  MZ        5000x1                  20000  single              
  Y         5000x2387            47740000  single              

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

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);

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

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

[CMZ, peaks_CMZ] = mspalign(peaks);

Повторите процесс визуализации с выровненным 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])

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

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

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

После сглаживания и повторной дискретизации набора данных LC/MS в регулярную сетку можно извлечь наиболее информативные спектры, рассмотрев локальные максимумы хроматограммы Total Ion (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., «The Peptide Atlas Project», Nucleic Acids Research, 34:D655-8, 2006.

[2] Purvine, S., Kolker, N. and Kolker, E., «Spectral Quality Assessment for High-Through Tandem Mass Spectromics Proteomics», 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. вычислительная биология и биоинформатика, 3 (3): 208-219, 2006.

[6] Sakoe, H. and Chiba s., «Dynamic programming algorithm optimization for spoken word recognition», IEEE Trans. Acoustics, Speech and Signal Processing, ASSP-26 (1): 43-9, 1978.