В этом примере показано, как улучшить качество необработанных данных о масс-спектрометрии. В частности, этот пример иллюстрирует типичные шаги для предварительной обработки белка улучшенный поверхностью лазер desorption/ionization-time спектров массы рейса (SELDI-TOF).
Данные о масс-спектрометрии могут храниться в различных форматах. Если данные хранятся в текстовых файлах с двумя столбцами (масса/заряд (M/Z) отношения и соответствующие значения интенсивности), можно использовать одну из следующих функций MATLAB® I/O: importdata
dlmread
, или textscan
. В качестве альтернативы, если данные хранятся в отформатированных файлах JCAMP-DX, можно использовать функциональный jcampread
. Если данные содержатся в электронной таблице рабочей книги Excel®, можно использовать функциональный xlsread
.If данные хранятся в 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
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)
msalign
| msbackadj
| msheatmap
| msnorm
| mspeaks
| msresample
| mssgolay
| msviewer