识别重要的蛋白质特性和分类资料
这个例子展示了如何将质谱数据和使用一些统计工具来寻找潜在的疾病诊断标记和蛋白质组模式。
介绍
血清蛋白质组学模式诊断可以用来区分样本和无疾病患者。配置文件生成模式应用表面增强激光解吸电离(SELDI)蛋白质质谱分析。这种技术有潜力改善癌症病理临床诊断测试。目标是选择一组减少测量或“特征”,可以用来区分癌症病人和控制。这些特性将在特定的质量/电荷离子强度水平值。
数据进行预处理
卵巢癌的数据集在这个例子中是FDA-NCI临床蛋白质组学数据库程序。生成数据集使用WCX2蛋白质阵列。数据集包括95控制和121卵巢癌。这个数据集的详细描述,请参见[1]和[4]。
这个示例假设您已经有了数据预处理OvarianCancerQAQCdataset.mat。然而,如果你没有数据文件,您可以通过下面的步骤创建的例子批处理的光谱使用顺序和并行计算。
或者,您可以运行脚本msseqprocessing.m。
目录(fullfile (matlabroot,“例子”,“bioinfo”,“主要”))%确保支持文件的搜索路径金宝app类型msseqprocessing
% MSSEQPROCESSING脚本创建OvarianCancerQAQCdataset。垫(用于% CANCERDETECTDEMO)。在运行这个文件之前初始化变量%“库”你放置你质谱仪%文件的完整路径。例如:% %库= ' F: / MassSpecRepository / OvarianCD_PostQAQC / ';% %或% %库= ' / home /用户名/ MassSpecRepository / OvarianCD_PostQAQC / ';% %的近似时间执行(奔腾4,4 ghz) 18分钟。如果你有并行计算工具箱指BIODISTCOMPDEMO看到%如何加速这一分析。% 2003 - 2008版权MathWorks公司repositoryC =(库“癌症/”);repositoryN =(库“正常/ ');filesCancer = dir ([repositoryC ' * . txt ']);NumberCancerDatasets =元素个数(filesCancer); fprintf('Found %d Cancer mass-spectrograms.\n',NumberCancerDatasets) filesNormal = dir([repositoryN '*.txt']); NumberNormalDatasets = numel(filesNormal); fprintf('Found %d Control mass-spectrograms.\n',NumberNormalDatasets) files = [ strcat('Cancer/',{filesCancer.name}) ... strcat('Normal/',{filesNormal.name})]; N = numel(files); % total number of files fprintf('Total %d mass-spectrograms to process...\n',N) [MZ,Y] = msbatchprocessing(repository,files); disp('Finished; normalizing and saving to OvarianCancerQAQCdataset.mat.') Y = msnorm(MZ,Y,'QUANTILE',0.5,'LIMITS',[3500 11000],'MAX',50); grp = [repmat({'Cancer'},size(filesCancer));... repmat({'Normal'},size(filesNormal))]; save OvarianCancerQAQCdataset.mat Y MZ grp
从脚本的预处理步骤和上面列出的例子是为了说明一组代表可能的预处理过程。使用不同的步骤或参数可能导致不同的结果和可能的改进这个例子。
加载数据
一旦预处理数据,您可以加载到MATLAB。
负载OvarianCancerQAQCdataset谁
类属性名称大小字节MZ 120000 x1 15000双Y 15000 x216 15000双grp 216 x1 25056细胞
有三个变量:MZ,Y,grp。MZ是质量/电荷向量,Y所有216名患者的强度值(控制和癌症),然后呢grp持有这些样品代表的索引信息,癌症患者和哪些代表正常的病人。
初始化一些变量,将通过使用示例。
N =元素个数(grp);%的样本数量Cidx = strcmp (“癌症”、grp);%逻辑索引向量对癌症样本Nidx = strcmp (“正常”、grp);%逻辑索引向量正常样本Cvec =找到(Cidx);%指数向量对癌症样本Nvec =找到(Nidx);%索引向量正常样本xAxisLabel =“质量/电荷(M / Z)”;% x标签的情节yAxisLabel =离子强度的;% y标签的情节
可视化的一些样品
你可以画出一些数据集到图窗口直观地比较两组资料;在本例中5色从癌症患者(蓝色)和五个控制病人显示(绿色)。
图;持有在;hC =情节(MZ, Y (:, Cvec (1:5)),“b”);hN =情节(MZ, Y (:, Nvec (1:5)),‘g’);包含(xAxisLabel);ylabel (yAxisLabel);轴([2000 12000 5 60])传说([环(1),hC (1)), {“控制”,卵巢癌的})标题(多个样品谱图的)
放大该地区从8500年到8700年M / Z显示了一些山峰可能对分类有用的数据。
轴([8450、8700、1、7])
可视化整个数据集的另一种方法是看组平均信号控制和癌症样本。您可以绘制组每组的平均和信封。
mean_N =意味着(Y (:, Nidx), 2);%组织控制样本的平均水平max_N = max (Y (:, Nidx), [], 2);%控制样本的信封min_N = min (Y (:, Nidx), [], 2);%控制样品的底部信封mean_C =意味着(Y (:, Cidx), 2);%组织癌症样本的平均水平max_C = max (Y (:, Cidx), [], 2);%控制样本的信封min_C = min (Y (:, Cidx), [], 2);%控制样品的底部信封图;持有在;hC =情节(MZ mean_C,“b”);hN =情节(MZ mean_N,‘g’);gC =情节(MZ, [max_C min_C),“b——”);gN =情节(MZ, [max_N min_N),“g——”);包含(xAxisLabel);ylabel (yAxisLabel);轴([8450、8700、1、7])传说([hN, hC, gN (1) gC (1)), {对照组Avg。,“Avg卵巢癌组”,…“控制信封”,卵巢癌信封的},…“位置”,“西北”)标题(“组平均和组信封”)
观察到明显没有单一的功能,能辨别两组。
排名关键特性
发现重要特性的一种简单的方法是假设每个M / Z值是独立的,计算一个双向的学习任务。rankfeatures返回一个指数最重要的M / Z值,例如100年指数排名的检验统计量的绝对值。该特征选择方法也被称为过滤方法,学习算法不涉及到如何选择的特性。
[壮举,统计]= rankfeatures (Y, grp,“标准”,的tt,“数量”,100);
第一个输出的rankfeatures可以用来提取M / Z值的显著特征。
sig_Masses = MZ(成绩);sig_Masses (1:7)%显示前七
ans = 1.0 e + 03 * 8.1009 8.1016 8.1024 8.1001 8.1032 7.7366 - 7.7359
第二输出rankfeatures是一个向量的检验统计量的绝对值。你可以画出它的光谱yyaxis。
图;yyaxis左情节(MZ (mean_N mean_C]);ylim ([20]) xlim((7950、8300))标题(“重要的M / Z值”)包含(xAxisLabel);ylabel (yAxisLabel);yyaxis正确的情节(MZ,统计);ylim ([22]) ylabel (检验统计量的);传奇({对照组Avg。,卵巢癌组Avg。,“测试数据”})
注意,有显著的地区高M / Z值但低强度(~ 8100 Da)。其他方法来测量类可分性是可用的rankfeatures基于熵等Bhattacharyya或经验下的面积接受者操作特征(ROC)曲线。
盲分类使用线性判别分析(LDA)
既然你已经确定了一些重要的特性,您可以使用这些信息对癌细胞和正常样本进行分类。由于小数量的样品,您可以运行使用20%交叉验证坚持更好地估计分类器的性能。cvpartition允许你设置不同类型的系统的训练和测试指标评价方法,如抵抗,K-fold Leave-M-Out。
per_eval = 0.20;%交叉验证培训规模rng (“默认”);%初始化随机发生器相同的状态%用于生成发布的例子简历= cvpartition (grp),“坚持”per_eval)
简历=抵抗交叉验证分区NumObservations: 216 NumTestSets: 1 TrainSize: 173 TestSize: 43
观察到的特性选择只有从训练子集和执行验证测试子集。classperf允许您跟踪多个验证。
cp_lda1 = classperf (grp);%初始化CP对象为k = 1:10% 10倍交叉验证运行简历=重新分区(简历);壮举= rankfeatures (Y(:,培训(cv)、grp(培训(cv)),“数量”,100);c =分类(Y(功绩,测试(cv)), Y(成绩、培训(cv)), grp(培训(cv)));classperf (cp_lda1 c测试(cv));%更新CP对象与当前验证结束
循环后你可以评估整体性能的CP盲分类使用任何属性对象,如错误率,灵敏度、特异性等。
cp_lda1
标签:“描述:“ClassLabels: {2 x1细胞}GroundTruth: x1双[216]NumberOfObservations: 216 ControlClasses: 2 TargetClasses: 1 ValidationCounter: 10 SampleDistribution: x1双[216]ErrorDistribution: x1双[216]SampleDistributionByClass: [2 x1双]ErrorDistributionByClass: [2 x1双]CountingMatrix: [3 x2双]CorrectRate: 0.8488 ErrorRate: 0.1512 LastCorrectRate: 0.8837 LastErrorRate: 0.1163 InconclusiveRate: 0 ClassifiedRate: 1敏感性:0.8208特异性:0.8842 PositivePredictiveValue: 0.8995 NegativePredictiveValue: 0.7962 PositiveLikelihood: 7.0890 NegativeLikelihood: 0.2026患病率:0.5581 DiagnosticTable: [2 x2双]
这个天真的特征选择方法可以改善通过消除一些功能基于区域信息。例如,“NWEIGHT”rankfeatures比邻近的M / Z特征的检验统计量,这样可以纳入其他重要M / Z值的子集选择功能
cp_lda2 = classperf (grp);%初始化CP对象为k = 1:10% 10倍交叉验证运行简历=重新分区(简历);壮举= rankfeatures (Y(:,培训(cv)、grp(培训(cv)),“数量”,100,“NWEIGHT”5);c =分类(Y(功绩,测试(cv)), Y(成绩、培训(cv)), grp(培训(cv)));classperf (cp_lda2 c测试(cv));%更新CP对象与当前验证结束cp_lda2.CorrectRate%的平均正确分类率
ans = 0.9023
PCA - LDA降低数据的维数
Lilien等人[2]中给出一个算法来降低数据维数,利用主成分分析(PCA),然后LDA用于分类组。在这个例子中2000 M / Z的最重要的特征空间映射到150主成分
cp_pcalda = classperf (grp);%初始化CP对象为k = 1:10% 10倍交叉验证运行简历=重新分区(简历);%选择2000年最重要的特性。壮举= rankfeatures (Y(:,培训(cv)、grp(培训(cv)),“数量”,2000);%主成分分析降低维数P = pca (Y(成绩、培训(cv)));%项目到PCA空间x = Y(功绩,:)* P (: 1:15);%利用LDAc =分类(x(测试(cv):), x(培训(cv):), grp(培训(cv)));classperf (cp_pcalda c测试(cv));结束cp_pcalda.CorrectRate%的平均正确分类率
ans = 0.9814
随机搜索特征子集的选择
特征选择也可以强化了分类,这种方法通常被称为包装器选择方法。随机搜索特征选择生成随机子集的特性和与学习算法独立评估其质量。后来,它选择的最常见的良好特性。李等人在[3]这个概念应用于蛋白表达模式的分析。的randfeatures函数允许您搜索功能的一个子集用LDA或再分类器在随机特性的子集。
注意:下面的例子是计算密集型,因此它被禁用的例子。也为更好的结果,你应该增加池大小和严格分类器的默认值randfeatures。类型帮助randfeatures为更多的信息。
如果0% < = = 1更改为启用。这可以用大量的时间来完成。简历=重新分区(简历);[壮举,fCount] = randfeatures (Y(:,培训(cv)、grp(培训(cv)),…“分类”,“哒”,“PerformanceThreshold”,0.90);其他的负载randFeatCancerDetect结束
评估选择的质量特性与评价集
第一个输出randfeatures是一组有序的MZ值的指数。第一项发生最频繁的子集分类好。第二个输出的实际计算每个值被选中的次数。您可以使用嘘看看这个分布。
图;嘘(fCount,马克斯(fCount) + 1);
你会发现最值出现在最一次选择的子集。放大的细节给出了一个更好的主意更频繁地选择值。
轴(80 0 100 [0])
只有少数的价值观选择的10倍以上。你可以想象这些通过阻止阴谋显示最常选择的特性。
图;持有在;sigFeats = fCount;sigFeats (sigFeats < = 10) = 0;情节(MZ (mean_N mean_C]);茎(MZ (sigFeats > 0), sigFeats (sigFeats > 0),“r”);轴([2000、12000、80])传说({对照组Avg。,卵巢癌组Avg。,的显著特征},…“位置”,“西北”)包含(xAxisLabel);ylabel (yAxisLabel);
这些特性似乎聚集在几个组。你可以进一步调查有多少通过运行以下特性是重要的实验。最常选择的特征用于分类数据,然后使用最频繁的两个选择特性等等,直到所有的特性选择使用10倍以上。您可以看到如果添加更多特性提高了分类器。
nSig = (fCount > 10)之和;nSig cp_rndfeat = 0(20日);为我= 1:nSig为j = 1:20简历=重新分区(简历);P = pca (Y(壮举(1:我)、培训(cv)) ");x = Y(壮举(1:我):)“* P;c =分类(x(测试(cv):), x(培训(cv):), grp(培训(cv)));cp = classperf (grp c测试(cv));cp_rndfeat (j, i) = cp.CorrectRate;%的平均正确分类率结束结束图的阴谋(1:nSig,[马克斯(cp_rndfeat);意味着(cp_rndfeat)]);传奇({“最佳CorrectRate”,“意思是CorrectRate”},“位置”,“东南”)
从这张图可以看到,三功能有时可能得到完美的分类。你还会注意到的最大的意思正确速度发生少量的特性,然后逐渐减少。
[bestAverageCR, bestNumFeatures] = max(平均(cp_rndfeat));
现在您可以可视化功能,给最好的平均分类。你可以看到,这些实际上只有三个峰对应数据。
图;持有在;sigFeats = fCount;sigFeats (sigFeats < = 10) = 0;ax_handle =情节(MZ, [mean_N mean_C]);茎(MZ(壮举(1:bestNumFeatures)), sigFeats(壮举(1:bestNumFeatures)),“r”);轴([7650、8850、80])传说({对照组Avg。,卵巢癌组Avg。,的显著特征})包含(xAxisLabel);ylabel (yAxisLabel);
选择统计学习算法
有许多分类工具MATLAB®,您还可以使用分析蛋白质组学数据。其中包括支持向量机(金宝appfitcsvm),再邻居(fitcknn),神经网络(深度学习工具箱™),和分类树(fitctree)。特征选择,您还可以使用连续的特征子集选择(sequentialfs)或优化利用遗传算法的随机搜索方法(全局优化工具箱)。例如,请参见遗传算法搜索特征质谱数据。
引用
[1]康拉德,T P, V Fusaro,罗斯,D约翰·拉贾帕克萨,B希特,S M斯坦伯格,et al。”高分辨率为卵巢癌血清蛋白质组功能检测。“Endocrine-Related癌症,2004年6月,163 - 78。
[2]Lilien,瑞安·H。,Hany Farid, and Bruce R. Donald. “Probabilistic Disease Classification of Expression-Dependent Proteomic Data from Mass Spectrometry of Human Serum.” Journal of Computational Biology 10, no. 6 (December 2003): 925–46.
[3]李,L。,D. M. Umbach, P. Terry, and J. A. Taylor. “Application of the GA/KNN Method to SELDI Proteomics Data.” Bioinformatics 20, no. 10 (July 1, 2004): 1638–40.
[4]Petricoin,伊曼纽尔F,阿里•M Ardekani本希特,彼得·J·莱文,文森特Fusaro,赛斯M斯坦伯格,戈登•米尔斯,et al。”使用血清蛋白质组学模式识别卵巢癌。“《柳叶刀》359年,没有。9306(2002年2月):572 - 77。
