ans=
1275
可以将检验结果包括t得分,p值,pFDRs,q值和BHFDR校正p-values存为DataMatrix对象。
testResults=[tscorespvaluesCorrpFDRqvaluespvaluesBH];
可以用DataMatrix对象的colnames方法将列名更新为BHFDR校正的p值。
testResults=colnames（testResults,5,{'FDR_BH'}）;
可以用sortrows方法对p-值pvaluesCorr进行排序。
testResults=sortrows（testResults,2）;
下面显示检验结果的前23个基因，注意，由于替换测验和bootstrap结果的差异，最终结果可能和这里显示的有所不同。
testResults（1:
23,:
）
ans=
t-scoresp-valuesFDRq-values
PAFAH1B310.2118.1924e-0091.8315e-0051.8315e-005
RAB31-13.7022.8016e-0083.1316e-0053.1316e-005
LRP1-9.07421.6453e-0070.000122615.3953e-005
KIAA0367-8.93981.7404e-0079.7272e-0055.3953e-005
ID2B-8.89661.7782e-0077.9509e-0055.3953e-005
FBL8.86681.8071e-0076.7333e-0055.3953e-005
PLEC1-8.79611.8858e-0076.0228e-0055.3953e-005
PEA15-8.76371.9306e-0075.3953e-0055.3953e-005
HNRPA18.45322.6546e-0076.5942e-0056.1939e-005
ID2B-8.41942.7705e-0076.1939e-0056.1939e-005
ITGAV-8.2653.5525e-0077.2201e-0056.4346e-005
KIR2DL1-8.13553.8039e-0077.0867e-0056.4346e-005
SPARCL1-8.12463.8603e-0076.6387e-0056.4346e-005
PTOV17.97385.5184e-0078.8123e-0056.4346e-005
RBMX7.91866.0384e-0078.9997e-0056.4346e-005
H3F3A7.88166.2271e-0078.7009e-0056.4346e-005
DHX97.84396.6235e-0078.7103e-0056.4346e-005
NAP1L17.82416.6698e-0078.2839e-0056.4346e-005
C5orf137.80376.6868e-0077.868e-0056.4346e-005
ZAP707.79776.6897e-0077.4778e-0056.4346e-005
MS4A1-7.76076.7364e-0077.1715e-0056.4346e-005
GPM6B-7.76056.737e-0076.8461e-0056.4346e-005
RNF113A-7.75486.7548e-0076.5658e-0056.4346e-005
FDR_BH
PAFAH1B34.7172e-005
RAB318.0657e-005
LRP10.00013896
KIAA03670.00013896
ID2B0.00013896
FBL0.00013896
PLEC10.00013896
PEA150.00013896
HNRPA10.00015953
ID2B0.00015953
ITGAV0.00016573
KIR2DL10.00016573
SPARCL10.00016573
PTOV10.00016573
RBMX0.00016573
H3F3A0.00016573
DHX90.00016573
NAP1L10.00016573
C5orf130.00016573
ZAP700.00016573
MS4A10.00016573
GPM6B0.00016573
RNF113A0.00016573
一个基因被认定是在两组样本间差异表达的应该具有统计学和生物学意义，本文比较MD和Mglio肿瘤样本之间的基因表达值之比，因此上调基因表明在MD中高表达而下调基因在Mglio中高表达。
用p-值的–log10对生物学效应可以画成火山图，注意，通过火山图用户界面UI，可以交互式地改变p-值的阈值和变化倍数的限定并输出差异表达基因。
diffStruct=mavolcanoplot（MDData,MglioData,pvaluesCorr）
diffStruct=
Name:
'DifferentiallyExpressed'
PVCutoff:
0.0500
FCThreshold:
2
GeneLabels:
{116x1cell}
PValues:
[116x1bioma.data.DataMatrix]
FoldChanges:
[116x1bioma.data.DataMatrix]
按下Ctrl键的同时点击基因列表中的基因名可以在图中找到基因，如火山图所见，小脑颗粒神经元细胞特异的基因ZIC和NEUROD上调，而星形和少突胶质细胞分化基因如SOX2,PEA15,和ID2B,则下调。
确定差异表达基因数。
nDiffGenes=diffStruct.PValues.NRows
nDiffGenes=
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）来注释差异表达基因，从上面的分析结果中查找上调基因，从GeneOntologyCurrentAnnotations下载人类注释（gene_association.goa_human.gz），解压缩并存入当前目录。
找出上调基因号用于GO分析。
huGenes=rownames（expr_cns_gcrma_eb）;
fori=1:
nUpGenes
up_geneidx（i）=find（strncmpi（huGenes,up_genes{i},length（up_genes{i}））,1）;
end
用geneont函数加载GO数据库为MATLAB对象。
GO=geneont（'live',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};%Homosapiensgenelist
HGgo=[HGann.GOid];%associatedGOterms
确定分子功能相关的人类注释基因数。
numel（HGgenes）
ans=
74298
并非所有的HuGeneFL芯片上的5758个基因都有注释，通过比较基因符号列表和GO中的基因符号列表可以看出其是否已被注释，可以跟踪注释基因数和每个GO术语关联的上调基因数。
m=GO.Terms（end）.id;%getsthelasttermid
chipgenesCount=zeros（m,1）;%avectorofGOtermcountsfortheentirechip.
upgenesCount=zeros（m,1）;%avectorofGOtermcountsforup-regulatedgenes.
fori=1:
length（huGenes）
idx=strncmpi（HGgenes,huGenes{i},length（huGenes{i}））;%lookupgenes
goid=HGgo（idx）;
goid=getrelatives（GO,goid）;
%Updatethetally
chipgenesCount（goid）=chipgenesCount（goid）+1;
if（any（i==up_geneidx））
upgenesCount（goid）=upgenesCount（goid）+1;
end
end
可以用超几何分布确定统计显著的GO术语，对于每一个GO术语，p-值表示关联的注释基因由于偶然性获得的概率。
gopvalues=hygepdf（upgenesCount,max（chipgenesCount）,...
max（upgenesCount）,chipgenesCount）;
[dummy,idx]=sort（gopvalues）;
report=sprintf（'GOTermp-valuecountsdefinition\n'）;
fori=1:
10
term=idx（i）;
report=sprintf（'%s%s\t%-1.5f\t%3d/%3d\t%s...\n',...
report,char（num2goid（term））,gopvalues（term）,...
upgenesCount（term）,chipgenesCount（term）,...
GO（term）.Term.definition（2:
min（50,end）））;
end
disp（report）;
GOTermp-valuecountsdefinition
GO:
00305280.000448/489Playsaroleinregulatingtranscription;maybin...
GO:
00037000.000908/538ThefunctionofbindingtoaspecificDNAsequenc...
GO:
00036770.001658/583InteractingselectivelywithDNA（deoxyribonuclei...
GO:
00037050.004226/351FunctionstoinitiateorregulateRNApolymerase...
GO:
00431690.005195/245Interactingselectivelywithcations,chargedato...
GO:
00154580.005591/1...
GO:
00162490.005591/1Functionstolocatethepositionofionchannels...
GO:
00169740.005591/1...
GO:
00055090.009687/564Interactingselectivelywithcalciumions（Ca2+）....
GO:
00300200.010692/30Aconstituentoftheextracellularmatrixthaten...
可以研究显著的GO术语并选择与具体的分子功能相关的术语来构建一颗包括其祖先的子本体论树，用biograph函数来可视化显示它，也可为节点着色，本文中红色的节点是最显著的，兰色的节点最不显著，由于人类基因注释文件的频繁更新可能返回不同的GO术语。
fcnAncestors=GO（getancestors（GO,idx（1:
4）））
[cmaccrels]=getmatrix（fcnAncestors）;
BG=biograph（cm,get（fcnAncestors.Terms,'name'））
fori=1:
numel（acc）
pval=gopvalues（acc（i））;
color=[（1-pval）.^
（1）pval.^（1/8）pval.^（1/8）];
set（BG.Nodes（i）,'Color',color）;
end
view（BG）
GeneOntologyobjectwith8Terms.
Biographobjectwith8nodesand9edges.
从代谢通路中发现差异表达基因
通过K
升级会员 