功能。

引言

包含寡核苷酸或cDNA探針的微陣列可用來比較基因組尺度的基因表達(dá)譜，微陣列試驗(yàn)的重要目的在于確定不同條件下，如兩種不同的腫瘤類型，是否存在統(tǒng)計(jì)顯著的基因表達(dá)量的變化進(jìn)而確定差異表達(dá)基因的生物學(xué)功能。

本文利用一個(gè)公共數(shù)據(jù)集來說明計(jì)算過程，這個(gè)數(shù)據(jù)集包括42個(gè)胚胎中樞神經(jīng)系統(tǒng)腫瘤組織(CNS, Pomeroy et al. 2002)，樣本采用Affymetrix 公司出品的HuGeneFL基因芯片進(jìn)行雜交。

這些CNS數(shù)據(jù)集(CEL文件)可在CNS實(shí)驗(yàn)獲得，42個(gè)腫瘤樣本包括10個(gè)10 個(gè)髓母細(xì)胞瘤, 10個(gè)橫紋肌樣腦膜瘤, 10個(gè)膠質(zhì)瘤, 8個(gè)幕上原始神經(jīng)外胚層腫瘤和4個(gè)正常人小腦，CNS原始數(shù)據(jù)集用魯棒多芯片平均(RMA)和GC魯棒多芯片平均(GCRMA)進(jìn)行了預(yù)處理。

可以采用t檢驗(yàn)和假發(fā)現(xiàn)率(FDR)來檢測(cè)不同腫瘤類型間差異表達(dá)的基因，還可以探索與顯著上跳基因相關(guān)的基因本體論術(shù)語(yǔ)。

載入基因表達(dá)數(shù)據(jù)

用Load命令加載MAT文件cnsexpressiondata包含三個(gè)DataMatrix對(duì)象，expr_cns_rma, expr_cns_gcrma_mle, and expr_cns_gcrma_eb，分別儲(chǔ)存用RMA和GCRMA(MLE和EB)預(yù)處理的基因表達(dá)值。

load cnsexpressiondata

在每個(gè)DataMatrix對(duì)象中，每行對(duì)應(yīng)一個(gè)HuGeneFl芯片的探針集，每列對(duì)應(yīng)于一個(gè)樣本，行名是探針集的ID而列名為樣本名，本文用expr_cns_gcrma_eb示例，當(dāng)然也可以用其他對(duì)象。

調(diào)用get命令獲取DataMatrix對(duì)象的特征。

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對(duì)象expr_cns_gcrma_eb中的基因和樣本的數(shù)目。

[nGenes, nSamples] = size(expr_cns_gcrma_eb)

nGenes =

7129

nSamples =

42

可以用基因符號(hào)來代替探針集的ID用于標(biāo)記基因表達(dá)值，HuGeneFl芯片的基因符號(hào)在一個(gè)包含Java哈希表的MAT文件中。

load HuGeneFL_genesymbol_hashtable;

為hu6800genesymbol_hashtable變量創(chuàng)建一個(gè)基因表達(dá)值的基因符號(hào)的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中的行名設(shè)成基因符號(hào)。

expr_cns_gcrma_eb = rownames(expr_cns_gcrma_eb, ':', huGenes);

基因表達(dá)數(shù)據(jù)的過濾

首先除去沒有基因符號(hào)的表達(dá)數(shù)據(jù)，如標(biāo)成'---'的空符號(hào)。

expr_cns_gcrma_eb('---', :) = [];

在這個(gè)研究中很多基因沒有表達(dá)或在樣本間變化很小，這些基因需要用非特異性過濾除去。

用genelowvalfilter函數(shù)濾除表達(dá)量值很低的基因。

[mask, expr_cns_gcrma_eb] = genelowvalfilter(expr_cns_gcrma_eb);

用genevarfilter函數(shù)濾除樣本間方差很小的基因。

[mask, expr_cns_gcrma_eb] = genevarfilter(expr_cns_gcrma_eb);

確定過濾以后的基因數(shù)目。

