Created
December 7, 2015 22:24
-
-
Save rhliang/29a43b48bc2881d002a2 to your computer and use it in GitHub Desktop.
Simulates primer ID coverage given the number of molecules and the number of pIDs available.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## Script that simulates the labelling of RNA molecules by pIDs. | |
| simulate.labelling <- function(num.trials, num.molecules, num.pIDs) | |
| { | |
| labelling.trials <- | |
| sapply(1:num.trials, | |
| function (x) | |
| { | |
| pIDs.assigned <- sample(1:num.pIDs, num.molecules, replace=TRUE) | |
| num.distinct.pIDs <- length(unique(pIDs.assigned)) | |
| return(num.distinct.pIDs) | |
| }) | |
| } | |
| labellings <- list() | |
| num.trials <- 1000 | |
| num.RNA <- c(1000, 1995, 5012, 10000, 19952, 50119, 100000, 199527) | |
| num.pIDs <- c(36196, 60407, 129057, 237390, 447132, 1066811, 2074580, 4064816) | |
| for (i in 1:length(num.RNA)) | |
| { | |
| cat("Simulating ", num.trials, " trials with ", num.RNA[i], | |
| " RNA molecules and ", num.pIDs[i], " pIDs....\n", sep="") | |
| labellings[[as.character(num.RNA[i])]] <- | |
| simulate.labelling(num.trials, num.RNA[i], num.pIDs[i]); | |
| ## Get the level of coverage which we saw in 90% of the trials, i.e. | |
| ## the coverage achieved with 0.9 probability. | |
| pIDs.used <- quantile(labellings[[as.character(num.RNA[i])]], 0.1) | |
| cat("Unique pIDs used (coverage achieved) in 90% of trials: ", | |
| pIDs.used, " (", pIDs.used/num.RNA[i], ")\n", sep="") | |
| } | |
| save.image(file=paste("Simulation ", date(), ".RData", sep="")) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment