#' Create a Kaplan Meier Estimator by using a TCGA Dataset
#'
#' \code{survival_adaptable()} plots a highly modifyable survival estimator using the survminer package as a basis
#'
#' @param Eset An Expression Set
#' @param survival For overall survival = "overall", for disease free survival: survival = "DFS"
#' @param clinical a character vector with one or more factors of the clinical patient data
#' @param mutation a character vector with one or more factors of the patients mutation data added to the ExpressionSet by add_mutation()
#' @param expression a character vector with one or more gene expression data
#' @param optimal calculate the optimal cutpoint using the max_stat ranking method, will overide the value value
#' @param value Defines the value to subdivide the gene expression groups.
#' Numeric: Devided into two groups
#' "q": 25 \% quantile, 25-75 \% quantile and 75 \% quantile
#' @param gene_signature For more than one gene, how the genes for the signature are averaged, either "median" or "mean". If left out, all gene expression combinations will be plotted individually
#' @param groupings You can group factors from the clinical data
#' @param exclude Factors to exclude from any of the given clinical, expression or mutation data. The order of the to be exlcuded variables does not matter.
#' @param show_only Works only if exclude is used. If you want to exclude specific values but only show the regression for either "gene expression", "clinical" data or "mutation data"
#' @param p.val Displays the p-Value on the graph, only possible if a single factor is analyzed
#' @param logrisk Displays how the p-value was calculated
#' @param xlabel User defined x-axis label
#' @param ylabel User defined y-axis label
#' @param plot_title Title of the Plot, if not stated, no title will be shown
#' @param legend_position Where should the legend be? doesn't work with theme_bw
#' @param legend_rows How many rows are used for the legend
#' @param plot_cutpoint Plot the graph how the optimal cutpoint was calculated, normal survival plot will be REPLACED by surv_cutpoint: Determine the optimal cutpoint for each variable using ’maxstat’
#' @param risk_table If the risk table is shown or not, use FALSE or TRUE
#' @param return_df Returns the Dataframe that is used for calculatinf the Cox Regression
#' @param return_fit Returns the fit of the Regression. [[1]] = Fit, [[2]] = p-value, [[3]] = Coxph
#' @param factor_list Returns the list of factors that are in the dataframe. These factors can be used e.g. for exclusen.
#' @param ... Additional variables that can be passed over to the ggsurvplot function of the survminer package
#'
#' @return The modified survival Estimator from the survminer package to be able to handle TCGA Data
#' @export
#'
#' @import survival
#'
#'
#' @examples
#' \dontrun{
#' Survival_adaptable(expression = c("FOXA2"), Eset = Eset,
#' optimal = TRUE, p.val=TRUE,
#' xlabel="Days", legend_position = "top", average = "median",
#' plot_cutpoint=F,risk_table=TRUE)
#' }
<- function (Eset,
survival="overall",
clinical = ,
mutation = ,
= ,
optimal = ,
value,
gene_signature = "no",
groupings= ,
exclude,
show_only=,
p.val = ,
logrisk = ,
= "jco",
xlabel,
ylabel,
plot_title = "",
legend_position,
legend_rows = 1,
plot_cutpoint=,
risk_table = ,
return_df=,
return_fit=,
factor_list=,
conf_int = F,
) {
(survival)
if ((xlabel)){
xlabel <- "time"
}
if ((ylabel)){
ylabel <- "Survival"
}
if ((clinical) & (mutation) & ()) {
("You have to define at least one gene, mutation or clinical feature you want to calculate survival estimators")
}
(survival=="overall") {
<- (Biobase::pData(Eset)$X_OS)
event <- Biobase::pData(Eset)$X_OS_IND
} (survival=="disease free"){
<- (Biobase::pData(Eset)$X_DFS)
event <- Biobase::pData(Eset)$X_DFS_IND
} {
(("You cannot use \"", survival, "\" as censoring, use \"overall\" or \"DFS\""), sep="")
}
# Check if variables exist in dataframe
not_existing_exprs <-
not_existing_clinical <-
not_existing_mutation <-
(clinical[1]!=){
exist_in_pData <- clinical[clinical %in% (Biobase::pData(Eset))]
((exist_in_pData)<(clinical)){
not_existing_clinical <- clinical!clinical %in% (Biobase::pData(Eset))]
}
} {
exist_in_pData <-
}
([1]!=){
exist_in_exprs <- [expression %in% (Biobase::exprs(Eset))]
((exist_in_exprs)<()){
not_existing_exprs <- ! %in% (Biobase::exprs(Eset))]
}
} {
exist_in_exprs <-
}
(mutation[1]!=){
exist_in_mutation <- mutation[mutation %in% (Biobase::pData(Eset))]
((exist_in_mutation)<(mutation)){
not_existing_mutation <- mutation!mutation %in% (Biobase::pData(Eset))]
}
} {
exist_in_mutation <-
}
((not_existing_mutation)!=0 | (not_existing_exprs)!=0 | (not_existing_clinical)!=0){
((not_existing_mutation)!=0){
not_existing_mutation_stop <- ((not_existing_mutation,collapse = " & "), ": Not added to the ExpressionSet by add_mutation(); ",sep="")
} {
not_existing_mutation_stop <-
}
((not_existing_mutation_stop,((not_existing_exprs,not_existing_clinical),collapse = " & "), ": not in the clinical and expression data", sep=""))
}
((exist_in_exprs)==1){
expression_data <- Biobase::exprs(Eset)[exist_in_exprs,]
} {
expression_data <- (Biobase::exprs(Eset)[exist_in_exprs,])
}
z <- (, event, expression_data, Biobase::pData(Eset),exist_in_pData], Biobase::pData(Eset),exist_in_mutation])
(z) <- ("time","event",exist_in_exprs,exist_in_pData,exist_in_mutation)
# remove NA
(i 1:((z))){
z <- z[z,i]!="[Not Available]",]
z <- z[z,i]!="[Not Applicable]",]
z <- z!(z,i]),]
}
((z)[1]==0){
("No data to do further analysis, maybe your clinical data factor was empty?")
}
# Groupings
(groupings[1]!=){
grouping_factors <- "groupings"
r <-
l <-
# paste(c(exist_in_exprs,exist_in_mutation,exist_in_pData)[!groupings])
(j 3:((z))){
grouping_indicator <- (z,j],(x) x==groupings)
i <- j-2
r[i]] <- (grouping_indicator,2,)
(i>1){
(i>2){
l <- r[i]] + l
} {
l <- r[i]] + r[i-1]]
}
}
}
z$groupings <- (l>0,("Group:",(groupings,collapse=", ")),"rest")
} {
grouping_factors <-
}
if (!(exclude)){
(i 1:(exclude)){
(j 3:((z))){
z <- z!z,j]==exclude[i], ]
}
}
(((z))==0){
("All data was excluded with your exclude settings, no further analysis possible.")
}
}
((exist_in_exprs)>0){
(optimal){
(!(value)){
("Value: When Optimal cutpoint is calculated, the value parameter will be ignored")
}
value_cutpoint <- survminer::surv_cutpoint(z, ="time", event="event", variables = exist_in_exprs)
value_sum <- (value_cutpoint)
value_single <- value_sum$cutpoint
value_combined <- (value_single)
(plot_cutpoint){
::(mfrow=(2,1))
(::(value_cutpoint, exist_in_exprs, = "jco"))
}
two_groups <- T
} ((value)){
("Please let either an optimal cutpoint be calculated by optimal=TRUE, use a user defined cutpoint by value=X or divide the groups into their quantiles by value=\"q\"")
} ((value)){
(plot_cutpoint){
("plot_cutpoint: The Optimal cutpoint cannot be plotted when optimal is not TRUE, therefore plot_cutpoint is ignored")
}
((value)==(exist_in_exprs)){
value_single <- value
} ((value)==1){
(i 1:(exist_in_exprs)){
value_single[i] <- value
}
} {
("The cutpoint value needs to be either as long as genes to be analyzed: \"value = c(gene1, gene2, gene3,...)\"\nor value needs to have one value for all genes \"value = X\"")
}
value_combined <- (value)
two_groups <- T
} (value=="q") {
two_groups <- F
}
(two_groups){
((exist_in_exprs)>1){
(gene_signature=="mean" | gene_signature=="median"){
df_matrix <- ((z,exist_in_exprs]), = (exist_in_exprs), byrow=)
z$ <- Biobase::rowMedians(df_matrix)
z$ <- (df_matrix)
(gene_signature=="mean"){
(i 1:(exist_in_exprs)){
z,("values_",exist_in_exprs[i],sep="")] <- z,exist_in_exprs[i]]
}
z$gene_signature <- ((z,exist_in_exprs],1,) > value_combined,
"High Expression",
"Low Expression")
} (gene_signature=="median"){
(i 1:(exist_in_exprs)){
z,("values_",exist_in_exprs[i],sep="")] <- z,exist_in_exprs[i]]
}
z$gene_signature <- ((z,exist_in_exprs],1,::) > value_combined,
"High Expression",
"Low Expression")
}
gene_signature_factor <- "gene_signature"
} (gene_signature=="no") {
gene_signature_factor <-
} (!(gene_signature)){
("Multiple genes can only be combined using \"median\" or \"mean\" as the parameter for `gene_signature`to use them as a gene expression signature. If you do not want a gene expression signature, leave the parameter out")
} {
gene_signature_factor <-
}
} {
(!(gene_signature) & (exist_in_exprs)==1){
("gene_signature: Gene expression signature cannot be calculated if only one gene there to be analyzed, parameter average will be ignored")
}
gene_signature_factor <-
}
(i 1:(exist_in_exprs)){
z,("values_",exist_in_exprs[i],sep="")] <- z,exist_in_exprs[i]]
z,exist_in_exprs[i]] <- (z,exist_in_exprs[i]] > value_single[i],
"High Expression",
"Low Expression")
}
} {
((exist_in_exprs)>1){
quantile_values <- (z,exist_in_exprs],2,(x) ::(x,(.25, .50, .75)))
(gene_signature=="mean" | gene_signature=="median"){
df_matrix <- ((z,exist_in_exprs]), = (exist_in_exprs), byrow=)
z$ <- Biobase::rowMedians(df_matrix)
z$ <- (df_matrix)
z$gene_signature <- (z,gene_signature] < (quantile_values[1,]), "Lower Quantile", (z,gene_signature] < (quantile_values[3,]), "Intermediate", "Upper Quantile"))
gene_signature_factor <- "gene_signature"
} {
gene_signature_factor <-
}
(i 1:(exist_in_exprs)){
z,exist_in_exprs[i]] <- (z,exist_in_exprs[i]] < quantile_values[1,i], "Lower Quantile", (z,exist_in_exprs[i]] < quantile_values[3,i], "Intermediate", "Upper Quantile"))
}
} {
gene_signature_factor <-
quantile_values <- ::(z,exist_in_exprs], (.25, .50, .75))
z,exist_in_exprs] <- (z,exist_in_exprs] < quantile_values[1], "Lower Quantile", (z,exist_in_exprs] < quantile_values[3], "Intermediate", "Upper Quantile"))
}
}
} {
gene_signature_factor <-
(optimal){
("Optimal: Optimal cutpoint cannot be calculated for clinical data or mutation data, parameter optimal will be ignored")
}
(!(value)){
("Value: Clinical and mutational data cannot be divided into groups by a user defined gene expression value, parameter value and optimal will be ignored")
}
(!(gene_signature)){
("gene_signature: A gene expression signature cannot be calculated if no genes are there to be analyzed, parameter gene_signature will be ignored")
}
}
(factor_list){
(((z))>3){
return_factor_list <- (z,(exist_in_exprs,exist_in_pData,exist_in_mutation,gene_signature_factor,grouping_factors)],function (x) ((x)))
} {
return_factor_list <- ((z)[3],": ", (((z,3])),collapse = ", "), sep="")
}
(return_df){
("return_df: Please use \"factor_list = FALSE\" to return the data frame")
}
(return_fit){
("return_fit: Please use \"factor_list = FALSE\" to return the fit")
}
(return_factor_list)
}
# Factorize
(((exist_in_exprs,exist_in_pData,exist_in_mutation,gene_signature_factor,grouping_factors))>1){
z[c(exist_in_exprs,exist_in_pData,exist_in_mutation,gene_signature_factor,grouping_factors)] <- (z,(exist_in_exprs,exist_in_pData,exist_in_mutation,gene_signature_factor,grouping_factors)],)
} ((gene_signature_factor)!=0){
z,3] <- (z,3])
}
(return_df){
(z)
}
((exist_in_pData)!=0 & show_only!="expression data" & show_only!="mutation data"){
pData_formula <- (("`", exist_in_pData, "`", sep=""),collapse = " + ")
} {
pData_formula <-
}
((exist_in_exprs)!=0 & show_only!="clinical" & show_only!="mutation data") {
exprs_formula <- (("`", exist_in_exprs, "`", sep=""),collapse = " + ")
# if(length(exist_in_pData!=0)){
# exprs_formula <- paste("+",exprs_formula)
# }
} {
exprs_formula <-
}
((exist_in_mutation)!=0 & show_only!="clinical" & show_only!="expression data" ){
mutation_formula <- (("`", exist_in_mutation, "`", sep=""),collapse = " + ")
# if(length(exist_in_exprs!=0)){
# mutation_formula <- paste("+",mutation_formula)
# }
} {
mutation_formula <-
}
# Calculate survfit
((gene_signature_factor)>0 & (exist_in_exprs)>1){
(groupings[1]!=){
<- "groupings + gene_signature"
} {
<- ((pData_formula,mutation_formula,"gene_signature"),collapse = " + ")
}
} {
(groupings[1]!=){
<- (("groupings",exprs_formula),collapse=" + ")
} {
<- ((pData_formula,mutation_formula,exprs_formula),collapse=" + ")
}
}
# p-value:
(((z))==3){
z,3] <- (z,3])
} {
z,3:((z))] <- (z,3:((z))],)
}
((exist_in_exprs)>0){
count_exps <- (z,exist_in_exprs])/(exist_in_exprs)
count_mut <- (z,exist_in_mutation])*(exist_in_mutation)
count_clin <- (z,exist_in_pData])*(exist_in_pData)
(count_exps==0) count_exps=1
(count_mut==0) count_mut=1
(count_clin==0) count_clin=1
((count_exps*count_mut*count_clin)<3){
(p.val == ){
p.val =
} {
p.val =
}
} {
p.val =
}
}
Survobject <- ::(("Surv(time = time, event = event) ~", ))
fit <- (survival::survfit,
( = Survobject,
= z,
="kaplan-meier",
conf.type="log"
)
)
(return_fit){
cox <- (survival::coxph,
( = Survobject,
= z)
)
sdf <- (survival::survdiff,
( = Survobject,
= z)
)
p.val <- 1 - ::(sdf$chisq, (sdf$n) - 1)
quantile_survival <- ::(fit, probs = (0.05,0.25,0.5,0.75,0.95))
((fit,cox,("P-Value =",p.val),quantile_survival))
}
# Calculate scaling of X-Axis
((z$,na.rm=)>5000){
breaks <- 2000
} ((z$,na.rm=)>500) {
breaks <- 500
} {
breaks <- 50
}
# if(length(colnames(z))>3){
# naming_factors <- sapply(z[,c(exist_in_exprs,exist_in_pData,exist_in_mutation,combined_factor)],function (x) levels(as.factor(x)))
# } else {
# naming_factors <- levels(factor(z[,3]))
# }
((gene_signature_factor)!=0){
expression_title <- ("gene signature of",(exist_in_exprs,collapse=" & "))
} {
expression_title <- exist_in_exprs
}
legend_title <- ((expression_title,exist_in_pData,exist_in_mutation),collapse=" & ")
p <- survminer::ggsurvplot(fit,
=z,
= plot_title,
surv.scale = ("percent"),
= ,
#xscale = "d_y",
= legend_position,
#legend.labs = naming_factors,
legend.title = legend_title,
# xlim = c(0,2000), # present narrower X axis, but not affect survival estimates.
break.time.by = breaks, # break X axis in time intervals by 500.
pval = p.val, # show p-value of log-rank test.
pval.method = ,
log.rank.weights = "n",
conf.int = conf_int, # show confidence intervals for point estimaes of survival curves.
xlab = xlabel,
ylab = ylabel,
risk.table = risk_table, # show risk table.
# risk.table.y.text.col = T, # colour risk table text annotations.
# risk.table.y.text = F, # show bars instead of names in text annotations in legend of risk table
#ncensor.plot = TRUE
#cumevents=T
#cumcensor=T
# tables.theme = theme_cleantable(),
# ggtheme = theme_bw() # Change ggplot2 theme
surv.median.line="hv",
)
p <- p + ggplot2::guides(colour = ggplot2::guide_legend( = legend_rows))
p
}
