Этот пример показывает, как управлять, предварительно обработать и визуализировать данные из Жидкостной хроматографии вместе с Масс-спектрометрией (LC/MS). Эти большие и высокие размерные наборы данных экстенсивно используются в протеомике и metabolomics исследовании. Визуализация комплексного пептида или смесей метаболита предоставляет интуитивный метод, чтобы оценить демонстрационное качество. Кроме того, методическое исправление и предварительная обработка могут привести к автоматизированному высокому анализу пропускной способности выборок, позволяющих точную идентификацию значительных метаболитов и определенных функций пептида в биологической выборке.
В типичном написанном через дефис эксперименте масс-спектрометрии белки и метаболиты получены от ячеек, тканей, или биологических жидкостей, расторгли и денатурировали в решении, и ферментативным образом переварили в смеси. Эти смеси затем разделяются или Высокоэффективной жидкостной хроматографией (HPLC), капиллярным электрофорезом (CE) или газовой хроматографией (GC) и связываются с методом идентификации масс-спектрометрии, таким как Масс-спектрометрия Ионизации Электроспрея (ESI-MS), матрица помогла ионизации (MALDI или TOF-MS SELDI), или тандемная масс-спектрометрия (MS/MS).
В данном примере мы используем тестовую выборку набор данных LC-ESI-MS с семью соединениями белка. Данные в этом примере (7MIX_STD_110802_1) от Репозитория данных сашими. Набор данных не распределяется с MATLAB®. Чтобы завершить этот пример, необходимо загрузить набор данных в локальную директорию или собственный репозиторий. Также можно попробовать другие наборы данных, доступные в других общедоступных базах данных для данных о выражении белка, таких как Атлас Пептида в Институте Системной биологии [1].
Большинство текущих массовых спектрометров может перевести или сохранить данные о приобретении с помощью mzXML схемы. Этот стандарт является XML (расширяемый язык разметки) - базирующийся формат общего файла, разработанный проектом сашими обратиться к проблемам, вовлеченным в представление наборов данных от различных производителей и от различных экспериментальных настроек в общую и расширяемую схему. файлы 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: 'false' DataProcessingDeisotoped: 'N/A' DataProcessingChargeDeconvoluted: 'N/A' DataProcessingSpotIntegration: 'N/A' NumberOfMSLevels: 2
Функция MZXMLREAD
читает XML-документ в структуру MATLAB. Поля scan
и index
помещаются в первый уровень выходной структуры для улучшенного доступа к спектральным данным. Остаток от mzXML дерева документов анализируется согласно спецификациям схемы. Этот набор данных LC/MS содержит 7 161 сканирование с двумя уровнями MS. Для этого примера вы будете использовать только первые сканирования уровня. Вторые спектры уровня обычно используются для идентификации пептида/белка и прибывают в более поздний этап в некоторых типах исследований рабочего процесса. MZXMLREAD
может отфильтровать желаемые сканирования, не загружая весь набор данных в память:
mzXML_struct = mzxmlread('7MIX_STD_110802_1.mzXML','LEVEL',1)
Starting to parse document... Building mzXML substructure... Building scan substructure... Building index substructure... DONE! 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 от каждого сканирования в массив ячеек (peaks]
) и их соответствующее время задержания в вектор-столбец (time
). Можно извлечь спектры определенного уровня путем установки параметра входа LEVEL
.
[peaks,time] = mzxml2peaks(mzXML_struct); numScans = numel(peaks)
numScans = 2387
Функция 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 14050896 cell peaks_fil 2387x1 308664 cell
Можно настроить 3-D обзор отфильтрованного peaks с помощью функции 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)
Общие методы в промышленности работают с пиковой информацией (иначе центроидные данные) вместо необработанных сигналов. Это может сохранить память, но некоторые важные детали не видимы, особенно когда необходимо осмотреть выборки с комплексными смесями. Чтобы далее анализировать этот набор данных, мы можем создать общую сетку в размерности массы/заряда. С тех пор не все сканирования имеют достаточно информации, чтобы восстановить исходный сигнал, мы используем пиковый метод передискретизации сохранения. Путем выбора соответствующих параметров для функции MSPPRESAMPLE
можно гарантировать, что разрешение спектров не потеряно, и что максимальные значения peaks коррелируют к исходной пиковой информации.
[MZ,Y] = msppresample(peaks,5000); whos('MZ','Y')
Name Size Bytes Class Attributes MZ 5000x1 20000 single Y 5000x2387 47740000 single
С этой матрицей ионной интенсивности, Y
, можно создать цветную карту тепла. Функция MSHEATMAP
автоматически настраивает шкалу палитры, используемую, чтобы показать статистически значительный peaks с горячими цветами и шумный peaks с холодными цветами. Алгоритм основан на кластеризирующемся значительном peaks и шумном peaks путем оценки смеси Gaussians с подходом Максимизации Ожидания. Кроме того, можно использовать параметр входа 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, реплицируют профили, как, сообщают, колеблются от 0,1% до 0,3% [3]. Такая изменчивость может быть вызвана несколькими факторами, например, плохой калибровкой детектора, низкого отношения сигнал-шум или проблем в пиковых алгоритмах экстракции. Функция 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])
MSPALIGN
вычисляет один CMZ для целого набора данных LC/MS. Это не может быть идеальным случаем для выборок с более комплексными смесями пептидов и/или метаболитов, чем набор данных, используемый в этом примере. В случае комплексных смесей можно выровнять каждый спектр к локальному набору спектров, которые содержат только информативный peaks (высокая интенсивность) с подобными временами задержания, в противном случае калибровочный процесс в областях с маленьким peaks (низкая интенсивность) может быть смещен другим peaks, который совместно использует подобные значения массы/заряда, но является в различные времена задержания. Чтобы выполнить более прекрасную калибровку, можно использовать функцию SAMPLEALIGN
. Эта функция является обобщением алгоритмов Ограниченного динамического времени, деформируясь (CDTW), обычно используемых в речи, обрабатывающей [6]. В следующем цикле for
мы поддерживаем буфер с интенсивностью предыдущих выровненных спектров (LAI). Ионная интенсивность спектров масштабируется с анонимной функцией SF
(в SAMPLEALIGN
), чтобы уменьшить расстояние между большим peaks. SAMPLEALIGN
уменьшает общее расстояние всех совпадающих точек и вводит разрывы по мере необходимости. В этом случае мы используем более прекрасный вектор MZ (FMZ), такой, что мы сохраняем правильное значение массы/заряда peaks как можно больше. Примечание: это может занять время, когда алгоритм CDTW выполняется 2,387 раз.
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. и др., "Проект Атласа Пептида", Исследование Нуклеиновых кислот, 34:D655-8, 2006.
[2] Purvine, S., Kolker, N. и Kolker, E., "Спектральная качественная оценка для тандемной протеомики масс-спектрометрии Высокой Пропускной способности", OMICS: журнал интегральной биологии, 8 (3):255-65, 2004.
[3] Kazmi, A.S., и др., "Выравнивание спектров массы высокого разрешения: Разработка эвристического подхода для metabolomics", Metabolomics, 2 (2):75-83, 2006.
[4] Jeffries, N., "Алгоритмы для выравнивания масс-спектрометрии протеомные данные", Биоинформатика, 21 (14):3066-3073, 2005.
[5] Ю, W., и др., "Несколько достигают максимума выравнивание в последовательном анализе данных: основанный на пробеле шкалы подход", Сделка IEEE®/ACM Вычислительная Биология и Биоинформатика, 3 (3):208-219, 2006.
[6] Sakoe, H. и Чиба s., "Динамическая оптимизация алгоритма программирования для распознавания произносимого слова", Сделка IEEE. Акустика, Речь и Обработка сигналов, ASSP-26 (1):43-9, 1978.