Загрузка данных

Данные о масс-спектрометрии могут храниться в различных форматах. Если данные хранятся в текстовых файлах с двумя столбцами (масса/заряд (M/Z) отношения и соответствующие значения интенсивности), можно использовать одну из следующих функций MATLAB® I/O: importdatadlmread, или textscan. В качестве альтернативы, если данные хранятся в отформатированных файлах JCAMP-DX, можно использовать функциональный jcampread. Если данные содержатся в электронной таблице рабочей книги Excel®, можно использовать функциональный xlsread.If данные хранятся в mzXML, отформатировал файлы, можно использовать функциональный mzxmlread, и наконец, если данные хранятся в SPC, отформатировал файлы, можно использовать tgspcread.

Этот пример использует спектрограммы, взятые из одного из рака яичника с низкой разрешающей способностью наборы данных NCI/FDA. Эти спектры были сгенерированы с помощью связывающего белок чипа WCX2, два с ручной демонстрационной обработкой и два с автоматизированным демонстрационным фармацевтом/процессором. Данные могут быть загружены с FDA-NCI Клинический Банк данных Программы Протеомики.

sample = importdata('mspec01.csv')

sample = 

  struct with fields:

          data: [15154x2 double]
      textdata: {'M/Z'  'Intensity'}
    colheaders: {'M/Z'  'Intensity'}

Отношения M/Z находятся в первом столбце data поле и ионная интенсивность находятся во втором.

MZ = sample.data(:,1);
Y  = sample.data(:,2);

Для лучшей манипуляции данных можно загрузить несколько спектрограмм и конкатенировать их в одну матрицу. Используйте dlmread функционируйте, чтобы читать, запятая разделила файлы значения. Примечание: Этот пример принимает, что отношения M/Z являются тем же самым для этих четырех файлов. Для наборов данных с различными отношениями M/Z используйте msresample создать универсальный вектор M/Z.

files = {'mspec01.csv','mspec02.csv','mspec03.csv','mspec04.csv'};

for i = 1:4
    Y(:,i) = dlmread(files{i},',',1,1); % skips the first row (header)
end

Используйте plot команда, чтобы смотреть загруженные спектрограммы.

plot(MZ,Y)
axis([0 20000 -20 105])
xlabel('Mass/Charge (M/Z)')
ylabel('Relative Intensity')
title('Four Low-Resolution Mass Spectrometry Examples')

Передискретизация спектров

Передискретизация данных о масс-спектрометрии имеет несколько преимуществ. Это гомогенизирует массу/заряд (M/Z) вектор, позволяя вам сравнить различные спектры под той же ссылкой и в том же разрешении. В наборах данных с высоким разрешением большой размер файлов приводит к в вычислительном отношении интенсивным алгоритмам. Однако спектры с высоким разрешением могут быть избыточными. Путем передискретизации можно десятикратно уменьшить сигнал в более управляемый вектор M/Z, сохранив информационное содержимое спектров. msresample функция позволяет вам выбирать новый вектор M/Z и также применяет фильтр сглаживания, который препятствует тому, чтобы высокочастотный шум свернулся в более низкие частоты.

Загрузите спектр с высоким разрешением, взятый из рака яичника с высоким разрешением набор данных NCI/FDA. Для удобства спектр включен в отформатированный MAT файл.

load sample_hi_res
numel(MZ_hi_res)

ans =

      355760

Downsample спектры к 10,000 M/Z указывает между 2 000 и 11,000. Используйте SHOWPLOT свойство создать индивидуально настраиваемый график, который позволяет вам следовать и оценить качество действия предварительной обработки.

[MZH,YH] = msresample(MZ_hi_res,Y_hi_res,10000,'RANGE',[2000 11000],'SHOWPLOT',true);

Изменение масштаба в уменьшаемую область показывает деталь процедуры субдискретизации.

axis([3875 3895 0 90])

Базовая коррекция

Данные о масс-спектрометрии обычно показывают различную базовую линию, вызванную химическим шумом в матрице или ионной перегрузкой. msbackadj функционируйте оценивает низкочастотную базовую линию, которая скрыта среди высокочастотного шума и peaks сигнала. Это затем вычитает базовую линию из спектрограммы.

Настройте базовую линию набора спектрограмм и покажите только вторую и его предполагаемое образование.

