A short script to use the nontarget pattern.match function to print useful information about clusters.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

94 lines
3.1 KiB

4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
  1. library("nontarget")
  2. library("purrr")
  3. library("funprog")
  4. library("enviPat")
  5. library("stringr")
  6. library("parallel")
  7. # Configurations #############################################################
  8. # file : decides which file to read in data from #############################
  9. file <- "/path/to/file"
  10. search_isos <- c("13C", "37Cl")
  11. # Minimum size of a cluster
  12. min_cluster_size <- 3
  13. # Number of cores to be used (will be adjusted if not possible)
  14. used_cores <- 3
  15. ##############################################################################
  16. # Read in the Table ##########################################################
  17. table <- read.table(file, header=TRUE, sep=",")
  18. # Organize the tables by number ##############################################
  19. fragments <- max(table[,"Spectra_Number"])
  20. # The algorithm below guarantees linear complexity of looking up data points.#
  21. # set use_cores to safe amount
  22. use_cores = min(use_cores, detectCores()-1)
  23. # minint is the lower bound of the interval. #################################
  24. # maxint is the upper bound of the interval that contains the fragment. ######
  25. # usable memoizes as to whether or not each fragment number exists. #####
  26. minint <- unlist(map(1:fragments, function(x) 0))
  27. maxint <- unlist(map(1:fragments, function(x) 0))
  28. usable <- unlist(map(1:fragments, function(x) FALSE))
  29. # Set all of the intervals ###################################################
  30. for (i in seq(1, nrow(table))) {
  31. fragment <- table[i, "Spectra_Number"]
  32. if (! usable[fragment]) {
  33. minint[fragment] <- i
  34. }
  35. maxint[fragment] <- max(maxint[fragment], i)
  36. usable[fragment] <- TRUE
  37. }
  38. for (i in 1:fragments) {
  39. if (maxint[i] - minint[i] + 1 < min_cluster_size)
  40. usable[i] = FALSE
  41. }
  42. getdata <- function(fragment, key) {
  43. if (! usable[fragment]) {
  44. stop(str_interp("Fragment $[d]{fragment} does not exist in the data set",
  45. list(fragment=fragment)))
  46. }
  47. return(table[minint[fragment]:maxint[fragment],key])
  48. }
  49. # Add all data frames as necessary for evaluation. ############################
  50. getdataframe <- function (fragment) {
  51. mz <- getdata(fragment, "mz")
  52. time <- getdata(fragment, "time")
  53. Intensity <- getdata(fragment, "Intensity")
  54. return(data.frame(mz=mz, time=time, Intensity=Intensity))
  55. }
  56. # Incomplete: cannot get pattern.search function to work per cluster. #########
  57. # Incomplete: Get the diagnostics for a cluster and turn them into a portable#
  58. # format. (such as through tidyverse) ########################################
  59. isos <- make.isos(isotopes,use_isotopes=search_isos,
  60. use_charges=rep(1, length(search_isos)))
  61. diagnostics <- function(fragment) {
  62. points <- getdataframe(fragment)
  63. ptrn <- pattern.search(points, isos)
  64. return(ptrn)
  65. }
  66. use <- Filter(function(x) usable[x], 1:fragments)
  67. results <- mclapply(use, diagnostics, mc.cores=use_cores)
  68. # Incomplete: Make the analysis more resiliant to different sorting.
  69. # Incomplete: How to process results. (is it actually supposed to be all negatives?)
  70. ##############################################################################