在Matlab中探索基因表达数据

1、在Matlab中探索基因表达数据本文利用Matlab及其生物信息学工具箱提供的函数识别差异表达基因并利用基因本体论确 定差异表达基因的生物学功能。引言包含寡核苷酸或cDNA探针的微阵列可用来比较基因组尺度的基因表达谱，微阵列试验的重 要目的在于确定不同条件下，如两种不同的肿瘤类型，是否存在统计显著的基因表达量的变 化进而确定差异表达基因的生物学功能。本文利用一个公共数据集来说明计算过程，这个数据集包括42个胚胎中枢神经系统肿瘤组 织（CNS, Pomeroy et al. 2002），样本采用Affymetrix公司出品的HuGeneFL基因芯片进行杂 交。这些CNS数据集（CEL文件）可在C

2、NS实验网站获得，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，分别储

3、存用 RMA 和 GCRMA（MLE 和 EB）预 处理的基因表达值。load cnsexpressiondata在每个DataMatrix对象中，每行对应一个HuGeneFl芯片的探针集，每列对应于一个样本， 行名是探针集的ID而列名为样本名，本文用expr_cns_gcrma_eb示例，当然也可以用其他对 象。调用get命令获取DataMatrix对象的特征。get（expr_cns_gcrma_eb）Name: CNS gene expression dataRowNames: 7129x1 cellColNames: 1x42 cellNRows: 7129NCols: 42NDims

4、: 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

5、, 1);for i =1:nGeneshuGenesi = hu6800genesymbol_hashtable.get(expr_cns_gcrma_eb.RowNamesi);end用DataMatrix的rownames方法将exprs_cns_gcrma_eb中的行名设成基因符号。expr_cns_gcrma_eb = rownames(expr_cns_gcrma_eb, :, huGenes);基因表达数据的过滤首先除去没有基因符号的表达数据，如标成-的空符号。expr_cns_gcrma_eb(-,:)=;在这个研究中很多基因没有表达或在样本间变化很小，这些基因需要用非特异性过

6、滤除去。用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)之间基因表达值的差 异了

7、。从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 cellColNames: 1x10 cellNRows: 5758NCols: 10NDims: 2ElementClass: singleMglioData = expr

8、_cns_gcrma_eb(:, Mglios);get(MglioData)Name: RowNames: 5758x1 cellColNames: 1x10 cellNRows: 5758NCols: 10NDims: 2ElementClass: single通常t检验是检测两组变量之间显著性差异的标准统计检验，对每个基因执行t检验以识别 MD样本和Mglio样本基因表达值的显著性差异，可以通过t得分的正态分位图和t得分及p 值的直方图来研究检验结果。pvalues, tscores = mattest(MDData, MglioData,.Showhist, true, Showplo

9、t, true);Histograms of t-test Results600500400t-scores200 -100 -Q I-10-50t-score在所有的检验情形下都存在两类误差，当一个非差异表达的基因被标记为差异表达的基因时 产生假阳性，而一个差异表达基因未能识别出来时产生假阴性，进行多重假设检验时，即用 基因表达数据对上千个基因同时检验员假设时，每个检验都存在其特殊的假阳性率，或加发 现率(FDR)，FDR定义为假阳性基因数与总阳性数的期望之比(Storey et al., 2003)。Storey-Tibshirani例程不仅计算FDR,还得出检验的q值，代表检验的最小FD

10、R，因为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

11、的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

12、或立方多项式拟合来估计模拟的空分布，也可以手动设置久值来估计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 =

13、mafdr(pvaluesCorr, BHFDR, true);sum(pvaluesBH &_bolJFEL3.9147Q54C-0KIPTQV1123A-D07iPTMA33 95M+6e-M7H3F3AC5drf133.7 5471T2C-D07漩E1.50y2JJB.LHKRPS111.51069555-006CCNG11.5246S72B-DO6SFRS3旺W006NSktld3.717061 9b-0D6二JUsSeaui-rtecGents.amIuwREA 5jl3.4733739E-008KliADW*4-3 71 M83e-D08D2BT.7 092391 e-008D26

14、I .35S65576-007SPARCL13.4316961B-D07却响异岩白昵-曲？PTFRZ167563213C-007HTR1/.67S8526&00?PCDFKC31 fO72ie-DWD0R12.243M25B-DD6RABBI誓 2&的丁 5dDO6Down Regul 朋eel阡 wk国s Fold changeCtesj改阳rt./ Volcano Plot按下Ctrl键的同时点击基因列表中的基因名可以在图中找到基因，如火山图所见，小脑颗粒 神经元细胞特异的基因ZIC和NEUROD上调，而星形和少突胶质细胞分化基因如SOX2, PEA15, 和ID2B,则下调。确定差异表达

15、基因数。nDiffGenes = diffStruct.PValues.NRowsnDiffGenes =116获得上调基因列表。up_geneidx = find(diffStruct.FoldChanges 0);up_genes = rownames(diffStruct.FoldChanges, up_geneidx);nUpGenes = length(up_geneidx)nUpGenes =75确定下调基因数。nDownGenes = sum(diffStruct.FoldChanges 0)nDownGenes =41用基因本体论(GO)注释上调基因可以通过基因本体论(GO)来

16、注释差异表达基因，从上面的分析结果中查找上调基因，从 Gene Ontology Current Annotations下载人类注释(gene association.goa human.gz)，解压缩 并存入当前目录。找出上调基因号用于GO分析。huGenes = rownames(expr_cns_gcrma_eb);for i = 1:nUpGenesup_geneidx(i) = find(strncmpi(huGenes, up_genesi, length(up_genesi), 1);end用geneont函数加载GO数据库为MATLAB对象。GO = geneont(live,

17、true);读人类注释文件，本文只看与分子功能相关的基因，故只需读Aspect域设为“F”信息，基 因符号和对应ID是我们感兴趣的域，在GO注释文件中分别为DB_Object_Symbol和GOid 域。HGann = goannotread(gene_association.goa_human,. Aspect,F,Fields,DB_Object_Symbol,GOid);创建人类基因及相应GO术语的列表。HGgenes = HGann.DB_Object_Symbol; % Homo sapiens gene listHGgo = HGann.GOid; % associated GO

18、terms确定分子功能相关的人类注释基因数。numel(HGgenes)ans =74298并非所有的HuGeneFL芯片上的5758个基因都有注释，通过比较基因符号列表和GO中的基 因符号列表可以看出其是否已被注释，可以跟踪注释基因数和每个GO术语关联的上调基因 数。m = GO.Terms(end).id; % gets the last term idchipgenesCount = zeros(m,1); % a vector of GO term counts for the entire chip.upgenesCount = zeros(m,1); % a vector of G

19、O term counts for up-regulated genes.for i = 1:length(huGenes)idx = strncmpi(HGgenes, huGenesi, length(huGenesi); % lookup genesgoid = HGgo(idx);goid = getrelatives(GO,goid);% Update the tallychipgenesCount(goid) = chipgenesCount(goid) + 1;if (any(i = up_geneidx)upgenesCount(goid) = upgenesCount(goi

20、d) +1;endend可以用超几何分布确定统计显著的GO术语，对于每一个GO术语，p-值表示关联的注释基 因由于偶然性获得的概率。gopvalues = hygepdf(upgenesCount,max(chipgenesCount),.max(upgenesCount),chipgenesCount);dummy, idx = sort(gopvalues);report = sprintf(GO Term p-value counts definitionn); for i = 1:10term = idx(i);report = sprintf(%s%st%-1.5ft%3d / %3

21、dt%s.n,.report, char(num2goid(term), gopvalues(term),. upgenesCount(term), chipgenesCount(term),.GO(term).Term.definition(2:min(50,end);end disp(report);GO Term p-value counts definitionGO:00305280.000448 / 489GO:00037000.000908 / 538GO:00036770.001658 / 583GO:00037050.004226 / 351GO:00431690.005195

22、 / 245GO:00154580.005591 / 1GO:00162490.005591 / 1GO:00169740.005591 / 1GO:00055090.009687 / 564GO:00300200.010692 / 30Plays a role in regulating transcription; may bin. The function of binding to a specific DNA sequenc. Interacting selectively with DNA (deoxyribonuclei. Functions to initiate or reg

23、ulate RNA polymerase .Interacting selectively with cations, charged ato.Functions to locate the position of ion channels .Interacting selectively with calcium ions (Ca2+).A constituent of the extracellular matrix that en.可以研究显著的GO术语并选择与具体的分子功能相关的术语来构建一颗包括其祖先的子本 体论树，用biograph函数来可视化显示它，也可为节点着色，本文中红色的节

24、点是最显著 的，兰色的节点最不显著，由于人类基因注释文件的频繁更新可能返回不同的GO术语。fcnAncestors = GO(getancestors(GO,idx(1:4)cm acc rels = getmatrix(fcnAncestors);BG = biograph(cm,get(fcnAncestors.Terms,name)for i=1:numel(acc)pval = gopvalues(acc(i);color = (1-pval).A(1) pval.A(1/8) pval.A(1/8);set(BG.Nodes(i),Color,color);endview(BG)Ge

25、ne Ontology object with 8 Terms.Biograph object with 8 nodes and 9 edges.bindingmolecular functionRNA palymerase II transcription factor activitytranscription regulator activitytranscription factor activitynucleic acid bindingDNA binding从代谢通路中发现差异表达基因通过KEGGs SOAP Web Service或将基因符号列表发送到KEGGs Web query tool可以从KEGG 代谢途径数据库查询差异表达基因的代谢途径信息。ReferencesPomeroy, S.L., Tamayo, P Gaasenbeek, M., Sturla, L.M.,