YB = msbackadj(MZ,Y,'WINDOWSIZE',500,'QUANTILE',0.20,'SHOWPLOT',2);

Спектральное выравнивание профилей

Miscalibration массового спектрометра приводит к изменениям отношения между наблюдаемым вектором M/Z и истинное время рейса ионов. Поэтому систематические сдвиги могут появиться в повторных экспериментах. Когда известный профиль peaks ожидается в спектрограмме, можно использовать функциональный msalign стандартизировать значения M/Z.

Чтобы выровнять спектрограммы, обеспечьте набор значений M/Z, где ссылочный peaks, как ожидают, появится. Можно также задать вектор с относительными весами, который используется выравнивающимся алгоритмом, чтобы подчеркнуть peaks с небольшой площадью.

P = [3991.4 4598 7964 9160]; % M/Z location of reference peaks
W = [60 100 60 100];         % Weight for reference peaks

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

msheatmap(MZ,YB,'MARKERS',P,'RANGE',[3000 10000])
title('Before Alignment')

Выровняйте набор спектрограмм к ссылочному данному peaks.

YA = msalign(MZ,YB,P,'WEIGHTS',W);
msheatmap(MZ,YA,'MARKERS',P,'RANGE',[3000 10000])
title('After Alignment')

Нормализация

В повторных экспериментах распространено найти систематические различия в общей сумме выделенных и ионизированных белков. msnorm функционируйте реализует несколько изменений типичной нормализации (или стандартизация) методы.

Например, один из многих методов, чтобы стандартизировать значения спектрограмм должен перемасштабировать максимальную интенсивность каждого сигнала к определенному значению, например, 100. Также возможно проигнорировать проблематичные области; например, в выборках сыворотки вы можете хотеть проигнорировать область малой массы (M/Z <1 000 дальтонов.).

YN1 = msnorm(MZ,YA,'QUANTILE',1,'LIMITS',[1000 inf],'MAX',100);
figure
plot(MZ,YN1)
axis([0 10000 -20 150])
xlabel('Mass/Charge (M/Z)')
ylabel('Relative Intensity')
title('Normalized to the Maximum Peak')

msnorm функция может также стандартизировать при помощи области под кривой (AUC) и затем перемасштабировать спектрограммы, чтобы иметь относительную интенсивность ниже 100.

YN2 = msnorm(MZ,YA,'LIMITS',[1000 inf],'MAX',100);
figure
plot(MZ,YN2)
axis([0 10000 -20 150])
xlabel('Mass/Charge (M/Z)')
ylabel('Relative Intensity')
title('Normalized Using the Area Under the Curve (AUC)')

Пиковое шумоподавление сохранения

Стандартизированные спектры обычно содержат смесь шума и сигнала. Некоторые приложения требуют, чтобы шумоподавление спектрограмм улучшило валидность и точность наблюдаемых значений массы/заряда peaks в спектрах. По той же причине шумоподавление также улучшает дальнейшие пиковые алгоритмы обнаружения. Однако важно сохранить резкость (или высокочастотные компоненты) пика как можно больше. Для этого можно использовать сглаживание Lowess (mslowess) и полиномиальные фильтры (mssgolay).

Сглаживайте спектрограммы с полиномиальным фильтром второго порядка.

YS = mssgolay(MZ,YN2,'SPAN',35,'SHOWPLOT',3);

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

axis([8000 9000 -1 8])

Пиковое открытие с шумоподавлением вейвлетов

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

P1 = mspeaks(MZ,YS,'DENOISING',false,'HEIGHTFILTER',2,'SHOWPLOT',1)

P1 =

  4x1 cell array

    {164x2 double}
    {171x2 double}
    {169x2 double}
    {147x2 double}

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

P2 = mspeaks(MZ,YN2,'BASE',12,'MULTIPLIER',10,'HEIGHTFILTER',1,'SHOWPLOT',1)

P2 =

  4x1 cell array

    {322x2 double}
    {370x2 double}
    {324x2 double}
    {295x2 double}

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

P3 = cellfun( @(x) x(x(:,1)>1500,:),P2,'UNIFORM',false)

P3 =

  4x1 cell array

    {81x2 double}
    {93x2 double}
    {57x2 double}
    {53x2 double}

Раскладывание: пиковое объединение иерархической кластеризацией

