1 / 25

Analýza dat v R

Analýza dat v R. Analýza reálneho príkladu - microarrays Eva Budinská 28-30. apríl 2011 Brno. Príklad. Dátový súbor: Z databáze GEO (GSE11543) 51 pacientov s nádorom kolorekta Namerané klinické parametre: Vek, pohlavie

tanek-romero
Download Presentation

Analýza dat v R

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Analýza dat v R Analýza reálneho príkladu - microarrays Eva Budinská 28-30. apríl 2011 Brno

  2. Príklad • Dátový súbor: • Z databáze GEO (GSE11543) • 51 pacientov s nádorom kolorekta • Namerané klinické parametre: • Vek, pohlavie • Mikrosatelitová nestabilita ( MSI – nestabilní (lepšie prežitie), MSS – stabilní (horšie prežitie)) • K dispozícii hodnoty aktivity génov reprezentovaných 7129 sondami, merané u každého pacienta • Cieľ: Zistiť, ktoré gény sú zasiahnuté mikrosatelitovou nestabilitou (zistiť teda, aké sú rozdiely medzi skupinami MSI a MSS v aktivite génov) • Dátové súbory: • clinical = read.table("clinical.csv", header=TRUE, sep=",") # klinické informácie o pacientoch • annotation = read.table("Annotation.csv", header=TRUE, sep=",", skip=12) #informácia o génoch (názov, chromozóm, funkcia...) • load("colon.Rdata") # hodnoty aktivity génov u pacientov

  3. Dátové súbory • clinical - klinické informácie o pacientoch ID MSS Sex DukesStage Age 1 William Carbines MSS Male B 78 2 Leontine Aubert MSS Female A 73 • Ivan Mineff MSS Female B 65 • colon– informácie o aktivite génov u pacientov (čím vyššie číslo, tým aktívnejší gén) William Carbines Leontine Aubert Ivan Mineff GABRA3 7.641879 7.616971 7.278516 OMD 6.818050 6.675985 6.266338 IFI44L 7.578773 8.198093 7.873286 • annotation – informácie o génoch a ich pozícii na chromozóme, stĺpce c(1,15,16) Chromosomal.Location Gene.Title GABRA3 chrXq28 gamma-aminobutyric acid (GABA) A receptor, alpha 3 OMD chr9q22.31 osteomodulin IFI44L chr1p31.1 interferon-induced protein 44-like

  4. Schéma analýzy • Cieľ: Zistiť ktoré gény sú rozdielne aktívne medzi skupinami MSI a MSS • Analýza: • Porovnanie aktivity je ekvivalentné porovnaniu priemerov číselných hodnôt aktivity jednotlivých génov medzi skupinami MSI a MSS – testovanie hypotéz Aplikujeme teda • Parametrický T-test • Neparametrický Mann-Whitney test • Hypotézy testujeme pre všetky gény! (teda 5311 krát!!!) – vzniká problém mnohonásobného porovnávania (u hladiny významnosti 5% máme 5311*0.05 = 265 falošne pozitívnych génov). Musíme teda vykonať korekciu p-hodnôt z testovania • Nonferonniho korekcia • FDR • Ďalej urobíme zhlukovanie génov, ktoré sme našli ako rozdielne aktívne medzi MSI a MSS, aby sme našli funkčné skupiny • Vykreslíme heatmapu!

  5. Načítanie dát • load("colon.Rdata") • annotation=read.table(„annotation.csv", header=TRUE, sep=",", rownames=1) • clinical = read.table("clinical.csv", header=TRUE, sep=",",rownames=1) Skontrolujeme načítanie dát: • dim(colon) • dim(annotation) • dim(clinical) • colon[1:5,1:5] • annotation[1:5,] • clinical[1:5,1:3]

  6. Schéma analýzy • Cieľ: Zistiť ktoré gény sú rozdielne aktívne medzi skupinami MSI a MSS • Analýza: • Porovnanie aktivity je ekvivalentné porovnaniu priemerov číselných hodnôt aktivity jednotlivých génov medzi skupinami MSI a MSS – testovanie hypotéz Aplikujeme teda • Parametrický T-test • Neparametrický Mann-Whitney test • Hypotézy testujeme pre všetky gény! (teda 5311 krát!!!) – vzniká problém mnohonásobného porovnávania (u hladiny významnosti 5% máme 5311*0.05 = 265 falošne pozitívnych génov). Musíme teda vykonať korekciu p-hodnôt z testovania • Nonferonniho korekcia • FDR • Ďalej urobíme zhlukovanie génov, ktoré sme našli ako rozdielne aktívne medzi MSI a MSS, aby sme našli funkčné skupiny • Vykreslíme heatmapu!

  7. Testovanie hypotéz • informácia o aktivite génov je uložená v súbore colon • gény sú uložené ako riadky, preto sa na ne budeme odvolávať pomocou colon[i, ] • informácia o príslušnosti pacientov k skupine MSI a MSS je v súbore clinical • uložme si informáciu o MSI a MSS do premennejgroup group = clinical$MSS • Funkcia preT-test: t.test() • Možný spôsob vstupných hodnôt • t.test(hodnoty v skupine 1, hodnoty v skupine 2) • t.test(hodnoty~skupina) Konkrétne pre otestovanie i-tého génu (riadku) • t.test(colon[i,group=="MSS"], colon[i,group=="MSI"]) • t.test(colon[i,]~group)

  8. Testovanie hypotéz i=1 vysledok = t.test(colon[i,]~group) vysledok #zobrazime výsledok testu Ako vyextrahujeme p-hodnotu??? summary(vysledok) #ukáže nám všetky elementy výsledného objektu Je medzi nimi i hodnota p.value. Získame ju teda odvolaním sa na ňu cez $p.value: vysledok$p.value [1] 0.5377606 Nakreslíme si obrázok, aby sme videli rozdiely medzi skupinami boxplot(colon[i,]~group) Testová úloha: Pridajte do obrázku p-hodnotu zaokrúhlenú na 3 desatinné miesta a názov génu do titulku grafu. Výsledný graf spolu s príkazom ktorým ste ho vytvorili uložte do Wordu. [3 body] Nápoveda: použite paste()

  9. Výsledný boxplot

  10. Testovanie hypotéz II. • Chceme spočítať T-test pre každý gén. Aplikujeme teda funkciu apply na súbor colon po riadkoch, a z t.testu vyberieme len p-hodnoty: vsetky.geny.ttest = apply(colon, 1, FUN=function(x) t.test(x~ group)$p.value) • Chceme zistiť, ktoré gény sú štatisticky významné, a teda ktoré p-hodnoty sú menšie ako 0.05 vyznamne.geny.t.test = vsetky.geny.ttest<0.05 #podmienkou vytvoríme vektor false/true Úloha: Ako zistíme, koľko ich je? Pozrieme si ktoré gény sú to konkrétne: rownames(annotation)[vyznamne.geny.t.test]

  11. Schéma analýzy • Cieľ: Zistiť ktoré gény sú rozdielne aktívne medzi skupinami MSI a MSS • Analýza: • Porovnanie aktivity je ekvivalentné porovnaniu priemerov číselných hodnôt aktivity jednotlivých génov medzi skupinami MSI a MSS – testovanie hypotéz Aplikujeme teda • Parametrický T-test • Neparametrický Mann-Whitney test • Hypotézy testujeme pre všetky gény! (teda 5311 krát!!!) – vzniká problém mnohonásobného porovnávania (u hladiny významnosti 5% máme 5311*0.05 = 265 falošne pozitívnych génov). Musíme teda vykonať korekciu p-hodnôt z testovania • Nonferonniho korekcia • FDR • Ďalej urobíme zhlukovanie génov, ktoré sme našli ako rozdielne aktívne medzi MSI a MSS, aby sme našli funkčné skupiny • Vykreslíme heatmapu!

  12. Mnohonásobné porovnanie • Keď testujeme 5311 krát a používame hladinu významnosti 5%, tak môžeme čakať 0.05*5311 = 265 falošne pozitívnych génov. • Korekcia p-hodnôt na tento problém, funkcia p.adjust() • Bonferroni: bonferroni.p.hodnoty = p.adjust(vsetky.geny.ttest, method="bonferroni") • FDR: fdr.p.hodnoty = p.adjust(vsetky.geny.ttest, method="fdr") Koľko máme významných génov teraz?? sum(fdr.p.hodnoty<0.05) sum(bonferroni.p.hodnoty<0.05) Testová úloha: Nájdite funkciu pre neparametrický Wilcoxon test[1 bod] a zanalyzujte dáta touto funkciou spolu s korekciou p-hodnôt na mnohonásobné porovnanie metódou FDR. Zapíšte koľko génov je významných na 0.05 hladine významnosti, a to po korekcii na mnohonásobné porovnávanie. [2 body] Ktorý gén má najmenšiu p-hodnotu (upravenú FDR)? – napíšte meno a zistite na ktorom chromozóme sa nachádza. [2 body]

  13. Schéma analýzy • Cieľ: Zistiť ktoré gény sú rozdielne aktívne medzi skupinami MSI a MSS • Analýza: • Porovnanie aktivity je ekvivalentné porovnaniu priemerov číselných hodnôt aktivity jednotlivých génov medzi skupinami MSI a MSS – testovanie hypotéz Aplikujeme teda • Parametrický T-test • Neparametrický Mann-Whitney test • Hypotézy testujeme pre všetky gény! (teda 7129 krát!!!) – vzniká problém mnohonásobného porovnávania (u hladiny významnosti 5% máme 7129*0.05 = 356 falošne pozitívnych génov). Musíme teda vykonať korekciu p-hodnôt z testovania • Nonferonniho korekcia • FDR • Ďalej urobíme zhlukovanie génov, ktoré sme našli ako rozdielne aktívne medzi MSI a MSS, aby sme našli funkčné skupiny • Vykreslíme heatmapu!

  14. Zhlukovanie • T-testom sme vybrali množinu 320 génov, ktoré sú štatisticky významné medzi skupinami MSI a MSS na hladine významnosti 0.05 po korekcii na FDR. • Chceme teraz vedieť, ktoré z týchto génov tvoria spoločné zhluky • Takisto chceme vedieť, či vybrané gény naozaj dobre rozdeľujú pacientov na skupiny MSI a MSS • Urobíme preto hierarchické zhlukovanie jednak génov a jednak pacientov. • hclust – funkcia pre hierarchické zhlukovanie z balíka cluster • parametre: d – matica vzdialeností medzi objektami method – algoritmus zhlukovania, na výber: "ward","single","complete", "average", "mcquitty", "median", "centroid" unit – ako sa majú zobraziť objekty v dendrograme labels – názvy listov stromu (defaultne názvy alebo čísla riadkov dát)

  15. Metrika vzdialeností I. • Ak chceme zhlukovať, najprv musíme definovať maticu vzdialeností d, ktorá je základným vstupom do funkcie hclust. Použijeme dve: • Euklideovskú vzdialenosť – funkcia dist() dist(x,method ) parametre: x – dátová tabuľka, objekty medzi ktorými počítame vzdialenosť musia byť v riadkoch method – metrika vzdialenosti euclid.pacient = dist(t(colon[fdr.p.hodnoty<0.05,]), method="euclidean")# vzdialenosť medzi vzorkami euclid.geny = dist(colon[fdr.p.hodnoty<0.05,], method="euclidean")# vzdialenosť medzi vzorkami

  16. Metrika vzdialeností II. • Metrika založená na Spearmanovej korelácii – cor() cor(x, method) – vypočíta koreláciu medzi stĺpcami matice matice x – matica method – "pearson" alebo "spearman" cor.pacient = cor(colon[fdr.p.hodnoty<0.05,], method="spearman")#počítame koreláciu medzi pacientami Potrebujeme ale vzdialenosť!!! cor.pacient.dist = as.dist(1-cor.pacient)# prekonvertujeme korelácie na maticu vzdialeností cor.geny = cor(t(colon[fdr.p.hodnoty<0.05,]), method="spearman") #počítame koreláciu medzi génmi cor.geny.dist = as.dist(1-cor.geny) # prekonvertujeme na vzdialenosť

  17. Hierarchické zhlukovanie • Zhlukujeme pacientov algoritmom average linking na euklideovskej vzdialenosti: hclust.pac.eucl.avg = hclust(euclid.pacient, method="average") a na vzdialenosti odvodenej od korelácie hclust.pac.spear.avg = hclust(cor.pacient.dist, method="average") • Podobne zhlukujeme gény: na euklideovskej vzdialenosti: hclust.geny.eucl.avg = hclust(euclid.geny, method="average") a na vzdialenosti odvodenej od korelácie hclust.geny.spear.avg = hclust(cor.geny.dist, method="average")

  18. Dendrogram • Máme teda 4 objekty hierarchického zhlukovania hclust.pac.eucl.avg #pacienti na euklideovskej vzd. hclust.pac.spear.avg #pacienti na koreláciách hclust.geny.eucl.avg #gény na euklideovskej vzd. hclust.geny.spear.avg #gény na koreláciách • Ak chceme výsledky zobraziť, musíme vytvoriť dendrogramy dendr.hclust.pac.eucl.avg = as.dendrogram(hclust.pac.eucl.avg, labels=group) dendr.hclust.pac.spear.avg = as.dendrogram(hclust.pac.spear.avg) dendr.hclust.geny.eucl.avg = as.dendrogram(hclust.geny.eucl.avg) dendr.hclust.geny.spear.avg = as.dendrogram(hclust.geny.spear.avg) #gény na koreláciách

  19. Dendrogram - zobrazenie • Vykreslíme dendrogramy do jedného grafu pre pacientov: par(mfrow=c(2,1),cex=0.5) plot(dendr.hclust.pac.eucl.avg) plot(dendr.hclust.pac.spear.avg) savePlot("dendrogram.pacienti.png", type="png") a pre gény windows(28,14) par(mfrow=c(2,1), cex=0.5) plot(dendr.hclust.geny.eucl.avg) plot(dendr.hclust.geny.spear.avg) savePlot("dendrogram.geny.png", type="png")

  20. Schéma analýzy • Cieľ: Zistiť ktoré gény sú rozdielne aktívne medzi skupinami MSI a MSS • Analýza: • Porovnanie aktivity je ekvivalentné porovnaniu priemerov číselných hodnôt aktivity jednotlivých génov medzi skupinami MSI a MSS – testovanie hypotéz Aplikujeme teda • Parametrický T-test • Neparametrický Mann-Whitney test • Hypotézy testujeme pre všetky gény! (teda 7129 krát!!!) – vzniká problém mnohonásobného porovnávania (u hladiny významnosti 5% máme 7129*0.05 = 356 falošne pozitívnych génov). Musíme teda vykonať korekciu p-hodnôt z testovania • Nonferonniho korekcia • FDR • Ďalej urobíme zhlukovanie génov, ktoré sme našli ako rozdielne aktívne medzi MSI a MSS, aby sme našli funkčné skupiny • Vykreslíme heatmapu!

  21. Heatmapa • Kreslíme heatmapu Nastavíme farebnú škálu zelenej cez čiernu po červenú, kde zelená budú málo aktívne gény a červená veľmi aktívne gény. Čierna farba znamená priemernú expresiu génu. farby = colorRampPalette(c("green", "black","red"))(15) heatmap(colon[fdr.p.hodnoty<0.05,]) heatmap(colon[fdr.p.hodnoty<0.05,], col=farby) pridáme farby skupín MSI/MSS (parameter ColSideColors): heatmap(colon[fdr.p.hodnoty<0.05,], col=farby, ColSideColors=as.character(as.numeric(group)))

  22. Heatmapa a vlastné zhlukovanie • Ak chceme heatmapu na základe nášho vlastného zhlukovania, musíme definovať parametre Rowv (pre dendrogram riadkov - génov) a Colv (pre dendrogram stĺpcov-pacientov) • Pre pacientov použijeme dendrogram zo zhlukovania na spearmanovej korelácii dendr.hclust.pac.spear.avg • Pre gény použijeme dendrogram zo zhlukovania na euklideovskej vzdialenosti dendr.hclust.geny.eucl.avg: heatmap(colon[fdr.p.hodnoty<0.05,], col=farby, ColSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.pac.spear.avg, Rowv= dendr.hclust.geny.eucl.avg) Zadáme dendrogram z korelácií i u génov: heatmap(colon[fdr.p.hodnoty<0.05,], col=farby, ColSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.pac.spear.avg, Rowv=dendr.hclust.geny.spear.avg)

  23. Upravujeme heatmapu • Chceme vidieť názvy génov (annotation$Gene.Symbol), zadefinujeme preto parameter labRow a zmenšíme písmo v riadkoch – parameter cexRow: heatmap(colon[fdr.p.hodnoty<0.05,], col=farby, ColSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.pac.spear.avg, Rowv= dendr.hclust.geny.spear.avg, labRow=annotation$Gene.Symbol[fdr.p.hodnoty<0.05], cexRow=0.5) Testová úloha: Zmeňte v matici označenie pacientov na označenie skupín MSI/MSS. Zapíšte príkaz a pripojte výsledný graf. [1 bod] • Ak nechceme zhlukovanie riadkov alebo stĺpcov, zadáme hodnotu NA do parametrov Colv alebo Rowv heatmap(colon[fdr.p.hodnoty<0.05,], col=farby, ColSideColors=as.character(as.numeric(group)), Colv= NA, Rowv= dendr.hclust.geny.spear.avg, labRow=annotation$Gene.Symbol[fdr.p.hodnoty<0.05], cexRow=0.5)

  24. Upravujeme heatmapu II. Heatmapa je zle čitateľná, uložíme ju preto príkazom png a zadefinujeme väčšiu veľkosť: png("heatmapa.png", height=1200, width=520) heatmap(colon[fdr.p.hodnoty<0.05,], col=farby, ColSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.pac.spear.avg, Rowv= dendr.hclust.geny.spear.avg, labRow=annotation$Gene.Symbol[fdr.p.hodnoty<0.05], cexRow=0.5) dev.off() Toto by bol štandardný postup, ale nefunguje to, pretože heatmapa tejto funkcie je zadefinovaná cez layout a nie je tu možné meniť ho a teda a ni zväčšovať ju. Ak sa Vám nechce prerábať zdrojový kód, skúste funkciu heatmap.2 z gplots, tá už má viac možností!

  25. Ďalšie testové úlohy • Okrem troch predchádzajúcich testových úloh vyriešte nasledovné: • Pridajte do súboru clinical premennú ktorá bude kategorizovať expresnú hodnotu génu “MSH6“ na dve skupiny podľa priemernej hodnoty a nazvite ju MSH6kat. [4 body]. • V súbore clinical prekódujte premennú Age na skupiny <50, 50-70, >70. [2 body] • Poznámka: Ku každej úlohe pripojte kód, ak je výsledkom graf, tak i graf. • Celkový počet bodov za príklady 15 • Známkovanie: • viac ako 13 A • 12-10 B • 9-8 C • 7-6 D • 5-4 E • 3 a menej F Bonus: Účasť na celej výuke (prípadne ospravedlnenie neprítomnosti) + 2 body

More Related