在Matlab中探索基因表达数据
- 格式：docx
- 大小：147.57 KB
- 文档页数：11
在Matlab中探索基因表达数据本文利用Matlab及其生物信息学工具箱提供的函数识别差异表达基因并利用基因本体论确定差异表达基因的生物学功能。
引言包含寡核苷酸或cDNA探针的微阵列可用来比较基因组尺度的基因表达谱，微阵列试验的重要目的在于确定不同条件下，如两种不同的肿瘤类型，是否存在统计显著的基因表达量的变化进而确定差异表达基因的生物学功能。
本文利用一个公共数据集来说明计算过程，这个数据集包括42个胚胎中枢神经系统肿瘤组织（CNS, Pomeroy et al. 2002），样本采用Affymetrix 公司出品的HuGeneFL基因芯片进行杂交。
这些CNS数据集（CEL文件）可在CNS实验网站获得，42个肿瘤样本包括10个10 个髓母细胞瘤, 10个横纹肌样脑膜瘤, 10个胶质瘤, 8个幕上原始神经外胚层肿瘤和4个正常人小脑，CNS原始数据集用鲁棒多芯片平均（RMA）和GC鲁棒多芯片平均（GCRMA）进行了预处理。
可以采用t检验和假发现率（FDR）来检测不同肿瘤类型间差异表达的基因，还可以探索与显著上跳基因相关的基因本体论术语。
载入基因表达数据用Load命令加载MAT文件cnsexpressiondata包含三个DataMatrix对象，expr_cns_rma, expr_cns_gcrma_mle, and expr_cns_gcrma_eb，分别储存用RMA和GCRMA（MLE和EB）预处理的基因表达值。
load cnsexpressiondata在每个DataMatrix对象中，每行对应一个HuGeneFl芯片的探针集，每列对应于一个样本，行名是探针集的ID而列名为样本名，本文用expr_cns_gcrma_eb示例，当然也可以用其他对象。
调用get命令获取DataMatrix对象的特征。
get(expr_cns_gcrma_eb)Name: 'CNS gene expression data'RowNames: {7129x1 cell}ColNames: {1x42 cell}NRows: 7129NCols: 42NDims: 2ElementClass: 'single'确定DataMatrix对象expr_cns_gcrma_eb中的基因和样本的数目。
[nGenes, nSamples] = size(expr_cns_gcrma_eb)nGenes =7129nSamples =42可以用基因符号来代替探针集的ID用于标记基因表达值，HuGeneFl芯片的基因符号在一个包含Java哈希表的MAT文件中。
load HuGeneFL_genesymbol_hashtable;为hu6800genesymbol_hashtable变量创建一个基因表达值的基因符号的cell矩阵。
huGenes = cell(nGenes, 1);for i =1:nGeneshuGenes{i} = hu6800genesymbol_hashtable.get(expr_cns_gcrma_eb.RowNames{i});end用DataMatrix的rownames方法将exprs_cns_gcrma_eb中的行名设成基因符号。
expr_cns_gcrma_eb = rownames(expr_cns_gcrma_eb, ':', huGenes);基因表达数据的过滤首先除去没有基因符号的表达数据，如标成'---'的空符号。
expr_cns_gcrma_eb('---', :) = [];在这个研究中很多基因没有表达或在样本间变化很小，这些基因需要用非特异性过滤除去。
用genelowvalfilter函数滤除绝对表达量值很低的基因。
[mask, expr_cns_gcrma_eb] = genelowvalfilter(expr_cns_gcrma_eb);用genevarfilter函数滤除样本间方差很小的基因。
[mask, expr_cns_gcrma_eb] = genevarfilter(expr_cns_gcrma_eb);确定过滤以后的基因数目。
nGenes = expr_cns_gcrma_eb.NRowsnGenes =5758识别差异基因表达现在可以比较一下CNS髓母细胞瘤(MD)和非神经源恶性胶质瘤(Mglio)之间基因表达值的差异了。
从42个样本中提取10个MD和10个Mglio样本数据。
MDs = strncmp(expr_cns_gcrma_eb.ColNames,'Brain_MD', 8);Mglios = strncmp(expr_cns_gcrma_eb.ColNames,'Brain_MGlio', 11);MDData = expr_cns_gcrma_eb(:, MDs);get(MDData)Name: ''RowNames: {5758x1 cell}ColNames: {1x10 cell}NRows: 5758NCols: 10NDims: 2ElementClass: 'single'MglioData = expr_cns_gcrma_eb(:, Mglios);get(MglioData)Name: ''RowNames: {5758x1 cell}ColNames: {1x10 cell}NRows: 5758NCols: 10NDims: 2ElementClass: 'single'通常t检验是检测两组变量之间显著性差异的标准统计检验，对每个基因执行t检验以识别MD样本和Mglio样本基因表达值的显著性差异，可以通过t得分的正态分位图和t得分及p 值的直方图来研究检验结果。
[pvalues, tscores] = mattest(MDData, MglioData,...'Showhist', true', 'Showplot', true);在所有的检验情形下都存在两类误差，当一个非差异表达的基因被标记为差异表达的基因时产生假阳性，而一个差异表达基因未能识别出来时产生假阴性，进行多重假设检验时，即用基因表达数据对上千个基因同时检验员假设时，每个检验都存在其特殊的假阳性率，或加发现率（FDR），FDR定义为假阳性基因数与总阳性数的期望之比（Storey et al., 2003）。
Storey-Tibshirani例程不仅计算FDR,还得出检验的q值，代表检验的最小FDR，因为FDR的估计依赖于多重检验的真实的空分布，这是未知的，通过交换基因表达数据矩阵的列的替换方法可用来估计真实的空分布（Storey et al., 2003, Dudoit et al., 2003），鉴于样本数量，考虑所有的替换是不可行的通常在大样本下只考虑一个随机子集，本文利用统计工具箱中的nchoosek函数找出所有可能的替换数。
all_possible_perms = nchoosek(1:MDData.NCols+MglioData.NCols, MDData.NCols);size(all_possible_perms, 1)ans =184756调用mattest的PERMUTE选项进行替换的t检验，通过对基因表达数据矩阵MDData和MglioData的列进行10，000次替换来算出p值。
（Dudoit et al., 2003）pvaluesCorr = mattest(MDData, MglioData, 'Permute', 10000);以0.05为p值的阈值确定具有统计显著的差异表达的基因数目，由于替换检验结果不同可能会得到不同的基因数。
cutoff = 0.05;sum(pvaluesCorr < cutoff)ans =2007用mafdr对每个检验估计FDR和q值，pi0是研究中真的空假设的总的分数，可以通过bootstrap或立方多项式拟合来估计模拟的空分布，也可以手动设置λ值来估计pi0。
figure;[pFDR, qvalues] = mafdr(pvaluesCorr, 'showplot', true);确定q值低于设定阈值的基因数，同样，由于替换和bootstrap结果的不同可能得出不同的基因数。
sum(qvalues < cutoff)ans =1925如果大量基因具有低的FDR表明两组样本具有生物学差异。
通过将mafdr函数的输入参数BHFDR设为真可以调用Benjamini-Hochberg (BH)例程（Benjamini et al, 1995）估计FDR调节p值。
pvaluesBH = mafdr(pvaluesCorr, 'BHFDR', true);sum(pvaluesBH < cutoff)ans =1275可以将检验结果包括t得分, p值, pFDRs, q值和BH FDR校正p-values存为DataMatrix对象。
testResults = [tscores pvaluesCorr pFDR qvalues pvaluesBH];可以用DataMatrix对象的colnames方法将列名更新为BH FDR校正的p值。
testResults = colnames(testResults, 5, {'FDR_BH'});可以用sortrows方法对p-值pvaluesCorr进行排序。
testResults = sortrows(testResults, 2);下面显示检验结果的前23个基因，注意，由于替换测验和bootstrap 结果的差异，最终结果可能和这里显示的有所不同。
testResults(1:23, :)ans =t-scores p-values FDR q-valuesPAFAH1B3 10.211 8.1924e-009 1.8315e-005 1.8315e-005 RAB31 -13.702 2.8016e-008 3.1316e-005 3.1316e-005 LRP1 -9.0742 1.6453e-007 0.00012261 5.3953e-005 KIAA0367 -8.9398 1.7404e-007 9.7272e-005 5.3953e-005 ID2B -8.8966 1.7782e-007 7.9509e-005 5.3953e-005 FBL 8.8668 1.8071e-007 6.7333e-005 5.3953e-005 PLEC1 -8.7961 1.8858e-007 6.0228e-005 5.3953e-005 PEA15 -8.7637 1.9306e-007 5.3953e-005 5.3953e-005 HNRPA1 8.4532 2.6546e-007 6.5942e-005 6.1939e-005 ID2B -8.4194 2.7705e-007 6.1939e-005 6.1939e-005 ITGAV -8.265 3.5525e-007 7.2201e-005 6.4346e-005 KIR2DL1 -8.1355 3.8039e-007 7.0867e-005 6.4346e-005 SPARCL1 -8.1246 3.8603e-007 6.6387e-005 6.4346e-005 PTOV1 7.9738 5.5184e-007 8.8123e-005 6.4346e-005 RBMX 7.9186 6.0384e-007 8.9997e-005 6.4346e-005 H3F3A 7.8816 6.2271e-007 8.7009e-005 6.4346e-005 DHX9 7.8439 6.6235e-007 8.7103e-005 6.4346e-005 NAP1L1 7.8241 6.6698e-007 8.2839e-005 6.4346e-005 C5orf13 7.8037 6.6868e-007 7.868e-005 6.4346e-005 ZAP70 7.7977 6.6897e-007 7.4778e-005 6.4346e-005 MS4A1 -7.7607 6.7364e-007 7.1715e-005 6.4346e-005 GPM6B -7.7605 6.737e-007 6.8461e-005 6.4346e-005 RNF113A -7.7548 6.7548e-007 6.5658e-005 6.4346e-005FDR_BHPAFAH1B3 4.7172e-005RAB31 8.0657e-005LRP1 0.00013896KIAA0367 0.00013896ID2B 0.00013896FBL 0.00013896PLEC1 0.00013896PEA15 0.00013896HNRPA1 0.00015953ID2B 0.00015953ITGAV 0.00016573KIR2DL1 0.00016573SPARCL1 0.00016573PTOV1 0.00016573RBMX 0.00016573H3F3A 0.00016573DHX9 0.00016573NAP1L1 0.00016573C5orf13 0.00016573ZAP70 0.00016573MS4A1 0.00016573GPM6B 0.00016573RNF113A 0.00016573一个基因被认定是在两组样本间差异表达的应该具有统计学和生物学意义，本文比较MD和Mglio肿瘤样本之间的基因表达值之比，因此上调基因表明在MD中高表达而下调基因在Mglio中高表达。