nGenes = expr_cns_gcrma_eb.NRows

nGenes =

5758

識(shí)別差異基因表達(dá)

現(xiàn)在可以比較一下CNS髓母細(xì)胞瘤(MD)和非神經(jīng)源惡性膠質(zhì)瘤(Mglio)之間基因表達(dá)值的差異了。

從42個(gè)樣本中提取10個(gè)MD和10個(gè)Mglio樣本數(shù)據(jù)。

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檢驗(yàn)是檢測(cè)兩組變量之間顯著性差異的標(biāo)準(zhǔn)統(tǒng)計(jì)檢驗(yàn)，對(duì)每個(gè)基因執(zhí)行t檢驗(yàn)以識(shí)別MD樣本和Mglio樣本基因表達(dá)值的顯著性差異，可以通過t得分的正態(tài)分位圖和t得分及p值的直方圖來研究檢驗(yàn)結(jié)果。

[pvalues, tscores] = mattest(MDData, MglioData,...

'Showhist', true', 'Showplot', true);

在所有的檢驗(yàn)情形下都存在兩類誤差，當(dāng)一個(gè)非差異表達(dá)的基因被標(biāo)記為差異表達(dá)的基因時(shí)產(chǎn)生假陽(yáng)性，而一個(gè)差異表達(dá)基因未能識(shí)別出來時(shí)產(chǎn)生假陰性，進(jìn)行多重假設(shè)檢驗(yàn)時(shí)，即用基因表達(dá)數(shù)據(jù)對(duì)上千個(gè)基因同時(shí)檢驗(yàn)員假設(shè)時(shí)，每個(gè)檢驗(yàn)都存在其特殊的假陽(yáng)性率，或加發(fā)現(xiàn)率(FDR)，F(xiàn)DR定義為假陽(yáng)性基因數(shù)與總陽(yáng)性數(shù)的期望之比(Storey et al., 2003)。Storey-Tibshirani例程不僅計(jì)算FDR,還得出檢驗(yàn)的q值，代表檢驗(yàn)的zui小FDR，因?yàn)镕DR的估計(jì)依賴于多重檢驗(yàn)的真實(shí)的空分布，這是未知的，通過交換基因表達(dá)數(shù)據(jù)矩陣的列的替換方法可用來估計(jì)真實(shí)的空分布(Storey et al., 2003, Dudoit et al., 2003)，鑒于樣本數(shù)量，考慮所有的替換是不可行的通常在大樣本下只考慮一個(gè)隨機(jī)子集，本文利用統(tǒng)計(jì)工具箱中的nchoosek函數(shù)找出所有可能的替換數(shù)。

all_possible_perms = nchoosek(1:MDData.NCols+MglioData.NCols, MDData.NCols);

size(all_possible_perms, 1)

ans =

184756

調(diào)用mattest的PERMUTE選項(xiàng)進(jìn)行替換的t檢驗(yàn)，通過對(duì)基因表達(dá)數(shù)據(jù)矩陣MDData和MglioData的列進(jìn)行10，000次替換來算出p值。(Dudoit et al., 2003)

pvaluesCorr = mattest(MDData, MglioData, 'Permute', 10000);

以0.05為p值的閾值確定具有統(tǒng)計(jì)顯著的差異表達(dá)的基因數(shù)目，由于替換檢驗(yàn)結(jié)果不同可能會(huì)得到不同的基因數(shù)。

cutoff = 0.05;

sum(pvaluesCorr < cutoff)

ans =

2007

用mafdr對(duì)每個(gè)檢驗(yàn)估計(jì)FDR和q值，pi0 是研究中真的空假設(shè)的總的分?jǐn)?shù)，可以通過bootstrap或立方多項(xiàng)式擬合來估計(jì)模擬的空分布，也可以手動(dòng)設(shè)置λ值來估計(jì)pi0 。

figure;

[pFDR, qvalues] = mafdr(pvaluesCorr, 'showplot', true);

