I’ve been given a data set and code to go off for learning R, but after making some adjustments I can’t generate a PDF of my graphs even though no error messages are appearing in my console. I’m also getting a “+” sign before every line in my console, I found out you can click ‘esc’ on your keyboard and it fixes the issue, but I’m still not getting a PDF – I’m not sure if these issues are connected or not…
I know this is a bunch of code, but my non-coding brain thinks there could be a function I could include to fix this?
If someone can help I’d really appreciate it!
install.packages(c("BiocManager","pacman"))
install.packages("ggrastr","rms","viridis","ggplot2")
BiocManager::install(c("Biobase","survival","survcomp","survminer","impute"))
library(Biobase)
library(viridis)
library(ggrastr)
library(rms)
library(survival)
library(survminer)
library(impute)
workdir <- "~/Desktop/Honours/Data Analysis/SOMS4101_Material"
gitrepo <- paste0(workdir,"/chromatinAnthra-master")
download.file(url = "https://github.com/cancersysbio/chromatinAnthra/archive/refs/heads/master.zip", destfile = "master.zip")
unzip("master.zip")
rescale <- function(x, na.rm=FALSE, q=0.05) {
ma <- quantile(x, probs=1-(q/2), na.rm=na.rm)
mi <- quantile(x, probs=q/2, na.rm=na.rm)
x <- (x - mi) / (ma - mi)
return((x - 0.5) * 2)
}
mulvarrfs.res2<-mulvarrfs.res
probes= mulvarrfs.res2$probe[which(mulvarrfs.res2$gene=="KDM4B")]
X.b = cbind(pheno.df[,setdiff(colnames(pheno.df),"tmain")],exp_c=rescale(exprs(mergedCOMBAT)[probes[1],rownames(pheno.df) ]))
d <- datadist(X.b)
#options(datadist='dd')
assign('dd', datadist(X.b))
options(datadist="dd")
label = "anthra_vs_nonAnthra_v2"
plot_geneInfo = function(mergedCOMBAT,mulvarrfs.res2,pheno.df,genename="KDM4B",label){
#KDM4B
pdf(paste0("surv_",label,"_",genename,".pdf"),width = 8,height = 7)
#i=which(mergedCOMBAT@featureData@data$`Gene Symbol`==genename)
#probes_i = mergedCOMBAT@featureData@data$ID[i]
probes= mulvarrfs.res2$probe[which(mulvarrfs.res2$gene==genename)]
X.b = cbind(pheno.df[,setdiff(colnames(pheno.df),"tmain")],exp_c=rescale(exprs(mergedCOMBAT)[probes[1],rownames(pheno.df) ]))
d <- datadist(X.b)
#options(datadist='dd')
assign('dd', datadist(X.b))
options(datadist="dd")
cxmodel2 = cph(Surv(rfs.t,rfs.e)~exp_c*anthra+rcs(age,3)+er+pr+her2+lymphNodePos+t.stage+cohort+MKI67_4,data = X.b,x = T,y=T,surv = T)
print(cxmodel2)
print(anova(cxmodel2))
min_lim = median(X.b$exp_c)-sd(X.b$exp_c)
max_lim = median(X.b$exp_c)+sd(X.b$exp_c)
# low antrha green / low non anthra violet
# high anthra yellow / high non anthra blue
#################################### MODIFY THIS! ###########################################################
# PLOT 1
survplot(cxmodel2,
anthra,
exp_c=min_lim,
er=F,
pr=F,
her2=F,
conf.int = T,
adj.subtitle=F,
col=c("orange", "magenta"),
xlim=c(0,13), ##extent of X-Axis
lwd=3,
lty=c(1,2),
label.curves=F, #not changed
n.risk = T,
y.n.risk = 0.6,
time.inc=2,
col.fill = c("#0000FF33", "#FFA50033"),
ylab = ("Surivival Probability"),
xlab = "follow up time (Days)",
# PLOT 2
survplot(cxmodel2,
anthra,
exp_c=max_lim,
er=F,
pr=FALSE,
her2=F,
conf.int = T,
adj.subtitle=F,
col=viridis_pal()(4)[c(2,4)],
xlim=c(0,13),
lwd=3,
lty=1,
label.curves=F),
#X.b$exp_chemo = NA,
new_exp = survminer::surv_categorize( survminer::surv_cutpoint(X.b,time="rfs.t",event="rfs.e",variables = "exp_c")),
# PLOT 3
plot(survminer::surv_cutpoint(X.b,time="rfs.t",event="rfs.e",variables = "exp_c")),
X.b$exp_chemo[new_exp$exp_c=="low" & X.b$anthra=="TRUE"]<-"Low expression and anthracycline",
X.b$exp_chemo[new_exp$exp_c=="low" & X.b$anthra=="FALSE"]<-"Low expression and non anthracycline",
X.b$exp_chemo[new_exp$exp_c=="high" & X.b$anthra=="TRUE"]<-"High expression and anthracycline",
X.b$exp_chemo[new_exp$exp_c=="high" & X.b$anthra=="FALSE"]<-"High expression and non anthracycline",
# PLOT 4
survcomp::km.coxph.plot(Surv(rfs.t,rfs.e)~exp_chemo,data.s = X.b,
x.label="time (years)",y.label="OS",main.title="",show.n.risk = T,
n.risk.step=3,leg.text = levels(factor(X.b$exp_chemo)),
.col= viridis_pal()(4)[c(4,2,3,1)],xlim=c(0,13),.lwd=3),
X.b$exp_chemo_2 <- X.b$exp_chemo,
X.b$exp_chemo_2[X.b$exp_chemo %in% c("High expression and anthracycline","High expression and non anthracycline")]<-NA,
# PLOT 5
survcomp::km.coxph.plot(Surv(rfs.t,rfs.e)~exp_chemo_2,
data.s = X.b,
x.label="time (years)",
y.label="OS",main.title=" ",
show.n.risk = T,
n.risk.step=3,
leg.text = levels(factor(X.b$exp_chemo_2)),
.col= viridis_pal()(4)[c(3,1)],
xlim=c(0,13),
.lwd=3),
X.b$exp_chemo_3 <- X.b$exp_chemo,
X.b$exp_chemo_3[!X.b$exp_chemo %in% c("High expression and anthracycline","High expression and non anthracycline")]<-NA,
survcomp::km.coxph.plot(Surv(rfs.t,rfs.e)~exp_chemo_3,data.s = X.b,
x.label="time (years)",y.label="OS",main.title=" ",show.n.risk = T,
n.risk.step=3,leg.text = levels(factor(X.b$exp_chemo_3)),
.col= viridis_pal()(4)[c(4,2)],xlim=c(0,13),.lwd=3),
#table(X.b$exp_chemo_3,X.b$cohort)
#table(X.b$exp_chemo_3,X.b$X20)
cxmodel3 = cph(Surv(rfs.t,rfs.e)~exp_c*anthra+rcs(age,3)+er+pr+her2+lymphNodePos+t.stage+cohort+MKI67_4,data = X.b,x = T,y=T,surv = T),
hazards3 = predict(cxmodel2,X.b,se.fit = T,center.terms=T),
p = ggplot(Predict(cxmodel3,
exp_c,
anthra,
er=T,
pr=F,
her2=F,
t.stage=2,
lymphNodePos=F,
cohort="KAO")) + theme_classic() +
scale_color_viridis(discrete = T) +
geom_point(data =X.b,aes(x=exp_c,y=hazards3$linear.predictors)) +
xlim(dd$limits["Low:prediction","exp_c"],dd$limits["High:prediction","exp_c"]),
print(p),
dev.off(),
plot_geneInfo(mergedCOMBAT,mulvarrfs.res2,pheno.df,"KDM4B","anthra_vs_nonAnthra_v2"),
plot_geneInfo(mergedCOMBAT,mulvarrfs.res2,pheno.df,"KAT6B","anthra_vs_nonAnthra_v2"),
plot_geneInfo(mergedCOMBAT,mulvarrfs.res2,pheno.df,"BCL11A","anthra_vs_nonAnthra_v2"),
New contributor
Jooya Kalantar is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.