E-learning 2. Analýza dat 2.2. Analýza vysokopokryvných genomických dat 2.2.1. DNA mikročipy 2.2.1.9. Příklad
Jeden z cílů mikročipových experimentů je nalézt nové biologické markery, které by mohly pomoct lékařům určit co nejdříve správnou diagnózu. Zvláště v onkologii je důležité zachytit nemoc v začátcích, protože čím nižší je stupeň při diagnostice, tím je vyšší pravděpodobnost remise nebo alespoň částečné remise. Stádium rakoviny je nejsilnější prediktor přežití. Ve vyšších stupních se vyskytuje vysoká mortalita. Nalezené geny, které jsou více nebo méně exprimované v rakovinových tkáních pomohou ve vytváření rychlejší a přesnější diagnózy, která je s ohledem na léčbu velmi důležitá. Heterogenita některých typů nádorů vyžaduje rozdílnou léčbu, ale je těžké je rozlišit.
Ukázkový datový soubor
V tomto příkladě budeme srovnáva dvě skupiny kolorektálních karcinomů. Jedna skupina obsahuje vzorek mikrosatelitního stabilního kolorektálního nádoru (MSS) a druhá obsahuje mikrosatelitní nestabilní kolorektální nádor (MSI). Datový soubor pochází z databáze Gene Expression Omnibus. Celkem obsahuje 38 mikrosatelitně stabilních a 13 mikrosatelitně nestabilních kolorektálních nádorů, každý od jiného pacienta. Všechny mRNA vzorky byly hybridizovány na Affymetrix HuGeneFL čipu se 7129 sondových souborů.
Normalizovaný datový soubor je ke stažení zde.
Ke sledování analýzy je potřebné mít spuštěný R program a kopírovat příkazy do R konzole. Pokud nejste seznámeni s R a Bioconductorem, přečtěte si kapitolu o Softwaru.
Odlišná exprese
V první řadě musíme načíst data.
> load(“GSE11543.rdata”)
Ke zkontrolování, zda jsou data načtena správně použijeme následující příkaz, který zobrazuje seznam všech proměnných.
> ls()
[1] "Xnrm"
Samozřejmě máme pouze jednu proměnnou, nazvanou Xnrm. Tato proměnná je objektem třídy ExpressionSet, která uchovává informace o sondových sobuborech na mikročipovém sklíčku, vzorcích a intenzitách genové exprese. Zavolání proměnné vrátí následující sumarizaci:
> Xnrm
Loading required package: Biobase
Loading required package: tools
Welcome to Bioconductor
Vignettes contain introductory material. To view, type
'openVignette()'. To cite Bioconductor, see
'citation("Biobase")' and for packages 'citation(pkgname)'.
ExpressionSet (storageMode: lockedEnvironment)
assayData: 7129 features, 51 samples
element names: exprs
phenoData
sampleNames: GSM290693, GSM290694, ..., GSM290743 (51 total)
varLabels and varMetadata description:
sample: arbitrary numbering
featureData
featureNames: A28102_at, AB000114_at, ..., Z97074_at (7129 total)
fvarLabels and fvarMetadata description: none
experimentData: use 'experimentData(object)'
Annotation: hu6800
AssayData ukládají matici genové exprese, phenoData obsahují informace o vzorcích a všechny meta-popisy datového souboru, featureData ukládají informace o názvech sondových souborů a experimentData nám poskytují informace o typu mikročipového sklíčka vzorků, které se hybridizují. Pro více informací o ExpressionSet class, zadejte:
>?ExpressionSet
Matice intenzit genové exprese může být získána z tohoto objektu použitím exprs() funkce:
> data.exprs = exprs(Xnrm)
Nyní potřebujeme definovat skupinovou proměnnou, která poskytne informace jaký vzorek patří do jaké skupiny, které právě porovnáváme. Skupinová proměnná může být extrahována ze souboru ExpressionSet.
Pokud nejsou z nějakého důvodu poskytnuty skupinové informace, můžeme je nahrát z externího souboru.
> groups = read.csv2("./Pheno.11543.csv", sep=";", header=FALSE)
Soubor s anotacemi o sondách by měl být nahraný ze samostatného souboru:
> annotation = read.csv2("./Annotation.csv", sep=";", header=FALSE)
Nyní máme všechny údaje potřebné ke spuštění analýzy. V tomto příkladě použijeme pro analýzu odlišné exprese 4 metody.
1) Dvouvýběrový t-test
R funkce pro T-test je implementována v základním statistickém balíku.
> ?t.test
Skupiny k testování mohou být definovány buď každá jednou proměnnou:
> a = data.exprs[1,] [groups[,2]==”S”] #selecting expression values of the first probeset of the group MSS
> b = data.exprs[1,] [groups[,2]==”I”] #selecting expression values of the first probeset of the group MSI
> t.test(a,b)
nebo můžeme použít modelový postup
> t.test(data.exprs[1,]~groups[,2])
Protože budeme testovat každou z tísíců sond v našem datovém souboru, není rozumné stejný příkaz nahrávat do R konzole tisícekrát. Můžeme to zautomatizovat vytvořením cyklu. Řekněme, že chceme uložit všechny p-hodnoty do jednoho vektoru, a tak definujeme proměnnou
> result.p.val = c()
Dále vytvoříme “for cycle” a v každém cyklu přidáme do našeho vektoru p-hodnotu z každho sondového souboru
> for (i in c(1:nrow(data.exprs)))
{
result.p.val = c(result.p.val, t.test(data.exprs[i,]~groups[,2])$p.value)
}
Tímto způsobem můžeme vytvořit více vektorů ukládajících např. T-statistiky atd.
Nicméně existuje rychlejší způsob, jak získat výsledky ze všech srovnání, než je cyklus. Funkce apply nebo tapply jsou rychlejší a poskytují nám výsledky ze všech sondových souborů jen jedním příkazem:
> ?apply
> result.p.val = apply(data.exprs, 1, FUN=function(x) t.test(x~groups[,2])$p.value)
Nyní to musíme upravit na vícenásobný testovací problém blíže popsaný v kapitole….
K tomuto účelu existují v R a Bioconductoru různé funkce. Jedna z nich je funkce p.adjust nabízející 7 metod pro korekci p-hodnoty. K upravení našich p-hodnot použijeme Benjamini&Hochberg [1995] metodu pro kontrolu FDR (false discovery rate) a také Boneferroniho korekci, která kontroluje FWER (familywise error rate).
> result.p.val adjust.bonf = p.adjust(result.p.val, method="bonferroni")
> result.p.val.adjust.BH = p.adjust(result.p.val, method="BH")
K porovnání výsledků těchto dvou korekcí zobrazíme histogramy
> png("./Comparison.of.t.test.p.val.adjustment.png", width=1100, height=400)
> par(mfrow=c(1,3))
> hist(result.p.val, ylim=c(0,7000), xlab = "P.values", main ="P.values without correction", cex.main=2, cex.lab=1.5, cex.axis=1.5)
> hist(result.p.val.adjust.bonf, ylim=c(0,7000), xlab = "P.values", main = "P.values-Bonferroni correction", cex.main=2, cex.lab=1.5, cex.axis=1.5)
> hist(result.p.val.adjust.BH, ylim=c(0,7000), xlab = "P.values", main = "P.values-Benjamini & Hochberg correction", cex.main=2, cex.lab=1.5, cex.axis=1.5)
> dev.off()
null device
1
Na obrázku jsou tři histogramy, vlevo je histogram původních p-hodnot, uprostřed p-hodnoty po Bonferroniho korekci a napravo histogram p-hodnot po B&H korekci. Je viditělné, že mezi těmito dvěma úpravami jsou velké rozdíly. Bonferroni korekce je pro tak vysoký počet pozorování příliš konzervativní, k dalšímu výpočtu jsme tedy vybrali B&H korekci.
> res.t.test.final = cbind(rownames(data.exprs),result.p.val.adjust.BH)
T-test je založen na několika předpokladech. Pokud nemohou být potvrzeny, musíme použít parametrické testy. Pro naše data zkontrolujme předpoklad normality:
> shap.wilk.res.MSS = apply(data.exprs[,groups[,2]==”MSS”], 1, FUN=function(x) shapiro.test(x)$p.value)
> shap.wilk.res.MSI = apply(data.exprs[,groups[,2]==”MSI”], 1, FUN=function(x) shapiro.test(x)$p.value)
> SW.MSS.p.adjust = p.adjust(shap.wilk.res.MSS, method="BH")
> SW.MSI.p.adjust = p.adjust(shap.wilk.res.MSI, method="BH")
> length(which(SW.MSS.p.adjust<0.05))
[1] 1694
> length(which(SW.MSI.p.adjust<0.05))
[1] 162
Shapiro Wilkův test normality ukazuje, že na hladině významnosti α = 0,05 nemůžeme potvrdit normalitu dat pro 162 sondových souborů v MSI skupině a pro 1694 sondových souborů v MSS skupině.
V takovém případě musíme zvážit neparametrický přístup. Druhá metoda, kterou použijeme je tedy Wilcoxonův pořadový test. Zpracování je velmi podobné jako v případě t-testu, postačí nám odpovídající R funkce - wilcox.test.