按下Ctrl鍵的同時(shí)點(diǎn)擊基因列表中的基因名可以在圖中找到基因，如火山圖所見，小腦顆粒神經(jīng)元細(xì)胞特異的基因ZIC 和NEUROD 上調(diào)，而星形和少突膠質(zhì)細(xì)胞分化基因如SOX2, PEA15, 和ID2B,則下調(diào)。

確定差異表達(dá)基因數(shù)。

nDiffGenes = diffStruct.PValues.NRows

nDiffGenes =

116

獲得上調(diào)基因列表。

up_geneidx = find(diffStruct.FoldChanges > 0);

up_genes = rownames(diffStruct.FoldChanges, up_geneidx);

nUpGenes = length(up_geneidx)

nUpGenes =

75

確定下調(diào)基因數(shù)。

nDownGenes = sum(diffStruct.FoldChanges < 0)

nDownGenes =

41

用基因本體論(GO)注釋上調(diào)基因

可以通過基因本體論(GO)來注釋差異表達(dá)基因，從上面的分析結(jié)果中查找上調(diào)基因，從Gene Ontology Current Annotations下載人類注釋(gene_association.goa_human.gz)，解壓縮并存入當(dāng)前目錄。

找出上調(diào)基因號(hào)用于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函數(shù)加載GO數(shù)據(jù)庫(kù)為MATLAB對(duì)象。

GO = geneont('live',true);

讀人類注釋文件，本文只看與分子功能相關(guān)的基因，故只需讀Aspect域設(shè)為“F”的信息，基因符號(hào)和對(duì)應(yīng)ID是我們感興趣的域，在GO注釋文件中分別為DB_Object_Symbol和GOid域。

HGann = goannotread('gene_association.goa_human',...

'Aspect','F','Fields',{'DB_Object_Symbol','GOid'});

創(chuàng)建人類基因及相應(yīng)GO術(shù)語(yǔ)的列表。

HGgenes = {HGann.DB_Object_Symbol}; % Homo sapiens gene list

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

確定分子功能相關(guān)的人類注釋基因數(shù)。

numel(HGgenes)

ans =

74298

并非所有的HuGeneFL芯片上的5758個(gè)基因都有注釋，通過比較基因符號(hào)列表和GO中的基因符號(hào)列表可以看出其是否已被注釋，可以跟蹤注釋基因數(shù)和每個(gè)GO術(shù)語(yǔ)關(guān)聯(lián)的上調(diào)基因數(shù)。

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

end

可以用超幾何分布確定統(tǒng)計(jì)顯著的GO術(shù)語(yǔ)，對(duì)于每一個(gè)GO術(shù)語(yǔ)，p-值表示關(guān)聯(lián)的注釋基因由于偶然性獲得的概率。

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術(shù)語(yǔ)并選擇與具體的分子功能相關(guān)的術(shù)語(yǔ)來構(gòu)建一顆包括其祖先的子本體論樹，用biograph函數(shù)來可視化顯示它，也可為節(jié)點(diǎn)著色，本文中紅色的節(jié)點(diǎn)是zui顯著的，蘭色的節(jié)點(diǎn)zui不顯著，由于人類基因注釋文件的頻繁更新可能返回不同的GO術(shù)語(yǔ)。

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.

<img alt="在Matlab中探索基因表達(dá)數(shù)據(jù)分析" 在matlab中探索基因表達(dá)數(shù)據(jù)分析"="" border="1" height="420" data-cke-saved-src="http://www.bio1000.com/uploads/allimg/120625/1615504c2-4.jpg" src="http://www.bio1000.com/uploads/allimg/120625/1615504c2-4.jpg" width="560" style="vertical-align: middle; border: 0px;">

從代謝通路中發(fā)現(xiàn)差異表達(dá)基因

通過KEGG's SOAP Web Service或?qū)⒒蚍?hào)列表發(fā)送到KEGG's Web query tool可以從KEGG代謝途徑數(shù)據(jù)庫(kù)查詢差異表達(dá)基因的代謝途徑信息。

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.