在Matlab中探索基因表达数据分析

互联网2013-09-06

7073

相关专题: DNA微阵列基因表达数据分析

本文利用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: 7129

NCols: 42

NDims: 2

ElementClass: 'single'

确定DataMatrix对象expr_cns_gcrma_eb中的基因和样本的数目。

[nGenes, nSamples] = size(expr_cns_gcrma_eb)

nGenes =

7129

nSamples =

可以用基因符号来代替探针集的ID用于标记基因表达值，HuGeneFl芯片的基因符号在一个包含Java哈希表的MAT文件中。

load HuGeneFL_genesymbol_hashtable;

为hu6800genesymbol_hashtable变量创建一个基因表达值的基因符号的cell矩阵。

huGenes = cell(nGenes, 1);

for i =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: {5758x1 cell}

ColNames: {1x10 cell}

NRows: 5758

NCols: 10

NDims: 2

ElementClass: 'single'

MglioData = expr_cns_gcrma_eb(:, Mglios);

get(MglioData)

Name: ''

RowNames: {5758x1 cell}

ColNames: {1x10 cell}

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定义为假阳性基因数与总阳性数的期望之比(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-values

PAFAH1B3 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-005

FDR_BH

PAFAH1B3 4.7172e-005

RAB31 8.0657e-005

LRP1 0.00013896

KIAA0367 0.00013896

ID2B 0.00013896

FBL 0.00013896

PLEC1 0.00013896

PEA15 0.00013896

HNRPA1 0.00015953

ID2B 0.00015953

ITGAV 0.00016573

KIR2DL1 0.00016573

SPARCL1 0.00016573

PTOV1 0.00016573

RBMX 0.00016573

H3F3A 0.00016573

DHX9 0.00016573

NAP1L1 0.00016573

C5orf13 0.00016573

ZAP70 0.00016573

MS4A1 0.00016573

GPM6B 0.00016573

RNF113A 0.00016573

一个基因被认定是在两组样本间差异表达的应该具有统计学和生物学意义，本文比较MD和Mglio肿瘤样本之间的基因表达值之比，因此上调基因表明在MD中高表达而下调基因在Mglio中高表达。

用p-值的–log10对生物学效应可以画成火山图，注意，通过火山图用户界面UI，可以交互式地改变p-值的阈值和变化倍数的限定并输出差异表达基因。

diffStruct = mavolcanoplot(MDData, MglioData, pvaluesCorr)

diffStruct =

Name: 'Differentially Expressed'

PVCutoff: 0.0500

FCThreshold: 2

GeneLabels: {116x1 cell}

PValues: [116x1 bioma.data.DataMatrix]

FoldChanges: [116x1 bioma.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 =

确定下调基因数。

nDownGenes = sum(diffStruct.FoldChanges < 0)

nDownGenes =

用基因本体论(GO)注释上调基因

可以通过基因本体论(GO)来注释差异表达基因，从上面的分析结果中查找上调基因，从Gene Ontology Current Annotations下载人类注释(gene_association.goa_human.gz)，解压缩并存入当前目录。

找出上调基因号用于GO分析。

huGenes = rownames(expr_cns_gcrma_eb);

for i = 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}; % Homo sapiens gene list

HGgo = [HGann.GOid]; % associated GO terms

确定分子功能相关的人类注释基因数。

numel(HGgenes)

ans =

74298

并非所有的HuGeneFL芯片上的5758个基因都有注释，通过比较基因符号列表和GO中的基因符号列表可以看出其是否已被注释，可以跟踪注释基因数和每个GO术语关联的上调基因数。

m = GO.Terms(end).id; % gets the last term id

chipgenesCount = zeros(m,1); % a vector of GO term counts for the entire chip.

upgenesCount = zeros(m,1); % a vector of GO term counts for up-regulated genes.

for i = 1:length(huGenes)

idx = strncmpi(HGgenes, huGenes{i}, length(huGenes{i})); % lookup genes

goid = HGgo(idx);

goid = getrelatives(GO,goid);

% Update the tally

chipgenesCount(goid) = chipgenesCount(goid) + 1;

if (any(i == up_geneidx))

upgenesCount(goid) = upgenesCount(goid) +1;

end

可以用超几何分布确定统计显著的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:10

term = idx(i);

report = sprintf('%s%st%-1.5ft%3d / %3dt%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 definition

GO:0030528 0.00044 8 / 489 Plays a role in regulating transcription; may bin...

GO:0003700 0.00090 8 / 538 The function of binding to a specific DNA sequenc...

GO:0003677 0.00165 8 / 583 Interacting selectively with DNA (deoxyribonuclei...

GO:0003705 0.00422 6 / 351 Functions to initiate or regulate RNA polymerase ...

GO:0043169 0.00519 5 / 245 Interacting selectively with cations, charged ato...

GO:0015458 0.00559 1 / 1 ...

GO:0016249 0.00559 1 / 1 Functions to locate the position of ion channels ...

GO:0016974 0.00559 1 / 1 ...

GO:0005509 0.00968 7 / 564 Interacting selectively with calcium ions (Ca2+)....

GO:0030020 0.01069 2 / 30 A constituent of the extracellular matrix that en...

可以研究显著的GO术语并选择与具体的分子功能相关的术语来构建一颗包括其祖先的子本体论树，用biograph函数来可视化显示它，也可为节点着色，本文中红色的节点是最显著的，兰色的节点最不显著，由于人类基因注释文件的频繁更新可能返回不同的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).^(1) pval.^(1/8) pval.^(1/8)];

set(BG.Nodes(i),'Color',color);

end

view(BG)

Gene Ontology object with 8 Terms.

Biograph object with 8 nodes and 9 edges.

从代谢通路中发现差异表达基因

通过KEGG's SOAP Web Service或将基因符号列表发送到KEGG's Web query tool可以从KEGG代谢途径数据库查询差异表达基因的代谢途径信息。

References

[1] Pomeroy, S.L., Tamayo, P., Gaasenbeek, M., Sturla, L.M., Angelo, M., McLaughlin, M.E., Kim, J.Y., Goumnerova, L.C., Black, P.M., Lau, C., Allen, J.C., Zagzag, D., Olson, J.M., Curran, T., Wetmore, C., Biegel, J.A., Poggio, T., Mukherjee, S., Rifkin, R., Califano, A., Stolovitzky, G., Louis, DN, Mesirov, J.P., Lander, E.S., and Golub, T.R. (2002). Prediction of central nervous system embryonal tumour outcome based on gene expression. Nature, 415(6870), 436-442.

[2] Storey, J.D., and Tibshirani, R. (2003). Statistical significance for genomewide studies. Proc.Nat.Acad.Sci., 100(16), 9440-9445.

[3] Dudoit, S., Shaffer, J.P., and Boldrick, J.C. (2003). Multiple hypothesis testing in microarray experiment. Statistical Science, 18, 71-103.

[4] Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. Royal Stat. Soc., B 57, 289-300