-
Notifications
You must be signed in to change notification settings - Fork 16
/
survival.R
162 lines (147 loc) · 5.65 KB
/
survival.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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
# MSomicsQC::surv_designation is a fucntion that creates a "surv_DF" attribute with elements:
# 1) t_death - defines column name in f_data
# 2) t_progress
# 3) ind_death
# 4) ind_progress
# 5) covariates
#
# Must appear together:
#
# t_death & ind_death (overall survival - OS option)
# OR
# t_death, t_progress, ind_death AND ind_progress (progression free survival - PFS option)
#
# covariates is always optional - two max for now
#
# ---------
# Once Lisa/Kelly makes the above, write function to create Kaplan-Meier Curves.
# Only account for categorical variables when making these curves (separate curves for each level of each
# covariate)
#
# Functions skeleton:
# plot_km <- function(omicsData,...Other figure options from pmartRqc plot functions...) - produce KM plot
# summary_km <- function(omicsData, percent = ?) - return something like summary.survfit where "percent" gives a
# survival prob and they want "time" variable returned
#' Basic survival analysis function
#'
#' Implements overall survival analysis or progression-free survival analysis, depending upon the datatypes supplied to
#' surv_designation, and return the "survfit" object
#'
#' @param omicsData A pmartR data object of any class, which has a `group_df` attribute that is usually created by the `group_designation()` function
#' @return if fitted survival analysis object is returned
#'
#' @examplesIf requireNamespace("pmartRdata", quietly = TRUE)
#' \dontrun{
#' library(MSomicsSTAT)
#' library(OvarianPepdataBP)
#'
#' # Basic analysis without covariates
#' attr(tcga_ovarian_pepdata_bp, "survDF") <- list(t_death = "survival_time",
#' ind_death = "vital_status")
#' sfit <- fit_surv(tcga_ovarian_pepdata_bp)
#' plot(sfit)
#'
#' # Add some covariate information
#' attr(tcga_ovarian_pepdata_bp, "survDF") <- list(
#' t_death = "survival_time",
#' ind_death = "vital_status",
#' covariates = "g__initial_pathologic_diagnosis_method_g1"
#' )
#' sfit <- fit_surv(tcga_ovarian_pepdata_bp)
#' plot(sfit, col = c(1, 2))
#' }
#'
#' @export
fit_surv <- function(omicsData) {
if (!requireNamespace("survival", quietly = TRUE)) {
stop("Please install the 'survival' package.")
}
t_death <- attr(omicsData, "survDF")$t_death
t_progress <- attr(omicsData, "survDF")$t_progress
ind_death <- attr(omicsData, "survDF")$ind_death
ind_progress <- attr(omicsData, "survDF")$ind_progress
covariates <- attr(omicsData, "survDF")$covariates
surv_df <- data.frame(time = omicsData$f_data[, t_death], status = 0)
surv_df$status[grep("dead", omicsData$f_data[, ind_death], ignore.case = TRUE)] <- 1
if (is.null(covariates)) {
sfit <- survfit(Surv(time, status) ~ 1, data = surv_df)
} else {
form <- "Surv(time,status)~"
ncovars <- length(covariates)
for (i in 1:ncovars) {
surv_df <- cbind(surv_df, omicsData$f_data[, covariates[i]])
form <- paste0(form, covariates[i])
if (i < ncovars) {
form <- paste0(form, "+")
}
}
colnames(surv_df)[ncol(surv_df) - c((ncovars - 1):0)] <- covariates
sfit <- survfit(as.formula(form), data = surv_df)
}
return(sfit)
}
#' Basic survival analysis plot
#'
#' Implements overall survival analysis or progression-free survival analysis, depending upon the datatypes supplied to
#' surv_designation, and plot the resulting Kaplan-Meier curve.
#'
#' @param omicsData A pmartR data object of any class, which has a `group_df` attribute that is usually created by the `group_designation()` function
#'
#' @return a Kaplan-Meier curve
#'
#' @examplesIf requireNamespace("pmartRdata", quietly = TRUE)
#' \dontrun{
#' library(MSomicsSTAT)
#' library(OvarianPepdataBP)
#' attr(tcga_ovarian_pepdata_bp, "survDF") <- list(t_death = "survival_time",
#' ind_death = "vital_status")
#' plot_km(omicsData = tcga_ovarian_pepdata_bp)
#'
#' # Add covariates to "survDF" attribute
#' attr(tcga_ovarian_pepdata_bp, "survDF") <- list(
#' t_death = "survival_time",
#' ind_death = "vital_status",
#' covariates = "age_at_initial_pathologic_diagnosis"
#' )
#' plot_km(omicsData = tcga_ovarian_pepdata_bp)
#' }
#'
#' @export
#'
plot_km <- function(omicsData) {
sfit <- fit_surv(omicsData)
plot(sfit)
}
#' Basic survival analysis summary
#'
#' Implements overall survival analysis or progression-free survival analysis, depending upon the datatypes supplied to
#' surv_designation, and gives a summary of the results.
#'
#' @param omicsData A pmartR data object of any class, which has a `group_df` attribute that is usually created by the `group_designation()` function
#' @param percent The percentile
#' @param ... extra arguments passed to regexpr if pattern is specified
#' @return if `percent` is provided then the time at which that probability of death is returned; else, the summary of the `survival` object is returned
#'
#' @examplesIf requireNamespace("pmartRdata", quietly = TRUE)
#' \dontrun{
#' library(OvarianPepdataBP)
#' attr(tcga_ovarian_pepdata_bp, "survDF") <- list(t_death = "survival_time",
#' ind_death = "vital_status")
#' # No percent is provided so the entire object is returned
#' summary_km(tcga_ovarian_pepdata_bp)
#'
#' # Percent is provided so corresponding time point is returned
#' summary_km(tcga_ovarian_pepdata_bp, .4)
#' }
#'
#' @export
#'
summary_km <- function(omicsData, percent = NULL, ...) {
# Summary object
sfit1 <- summary(fit_surv(omicsData, ...))
if (!is.null(percent)) {
row_of_interest <- which.min(abs(sfit1$surv - percent))
return(sfit1$time[row_of_interest])
}
return(sfit1)
}