О peaks, соответствующем подобным составным объектам, можно все еще сообщить с небольшими различиями в массе/заряде или дрейфами. Предположение, что эти четыре спектрограммы соответствуют сопоставимым биологическим/химическим выборкам, может быть полезно сравнить peaks от различных спектров, который требует пикового раскладывания (a.k.a. объединение пика). Решающая задача в раскладывании данных состоит в том, чтобы создать общий вектор ссылки массы/заряда (или интервалы). Идеально, интервалы должны собрать один пик из каждого сигнала и должны постараться не собирать несколько соответствующих peaks из того же сигнала в тот же интервал.

Этот пример использует иерархическую кластеризацию, чтобы вычислить общий вектор ссылки массы/заряда. Подход достаточен при использовании спектров с низкой разрешающей способностью; однако, для спектров с высоким разрешением или для наборов данных со многими спектрограммами, функциональным mspalign предоставляет другие масштабируемые методы, чтобы оценить общую ссылку массы/заряда и выполнить раскладывание данных.

Поместите весь peaks в единый массив и создайте вектор с индексом спектрограммы для каждого пика.

allPeaks = cell2mat(P3);
numPeaks = cellfun(@(x) length(x),P3);
Sidx = accumarray(cumsum(numPeaks),1);
Sidx = cumsum(Sidx)-Sidx;

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

distfun = @(x,y) (x(:,1)-y(:,1)).^2 + (x(:,2)==y(:,2))*10^6

tree = linkage(pdist([allPeaks(:,1),Sidx],distfun));
clusters = cluster(tree,'CUTOFF',75,'CRITERION','Distance');

distfun =

  function_handle with value:

    @(x,y)(x(:,1)-y(:,1)).^2+(x(:,2)==y(:,2))*10^6

Общий вектор ссылки массы/заряда (CMZ) найден путем вычисления центроидов для каждого кластера.

CMZ = accumarray(clusters,prod(allPeaks,2))./accumarray(clusters,allPeaks(:,2));

Точно так же максимальная пиковая интенсивность каждого кластера также вычисляется.

PR = accumarray(clusters,allPeaks(:,2),[],@max);

[CMZ,h] = sort(CMZ);
PR = PR(h);

figure
hold on
box on
plot([CMZ CMZ],[-10 100],'-k')
plot(MZ,YN2)
axis([7200 8500 -10 100])
xlabel('Mass/Charge (M/Z)')
ylabel('Relative Intensity')
title('Common Mass/Charge (M/Z) Locations found by Clustering')

Динамическое раскладывание программирования

samplealign функция позволяет вам использовать динамический алгоритм программирования, чтобы присвоить наблюдаемый peaks в каждой спектрограмме к общему вектору ссылки массы/заряда (CMZ).

Когда более простые подходы раскладывания используются, такие как округление значений массы/заряда или использование самого близкого соседнего квантования к CMZ вектор, тот же пик от различных спектров мой быть присвоенным различным интервалам из-за маленьких дрейфов, которые все еще существуют. Чтобы обойти эту проблему, размер интервала может быть увеличен с жертвой разрешения пика масс-спектрометрии. При помощи динамического раскладывания программирования вы сохраняете разрешение при минимизации проблемы присвоения подобных составных объектов от различных спектрограмм до различных пиковых местоположений.

PA = nan(numel(CMZ),4);
for i = 1:4
    [j,k] = samplealign([CMZ PR],P3{i},'BAND',15,'WEIGHTS',[1 .1]);
    PA(j,i) = P3{i}(k,2);
end

figure
hold on
box on
plot([CMZ CMZ],[-10 100],':k')
plot(MZ,YN2)
plot(CMZ,PA,'o')
axis([7200 8500 -10 100])
xlabel('Mass/Charge (M/Z)')
ylabel('Relative Intensity')
title('Peaks Aligned to the Common Mass/Charge (M/Z) Reference')

Используйте msviewer смотреть предварительно обработанные спектрограммы на данной области значений (например, между значениями 7600 и 8200).

r1 = 7600;
r2 = 8200;
range = MZ > r1 & MZ < r2;
rangeMarkers = CMZ(CMZ > r1 & CMZ < r2);

msviewer(MZ(range),YN2(range,:),'MARKERS',rangeMarkers,'GROUP',1:4)