在Matlab中探索基因表达数据分析.docx
- 文档编号:23148327
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:17
- 大小:130.37KB
在Matlab中探索基因表达数据分析.docx
《在Matlab中探索基因表达数据分析.docx》由会员分享,可在线阅读,更多相关《在Matlab中探索基因表达数据分析.docx(17页珍藏版)》请在冰豆网上搜索。
在Matlab中探索基因表达数据分析
在Matlab中探索基因表达数据分析
日期:
2012-06-25来源:
未知作者:
青岚点击:
547次
摘要:
本文利用Matlab及其生物信息学工具箱提供的函数识别差异表达基因并利用基因本体论确定差异表达基因的生物学功能。
引言包含寡核苷酸或cDNA探针的微阵列可用来比较基因组尺度的基因表达谱,微阵列试验的重要目的在于确定不同条件下,如两种不同的肿瘤类型,是否存在统计显著的基因表达量的变化进而确定差异表达基因的生物学功能。
本文利用一个公共数据集来说明计算
本文利用Matlab及其生物信息学工具箱提供的函数识别差异表达基因并利用基因本体论确定差异表达基因的生物学功能。
引言
包含寡核苷酸或cDNA探针的微阵列可用来比较基因组尺度的基因表达谱,微阵列试验的重要目的在于确定不同条件下,如两种不同的肿瘤类型,是否存在统计显著的基因表达量的变化进而确定差异表达基因的生物学功能。
本文利用一个公共数据集来说明计算过程,这个数据集包括42个胚胎中枢神经系统肿瘤组织(CNS,Pomeroyetal.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,andexpr_cns_gcrma_eb,分别储存用RMA和GCRMA(MLE和EB)预处理的基因表达值。
loadcnsexpressiondata
在每个DataMatrix对象中,每行对应一个HuGeneFl芯片的探针集,每列对应于一个样本,行名是探针集的ID而列名为样本名,本文用expr_cns_gcrma_eb示例,当然也可以用其他对象。
调用get命令获取DataMatrix对象的特征。
get(expr_cns_gcrma_eb)
Name:
'CNSgeneexpressiondata'
RowNames:
{7129x1cell}
ColNames:
{1x42cell}
NRows:
7129
NCols:
42
NDims:
2
ElementClass:
'single'
确定DataMatrix对象expr_cns_gcrma_eb中的基因和样本的数目。
[nGenes,nSamples]=size(expr_cns_gcrma_eb)
nGenes=
7129
nSamples=
42
可以用基因符号来代替探针集的ID用于标记基因表达值,HuGeneFl芯片的基因符号在一个包含Java哈希表的MAT文件中。
loadHuGeneFL_genesymbol_hashtable;
为hu6800genesymbol_hashtable变量创建一个基因表达值的基因符号的cell矩阵。
huGenes=cell(nGenes,1);
fori=1:
nGenes
huGenes{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.NRows
nGenes=
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:
{5758x1cell}
ColNames:
{1x10cell}
NRows:
5758
NCols:
10
NDims:
2
ElementClass:
'single'
MglioData=expr_cns_gcrma_eb(:
Mglios);
get(MglioData)
Name:
''
RowNames:
{5758x1cell}
ColNames:
{1x10cell}
NRows:
5758
NCols:
10
NDims:
2
ElementClass:
'single'
通常t检验是检测两组变量之间显著性差异的标准统计检验,对每个基因执行t检验以识别MD样本和Mglio样本基因表达值的显著性差异,可以通过t得分的正态分位图和t得分及p值的直方图来研究检验结果。
[pvalues,tscores]=mattest(MDData,MglioData,...
'Showhist',true','Showplot',true);
在所有的检验情形下都存在两类误差,当一个非差异表达的基因被标记为差异表达的基因时产生假阳性,而一个差异表达基因未能识别出来时产生假阴性,进行多重假设检验时,即用基因表达数据对上千个基因同时检验员假设时,每个检验都存在其特殊的假阳性率,或加发现率(FDR),FDR定义为假阳性基因数与总阳性数的期望之比(Storeyetal.,2003)。
Storey-Tibshirani例程不仅计算FDR,还得出检验的q值,代表检验的最小FDR,因为FDR的估计依赖于多重检验的真实的空分布,这是未知的,通过交换基因表达数据矩阵的列的替换方法可用来估计真实的空分布(Storeyetal.,2003,Dudoitetal.,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值。
(Dudoitetal.,2003)
pvaluesCorr=mattest(MDData,MglioData,'Permute',10000);
以0.05为p值的阈值确定具有统计显著的差异表达的基因数目,由于替换检验结果不同可能会得到不同的基因数。
cutoff=0.05;
sum(pvaluesCorr ans= 2007 用mafdr对每个检验估计FDR和q值,pi0是研究中真的空假设的总的分数,可以通过bootstrap或立方多项式拟合来估计模拟的空分布,也可以手动设置λ值来估计pi0。 figure; [pFDR,qvalues]=mafdr(pvaluesCorr,'showplot',true); 确定q值低于设定阈值的基因数,同样,由于替换和bootstrap结果的不同可能得出不同的基因数。 sum(qvalues ans= 1925 如果大量基因具有低的FDR表明两组样本具有生物学差异。 通过将mafdr函数的输入参数BHFDR设为真可以调用Benjamini-Hochberg(BH)例程(Benjaminietal,1995)估计FDR调节p值。 pvaluesBH=mafdr(pvaluesCorr,'BHFDR',true); sum(pvaluesBH 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 探索 基因 表达 数据 分析