Praktický příklad analýzy

Ukažme si na příkladu jak může probíhat analýza přežití v analýze genomických a proteomických dat.

Budeme pokračovat v práci se stejným příkladem experimentu porovnávajícího ER (estrogen receptor) pozitivní a ER negativní karcinomy prsu, který jsme používali v předcházejících kapitolách. Budeme používat výsledky z praktického příkladu analýzy z kapitoly "Porovnávání skupin" a to zejména:

- objekt top100 obsahující top 100 odlišně exprimovaných genů mezi skupinami ER+ a ER- (jak získat top100 najdete ZDE)

Načteme si nejdříve genové expresní data karcinomu prsu TransBig spolu s tabulkou klinických dat. Data lze stáhnout ZDE (nutno rozbalit archiv). Načteme také top 100 genů odlišně

> load(file="TransBIG/TransBIG-expression.rdata")
> clinical2 <- read.csv("TransBIG/TransBIG-Sample_info.csv")
> load("top100.rdata")

1) Chceme zjistit rozdíl v přežití mezi skupinami ER+ a ER-

Načteme funkci pro vykreslení Kaplan-Meierových křivek a balík pro analýzu přežití:

> source("surv.plot.R")
> library(survival)
Loading required package: splines
> rfs <- clinical2$Time.RFS..days
> rfs.event <- clinical2$Event.RFS
> surv.plot(survfit(Surv(rfs, rfs.event)~clinical2$ER.pos.1.yes.))

Odhadneme medián přežití:
> survfit(Surv(rfs, rfs.event)~clinical2$ER.pos.1.yes.)
Call: survfit(formula = Surv(rfs, rfs.event) ~ clinical2$ER.pos.1.yes.)

   4 observations deleted due to missingness
                          records n.max n.start events median 0.95LCL 0.95UCL
clinical2$ER.pos.1.yes.=0      75    75      75     38   4493    2192      NA
clinical2$ER.pos.1.yes.=1     174   174     174     79   5588    4273      NA

Je rozdíl v přežití statisticky významný?
> survdiff(Surv(rfs, rfs.event)~clinical2$ER.pos.1.yes.)
Call:
survdiff(formula = Surv(rfs, rfs.event) ~ clinical2$ER.pos.1.yes.)

n=249, 4 observations deleted due to missingness.

                            N Observed Expected (O-E)^2/E (O-E)^2/V
clinical2$ER.pos.1.yes.=0  75       38     31.8     1.208      1.67
clinical2$ER.pos.1.yes.=1 174       79     85.2     0.451      1.67

 Chisq= 1.7  on 1 degrees of freedom, p= 0.196

Není, i když je tam trend. Nyní to zkusíme se skupinami definovanými shlukováním, kterým jsme v kapitole Objevování skupin tvořili na datech heatmapu. Heatmapa sice vnitřně aplikovala shlukovací algoritmus, ale my to musíme udělat znovu abychom získali dendrogram, ze kterého pak vyextrahujeme příslušnost ke skupinám.

> eukl2 <- dist(X[,intersect(colnames(X),top100)])
> hclust2 <- hclust(eukl2)
> plot(hclust2)

> groups2 <- cutree(hclust2, k=2)# delime na dve skupiny
> head(groups2)
4000471 4000473 4000475 4000477 4000480 4000482
      1       2       2       2       2       1

Nyní analýza přežití s novými skupinami:
> surv.plot(survfit(Surv(rfs, rfs.event)~groups2))

Další úlohy k procvičení:
1) Udělejte podobnou analýzu s top200 a top500 geny
2) Poté udělejte podobnou analýzu s HER2

Nyní vše zkusíme spojit. Modelování pomocí Coxova modelu. Aplikujeme normalizaci, median clustering a IQR dělení.
Vypočítáme madián a IQR exprese každého genu.
> med<-median(X)
> iqr<-IQR(X)
> X2<-(X-med)/iqr

Má gen ESR1 nějaký vliv na dobu přežití do relapsu (RFS)?
> coxph(Surv(rfs, rfs.event)~X2[,"ESR1"])
Call:
coxph(formula = Surv(rfs, rfs.event) ~ X2[, "ESR1"])

                 coef exp(coef) se(coef)      z    p
X2[, "ESR1"] -0.00414     0.996   0.0231 -0.179 0.86

Likelihood ratio test=0.03  on 1 df, p=0.858  n= 253, number of events= 118

ESR1 nemá vliv na RFS (p = 0.86)

Nyní hledáme geny, které mají významný vliv na relepse-free-survival:
> surv.res<-c()
> for (i in c(1:ncol(X2)))
+ {
+ temp<-coxph(Surv(rfs, rfs.event)~X2[,i])
+ a<-summary(temp)
+ surv.res<-rbind(surv.res, c(colnames(X2)[i],a$coefficients))
+ }
> colnames(surv.res)<-c("gene",colnames(a$coefficients))