-
Notifications
You must be signed in to change notification settings - Fork 0
/
xtile.R
89 lines (67 loc) · 2.84 KB
/
xtile.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
library(survminer)
library(survival)
library(ggplot2)
myfile = list.files("Xtile")
dir = paste("./Xtile/", myfile, sep="")
n = length(dir)
# Take cut-off one by one and draw pictures one by one ---------------------------------------------------------
i <- 1
data = read.table(file = dir[i],header=T, check.names = FALSE)
dir[i]
gene_name <- colnames(data)[3]
exprSet <- data.frame(t(data))
## Set cutoff value
alpha <- 0.500 # Obtained by Xtile
risk_score <- t(as.matrix(exprSet[gene_name,]))
cut_off <- rep(as.numeric(quantile(exprSet[gene_name,],alpha)), dim(exprSet)[2])
# data$time <- data$time/12
data$riskscore <- ifelse(risk_score > cut_off, 'high','low')
table(data$risk_score)
fit <- survfit(Surv(time, status)~riskscore, data = data)
p <- ggsurvplot(fit, data = data,
conf.int = F,
# surv.median.line = "hv",
risk.table = TRUE,
tables.height= 0.25,
cumcensor = T,
legend = c(0.83,0.95),
# P value
pval = TRUE,
pval.size=6,
font.pval= c(14, "bold", "black"),
pval.coord = c(0.00, 0.05),
# legend
legend.title = '', # gene_name
legend.labs=c("High risk", "Low risk"),
font.legend= c(14, "plain", "black"),
# font.main = c(100, "bold", "black"),
# xlim = c(0,72), # present narrower X axis, but not affect
# survival estimates.
palette=c("red", "blue"),
font.x = c(14, "plain", "black"),
font.y = c(14, "plain", "black"),
font.tickslab = c(14, "plain", "black"),
xlab = "Time in years", # customize X axis label. year
break.time.by = 2
) # break X axis in time intervals by 500.
# Ìí¼ÓHRºÍCI
res_cox <- coxph(Surv(time, status) ~riskscore, data=data)
HR <- round(summary(res_cox)$conf.int[1],2)
ci_l <- round(summary(res_cox)$conf.int[3],2)
ci_r <- round(summary(res_cox)$conf.int[4],2)
p1 <- p
p1$plot = p1$plot +
ggplot2::annotate("text",x = 2.13, y = 0.3, label = paste("HR : ", HR), size = 5) +
ggplot2::annotate("text",x = 3.88, y = 0.2,
label = paste("(","95%CI : ", ci_l,"-",ci_r,")", sep = ""), size = 5)
library(tidyverse)
name <- dir[i]
name1 <- str_split_fixed(name, "./", 2);
name2 <- str_split_fixed(name1[2], "/", 2)
name3 <- str_split_fixed(name2[2], "[.]", 2);
name4 <- name3[1]
name5 <- str_c(name4,"_ex")
path <- paste("./Xtilefigure/", paste(name5,".pdf", sep=""), sep="")
pdf(path, width = 4.8, height = 6, onefile = FALSE)
p1
dev.off()