library(survival)
library(ggplot2)
library(survminer)
## Loading required package: ggpubr
## Loading required package: magrittr
library(grid)
library(mvtnorm)
library(stats4)
library(modeltools)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(strucchange)
## Loading required package: sandwich
library(sandwich)
library(party)
library(maxstat)
methy <- read.csv("TCGA_LAML_hMethyl450.csv", header=TRUE, sep=",")
head(methy)
## sampleID OS OS_IND cg10262052 rank_median cut_value
## 1 TCGA-AB-2802-03 365 1 -0.4375 Low Low
## 2 TCGA-AB-2803-03 792 1 -0.3023 High High
## 3 TCGA-AB-2804-03 2557 0 -0.3989 High High
## 4 TCGA-AB-2805-03 577 1 -0.3881 High High
## 5 TCGA-AB-2806-03 945 1 -0.4033 High High
## 6 TCGA-AB-2807-03 181 1 -0.3402 High High
stat <- maxstat.test(Surv(OS, OS_IND)~cg10262052, data=methy, smethod="LogRank", pmethod="exactGauss", abseps=0.01)
stat
##
## Maximally selected LogRank statistics using exactGauss
##
## data: Surv(OS, OS_IND) by cg10262052
## M = 3.3604, p-value = 0.01463
## sample estimates:
## estimated cutpoint
## -0.4242
plot(stat)
methygroup = ifelse(methy$cg10262052<=-0.4242, 0, 1)
methygroup = factor(methygroup)
levels(methygroup) = c("below -0.4242", "above -0.4242")
fit <- survfit(Surv(OS, OS_IND) ~ methygroup, data=methy)
fit
## Call: survfit(formula = Surv(OS, OS_IND) ~ methygroup, data = methy)
##
## n events median 0.95LCL 0.95UCL
## methygroup=below -0.4242 62 48 243 184 365
## methygroup=above -0.4242 117 68 608 486 854
plot(fit)
survdiff(Surv(OS, OS_IND)~methygroup, data=methy)
## Call:
## survdiff(formula = Surv(OS, OS_IND) ~ methygroup, data = methy)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## methygroup=below -0.4242 62 48 31 9.27 13.2
## methygroup=above -0.4242 117 68 85 3.39 13.2
##
## Chisq= 13.2 on 1 degrees of freedom, p= 0.000282
out=coxph(Surv(OS, OS_IND)~methygroup, data=methy)
out
## Call:
## coxph(formula = Surv(OS, OS_IND) ~ methygroup, data = methy)
##
## coef exp(coef) se(coef) z p
## methygroupabove -0.4242 -0.686 0.504 0.191 -3.59 0.00033
##
## Likelihood ratio test=12.2 on 1 df, p=0.000471
## n= 179, number of events= 116
ggsurvplot(fit, pval=TRUE, legent.title="cg10262052",
legend.labs=c("below cut-point", "above cut-point"),
xlab="Time(days)", break.time.by=500, main="TCGA LAML hMethyl 450 - cg10262052",
legend="right", palette=c("#00BFC4", "#F8766D"))
