功能。

引言

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

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

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

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

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

用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示例，當然也可以用其他對象。

調(diào)用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中的基因和樣本的數(shù)目。

[nGenes, nSamples] = size(expr_cns_gcrma_eb)

nGenes =

7129

nSamples =

42

可以用基因符號來代替探針集的ID用于標記基因表達值，HuGeneFl芯片的基因符號在一個包含Java哈希表的MAT文件中。

load HuGeneFL_genesymbol_hashtable;

為hu6800genesymbol_hashtable變量創(chuàng)建一個基因表達值的基因符號的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);

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

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

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

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

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

[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

識別差異基因表達

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

從42個樣本中提取10個MD和10個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檢驗是檢測兩組變量之間顯著性差異的標準統(tǒng)計檢驗，對每個基因執(zhí)行t檢驗以識別MD樣本和Mglio樣本基因表達值的顯著性差異，可以通過t得分的正態(tài)分位圖和t得分及p值的直方圖來研究檢驗結果。

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

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

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

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

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

cutoff = 0.05;

sum(pvaluesCorr < cutoff)

ans =

2007

用mafdr對每個檢驗估計FDR和q值，pi0 是研究中真的空假設的總的分數(shù)，可以通過bootstrap或立方多項式擬合來估計模擬的空分布，也可以手動設置λ值來估計pi0 。

figure;

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

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

確定差異表達基因數(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)來注釋差異表達基因，從上面的分析結果中查找上調(diào)基因，從Gene Ontology Current Annotations下載人類注釋(gene_association.goa_human.gz)，解壓縮并存入當前目錄。

找出上調(dià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ù)庫為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'});

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

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

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

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

numel(HGgenes)

ans =

74298

并非所有的HuGeneFL芯片上的5758個基因都有注釋，通過比較基因符號列表和GO中的基因符號列表可以看出其是否已被注釋，可以跟蹤注釋基因數(shù)和每個GO術語關聯(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)計顯著的GO術語，對于每一個GO術語，p-值表示關聯(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術語并選擇與具體的分子功能相關的術語來構建一顆包括其祖先的子本體論樹，用biograph函數(shù)來可視化顯示它，也可為節(jié)點著色，本文中紅色的節(jié)點是zui顯著的，蘭色的節(jié)點zui不顯著，由于人類基因注釋文件的頻繁更新可能返回不同的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.

<img alt="在Matlab中探索基因表達數(shù)據(jù)分析" 在matlab中探索基因表達數(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)差異表達基因

通過KEGG's SOAP Web Service或將基因符號列表發(fā)送到KEGG's Web query tool可以從KEGG代謝途徑數(shù)據(jù)庫查詢差異表達基因的代謝途徑信息。

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.