0
|
1 # author: Pauline Ribeyre
|
|
2
|
|
3
|
|
4 #############
|
|
5 # find ions #
|
|
6 #############
|
|
7
|
|
8 is_ion_present <- function(GC_results, rt, mz_list) {
|
|
9 # Looks for a ion in the pseudospectra.
|
|
10 #
|
|
11 # Args:
|
|
12 # GC_results: result of a test by runGC().
|
|
13 # rt: retention time of the ion to look for.
|
|
14 # mz_list: list of the ion's isotopes' m/z.
|
|
15 #
|
|
16 # Returns:
|
|
17 # The number of times the ion was found in the pseudospectra.
|
|
18
|
|
19 delta_rt <- ions_delta_rt
|
|
20 delta_mz <- ions_delta_mz
|
|
21
|
|
22 res = 0
|
|
23
|
|
24 pseudospectra <- GC_results$PseudoSpectra
|
|
25 for (peak in pseudospectra) { # check all the peaks
|
|
26 if (peak$rt >= rt - delta_rt && peak$rt <= rt + delta_rt) { # check if the retention time matches
|
|
27
|
|
28 ps <- peak$pspectrum
|
|
29 found_ion <- rep(NA, length(mz_list))
|
|
30
|
|
31 for (i in 1:length(mz_list)) { # check the presence of all the ions given in parameter
|
|
32
|
|
33 ion <- mz_list[i][[1]]
|
|
34 found_ion[i] <- sum((ps[,"mz"] >= ion - delta_mz) & (ps[,"mz"] <= ion + delta_mz))
|
|
35
|
|
36 } # end for
|
|
37
|
|
38 if (length(which(found_ion > 1)) > 0) {
|
|
39 stop("The ion was found more than 1 time in the same pseudospectrum.")
|
|
40 }
|
|
41
|
|
42 if (sum(found_ion) == length(mz_list)) {
|
|
43 res <- res + 1
|
|
44 }
|
|
45
|
|
46 } # end if
|
|
47 } # end for
|
|
48
|
|
49 return (res)
|
|
50
|
|
51 }
|
|
52
|
|
53
|
|
54 ##############
|
|
55 # duplicates #
|
|
56 ##############
|
|
57
|
|
58 count_duplicates_function <- function(GC_results) {
|
|
59 # Counts the number of ions that have been found in more than one pseudospectra (slow).
|
|
60 #
|
|
61 # Args:
|
|
62 # GC_results: result of a test by runGC().
|
|
63 #
|
|
64 # Returns:
|
|
65 # The number of duplicates.
|
|
66
|
|
67 delta_rt <- duplicates_delta_rt
|
|
68 delta_mz <- duplicates_delta_mz
|
|
69
|
|
70 pseudospectra <- GC_results$PseudoSpectra
|
|
71
|
|
72 # order by pseudospectrum rt
|
|
73 rts <- sapply(pseudospectra, function(x) x$rt)
|
|
74 rts_order <- order(rts)
|
|
75
|
|
76 duplicates <- NULL # record all the duplicates we find
|
|
77 already_counted <- NULL # record which duplicates have already been counted
|
|
78
|
|
79 for (pspec in 1:length(rts)) { # for each pseudospectrum
|
|
80
|
|
81 i <- rts_order[pspec]
|
|
82 this_rt <- rts[i]
|
|
83
|
|
84 # number of following pseudospectra to consider
|
|
85 k <- 0
|
|
86 while (pspec + k + 1 < length(rts_order) && rts[rts_order[pspec + k + 1]] < this_rt + delta_rt) { # check if the retention time matches
|
|
87 k <- k + 1
|
|
88 }
|
|
89
|
|
90 mzs <- lapply(pseudospectra[rts_order[pspec:(pspec + k)]], function(x) x$pspectrum)
|
|
91 mzs_sorted <- NULL
|
|
92 for (j in 1:(k + 1)) {
|
|
93 mzs_sorted <- rbind(mzs_sorted, t(sapply(1:length(mzs[[j]][,1]), function(n) c(rts[rts_order[pspec + j - 1]], mzs[[j]][,1][n]))))
|
|
94 }
|
|
95 mzs_sorted <- mzs_sorted[order(mzs_sorted[,2]),] # column 1: rt ; column 2: m/z
|
|
96
|
|
97 i <- 1
|
|
98 while (i <= nrow(mzs_sorted)) { # for each m/z
|
|
99
|
|
100 j <- 1
|
|
101 this_already_counted <- nrow(merge(data.frame(t(mzs_sorted[i,])), data.frame(already_counted))) > 0
|
|
102
|
|
103 if (is.null(already_counted) || !this_already_counted) {
|
|
104 this_mz <- mzs_sorted[i,2]
|
|
105 already_counted <- rbind(already_counted, mzs_sorted[i,])
|
|
106 same_ion <- 1 # count how many times this m/z appears
|
|
107
|
|
108 while (j <= nrow(mzs_sorted) - i && mzs_sorted[i+j,2] >= this_mz - delta_mz && mzs_sorted[i+j,2] <= this_mz + delta_mz) { # check if the m/z matches
|
|
109 this_already_counted <- nrow(merge(data.frame(t(mzs_sorted[i+j,])), data.frame(already_counted))) > 0
|
|
110
|
|
111 if (!this_already_counted) { # do not count a duplicate more than once
|
|
112 same_ion <- same_ion + 1
|
|
113 already_counted <- rbind(already_counted, mzs_sorted[i+j,]) # record which duplicates have already been counted
|
|
114 this_mz <- mzs_sorted[i+j,2]
|
|
115 }
|
|
116
|
|
117 j <- j + 1
|
|
118 }
|
|
119
|
|
120 # first/min m/z, last/max m/z, rt where we started, number of duplicates
|
|
121 if (same_ion > 1)
|
|
122 duplicates <- rbind(duplicates, c(mzs_sorted[i,2], this_mz, mzs_sorted[i,1], same_ion))
|
|
123 }
|
|
124 i <- i + j
|
|
125 }
|
|
126
|
|
127 }
|
|
128
|
|
129 if (is.null(duplicates))
|
|
130 duplicates <- c(0, 0, 0, 0)
|
|
131
|
|
132 return (duplicates)
|
|
133
|
|
134 }
|
|
135
|