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

Данные масс-спектрометрии могут храниться в различных форматах. Если данные хранятся в текстовых файлах с двумя столбцами (отношение масса/заряд (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, сохраняя информационное содержимое спектров. The msresample функция позволяет вам выбрать новый вектор M/Z, а также применяет фильтр antialias, который препятствует сворачиванию высокочастотного шума на более низкие частоты.

Загрузите спектр с высоким разрешением, взятый из набора данных 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])

Коррекция базового уровня

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

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

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

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

Выкалибрация масс-спектрометра приводит к изменениям зависимости между наблюдаемым вектором 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')

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

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

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

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

The 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}

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

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

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