Skip to content

Instantly share code, notes, and snippets.

@rhliang
Last active December 7, 2015 22:27
Show Gist options
  • Select an option

  • Save rhliang/111518a44dc83dbcafeb to your computer and use it in GitHub Desktop.

Select an option

Save rhliang/111518a44dc83dbcafeb to your computer and use it in GitHub Desktop.
Produce histograms to summarize the simulation results produced by pIDLabellingSimulation.r.
## This code assumes you've run pIDLabellingSimulation.r already. Either
## run this immediately after, or load the file created by it, e.g.
## load([file created by pIDLabellingSimulation.r])
symbol.scale <- 1;
axis.scale <- 1;
label.scale <- 1;
text.size <- 0.79;
for (i in 1:length(num.RNA))
{
cat("Outputting results for ", num.trials, " trials with ", num.RNA[i],
" RNA molecules and ", num.pIDs[i], " pIDs....\n", sep="")
curr.num.RNA <- num.RNA[i]
curr.num.pIDs <- num.pIDs[i]
curr.num.RNA.string <- as.character(curr.num.RNA)
## Make a histogram of the results for this number of RNA.
pdf(paste("Labelling", curr.num.RNA.string, "molecules.pdf", sep=""))
curr.data <- labellings[[as.character(curr.num.RNA)]]
curr.hist <- hist(curr.data, main=paste(curr.num.RNA.string,
"RNA molecules,", curr.num.pIDs, "pIDs"),
xlab="unique pIDs used",
cex.axis=axis.scale,
cex.lab=label.scale)
pIDs.used <- quantile(labellings[[as.character(num.RNA[i])]], 0.1)
abline(v=pIDs.used, lty="dashed")
text(pIDs.used,
max(curr.hist$counts)*0.6,
labels=expression(paste("pIDs used in" >= 90, "% of trials")),
pos=2,
cex=text.size)
text(pIDs.used,
max(curr.hist$counts)*0.57,
labels=paste("(", round(100*pIDs.used/num.RNA[i], digits=2),
"%-good coverage)", sep=""),
pos=2,
cex=text.size)
dev.off()
cat("Unique pIDs used (coverage achieved) in 90% of trials: ",
pIDs.used, " (", pIDs.used/num.RNA[i], ")\n", sep="")
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment