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

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

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

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

Выполните выборку спектров до 10000 точек M/Z между 2000 и 11000. Используйте SHOWPLOT для создания настраиваемого графика, позволяющего отслеживать и оценивать качество операции предварительной обработки.

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

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

axis([3875 3895 0 90])

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Уменьшение шума при сохранении пиковых значений

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

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

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

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

axis([8000 9000 -1 8])

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

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

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

P1 =

  4x1 cell array

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

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

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

P2 =

  4x1 cell array

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

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

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

P3 =

  4x1 cell array

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

Binning: пиковое слияние с помощью иерархической кластеризации

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

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

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

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

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

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

Динамическое программирование Binning

samplealign функция позволяет использовать алгоритм динамического программирования для назначения наблюдаемых пиков в каждой спектрограмме общему опорному вектору масса/заряд (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)