### R code from vignette source 'vignettes/DTA/inst/doc/DTA.Rnw' ################################################### ### code chunk number 1: Loading library ################################################### library(DTA) ################################################### ### code chunk number 2: Path to RData files ################################################### dataPath = system.file("data", package="DTA") ################################################### ### code chunk number 3: Loading datamat ################################################### load(file.path(dataPath, "datamat.RData")) colnames(datamat)[1:6] rownames(datamat)[1:6] datamat[1:6,1:3] ################################################### ### code chunk number 4: Loading phenomat ################################################### load(file.path(dataPath, "phenomat.RData")) head(phenomat) ################################################### ### code chunk number 5: Loading tnumber ################################################### load(file.path(dataPath, "tnumber.RData")) head(tnumber) ################################################### ### code chunk number 6: Total expression histogram ################################################### Totals = datamat[,which(phenomat[,"fraction"]=="T")] Total = apply(log(Totals),1,median) plotsfkt = function(){ hist(Total,breaks = seq(0,ceiling(max(Total)),1/4), cex.main=1.25,cex.lab=1,cex.axis=1, main="Histogram of log(Total)", xlab="gene-wise median of total samples") hist(Total[Total >= 5],breaks = seq(0,ceiling(max(Total)),1/4), col = "#08306B",add = TRUE) } DTA.plot.it(filename = "histogram_cut_off",plotsfkt = plotsfkt,saveit = TRUE,notinR = TRUE) ################################################### ### code chunk number 7: Loading reliable ################################################### load(file.path(dataPath, "reliable.RData")) head(reliable) ################################################### ### code chunk number 8: Estimation procedure ################################################### res = DTA.estimate(phenomat,datamat,tnumber,ccl = 150, mRNAs = 60000,reliable = reliable, correctedlabeling = TRUE,condition = "real_data",plots = TRUE,notinR = TRUE,folder = ".") ################################################### ### code chunk number 9: Entries of res ################################################### names(res) ################################################### ### code chunk number 10: Entries of res6 ################################################### names(res[["6"]]) ################################################### ### code chunk number 11: plot half-lives/sythesis rates ################################################### load(file.path(dataPath, "orfs_ribig.RData")) load(file.path(dataPath, "orfs_rpg.RData")) load(file.path(dataPath, "orfs_tf.RData")) load(file.path(dataPath, "orfs_stress.RData")) ################################################### ### code chunk number 12: plot half-lives/sythesis rates ################################################### plotsfkt = function(){ x=res[["6"]][["hl"]][c(reliable,orfs_rpg)] y=res[["6"]][["Rsr"]][c(reliable,orfs_rpg)] ellipsescatter(x,y, groups=list("Stress"=orfs_stress,"RiBi-genes"=orfs_ribig,"Rp-genes"=orfs_rpg,"TF"=orfs_tf), colors=c("darkred","darkgreen","darkblue","grey20"), xlim=c(0,150),ylim=c(0,600),xlab="half-lives",ylab="synthesis rates", cex.main=1.25,cex.lab=1,cex.axis=1)} DTA.plot.it(filename = "ellipse_scatter",plotsfkt = plotsfkt,saveit = TRUE,notinR = TRUE) ################################################### ### code chunk number 13: Estimation procedure (no bias correction) ################################################### res.nobias = DTA.estimate(phenomat,datamat,tnumber,ccl = 150, mRNAs = 60000,reliable = reliable,plots = TRUE,notinR = TRUE,folder = ".",bicor = FALSE,condition="no_bias_correction") ################################################### ### code chunk number 14: Creating simulation object ################################################### sim.object = DTA.generate(timepoints=rep(c(6,12),2)) ################################################### ### code chunk number 15: Entries of sim.object ################################################### names(sim.object) ################################################### ### code chunk number 16: Estimation procedure for simulated data ################################################### res.sim = DTA.estimate(ratiomethod = "bias",plots = TRUE,notinR = TRUE, folder = ".",simulation = TRUE,sim.object = sim.object, condition = "simulation") ################################################### ### code chunk number 17: Estimation procedure for simulated data ################################################### res.sim = DTA.estimate(ratiomethod = "bias",plots = TRUE,notinR = TRUE, folder = ".",simulation = TRUE,sim.object = sim.object,labeling = FALSE, rankpairs = FALSE,assessment = FALSE,correlation = FALSE,bicor = FALSE, condition = "simulation_no_bias_correction") ################################################### ### code chunk number 18: Loading tDTAdata ################################################### load(file.path(dataPath, "tDTAdatamat.RData")) load(file.path(dataPath, "tDTAphenomat.RData")) load(file.path(dataPath, "tDTAtnumber.RData")) load(file.path(dataPath, "tDTAreliable.RData")) ################################################### ### code chunk number 19: Estimation procedure ################################################### timecourse.res = DTA.dynamic.estimate(tDTAphenomat,tDTAdatamat,tDTAtnumber,ccl = 150, mRNAs = 60000,reliable = tDTAreliable,LtoTratio = rep(0.1,7),plots = TRUE,notinR = TRUE, folder = ".",condition = "timecourse") ################################################### ### code chunk number 20: Entries of timecourse.res ################################################### names(timecourse.res) ################################################### ### code chunk number 21: Entries of timecourse.res1 ################################################### names(timecourse.res[["1"]]) ################################################### ### code chunk number 22: Creating timecourse simulation ################################################### nrgenes = 1000 truesynthesisrates = rf(nrgenes,5,5)*18 steady = rep(1,nrgenes) shock = 1/pmax(rnorm(nrgenes,mean = 8,sd = 4),1) induction = pmax(rnorm(nrgenes,mean = 8,sd = 4),1) changes.mat = cbind(steady,shock,shock*induction) mu.values.mat = changes.mat*truesynthesisrates mu.breaks.mat = cbind(rep(12,nrgenes),rep(18,nrgenes)) truehalflives = rf(nrgenes,15,15)*12 truelambdas = log(2)/truehalflives changes.mat = cbind(steady,shock,shock*induction,steady) lambda.values.mat = changes.mat*truelambdas lambda.breaks.mat = cbind(rep(12,nrgenes),rep(18,nrgenes),rep(27,nrgenes)) ################################################### ### code chunk number 23: Creating timecourse simulation object ################################################### timecourse.sim.object = DTA.dynamic.generate(duration = 36,lab.duration = 6, nrgenes = nrgenes,mu.values.mat = mu.values.mat,mu.breaks.mat = mu.breaks.mat, lambda.values.mat = lambda.values.mat,lambda.breaks.mat = lambda.breaks.mat) ################################################### ### code chunk number 24: Entries of timecourse.sim.object ################################################### names(timecourse.sim.object) ################################################### ### code chunk number 25: Estimation procedure for simulated timecourse data ################################################### timecourse.res.sim = DTA.dynamic.estimate(plots = TRUE,notinR = TRUE,simulation = TRUE, sim.object = timecourse.sim.object,ratiomethod = "tls",folder = ".",condition = "simulation_timecourse", regression = FALSE,labeling = FALSE,rankpairs = FALSE,assessment = FALSE,correlation = FALSE,folds = FALSE) ################################################### ### code chunk number 26: sessionInfo ################################################### toLatex(sessionInfo())