Browse code

Added first support for elemental filters

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/Rdisop@70786 bc3139a8-67e5-0310-9ffc-ced21a209358

s.neumann authored on 26/10/2012 11:32:08
Showing 6 changed files

... ...
@@ -1,12 +1,16 @@
1
+2012-10-26  Steffen Neumann  <sneumann@ipb-halle.de>
2
+
3
+	* src/disop.cpp (decomposeIsotopes): Add support for element count filtering
4
+
1 5
 2012-04-21  Steffen Neumann  <sneumann@ipb-halle.de>
2 6
 
3 7
 	* R/elements.R (initializeCharges): Fixed electron mass
4 8
 
5 9
 2012-03-28  Steffen Neumann  <sneumann@ipb-halle.de>
6 10
 	* Added citations in package and documentation
7
-	
8
-2012-03-26  Steffen Neumann  <sneumann@ipb-halle.de>	
9
-	* src/disop.cpp: fixed missing UNPROTECT(), avoids warning 
11
+
12
+2012-03-26  Steffen Neumann  <sneumann@ipb-halle.de>
13
+	* src/disop.cpp: fixed missing UNPROTECT(), avoids warning
10 14
 	and possibly memory corruprtion and crashes
11 15
 
12 16
 2012-03-11  Steffen Neumann  <sneumann@ipb-halle.de>
... ...
@@ -21,9 +25,9 @@
21 25
 2012-02-23  Steffen Neumann  <sneumann@ipb-halle.de>
22 26
 	* Moved from the ancient embedded Rcpp to RcppClassic,
23 27
   	  thanks to the work of Dirk Eddelbuettel who did the porting
24
-	
28
+
25 29
 2010-11-04  Steffen Neumann  <sneumann@ipb-halle.de>
26
-	* Corrected bug that leads to wrong monoisotopic masses for molecules 
30
+	* Corrected bug that leads to wrong monoisotopic masses for molecules
27 31
 	  containing elements where the most abundant isotope is not the first one,
28 32
 	  discovered by Ralf Tautenhahn
29 33
 
... ...
@@ -51,10 +55,10 @@
51 55
 	can be useful for getMolecule("H3O+", elements=c(initializeCHNOPS(),initializeCharges()))
52 56
 	* R/Rdisop.R: fixed important bug which always initialised the
53 57
 	full PSE when supplying a special elements list
54
-	
58
+
55 59
 2008-04-23  sneumann  <sneumann@ipb-halle.de>
56 60
 
57
-	* R/elements.R: fixed some masses/abundances, 
61
+	* R/elements.R: fixed some masses/abundances,
58 62
 	                added more elements
59 63
 
60 64
 2008-02-15  sneumann  <sneumann@kons.ipb-sub.ipb-halle.de>
... ...
@@ -83,7 +87,7 @@
83 87
 
84 88
 2007-10-14  Steffen Neumann  <sneumann@ipb-halle.de>
85 89
 
86
-	* man: fixed some manpages 
90
+	* man: fixed some manpages
87 91
 
88 92
 2007-07-26  sneumann  <sneumann@ipb-halle.de>
89 93
 
... ...
@@ -101,7 +105,7 @@
101 105
 
102 106
 2007-07-03  sneumann  <sneumann@ipb-halle.de>
103 107
 
104
-	* src/disop.cpp: Clear error Message before running 
108
+	* src/disop.cpp: Clear error Message before running
105 109
 
106 110
 2007-06-14  sneumann  <sneumann@ipb-halle.de>
107 111
 
... ...
@@ -115,7 +119,7 @@
115 119
 2007-06-05  sneumann  <sneumann@ipb-halle.de>
116 120
 
117 121
 	* src/disop.cpp: Fixed stack imbalance for empty decompositions
118
-	* src/disop.cpp: improved Error Handling 
122
+	* src/disop.cpp: improved Error Handling
119 123
 
120 124
 2007-05-25  sneumann  <sneumann@ipb-halle.de>
121 125
 
... ...
@@ -125,7 +129,7 @@
125 129
 
126 130
 2007-05-15  Steffen Neumann  <sneumann@ipb-halle.de>
127 131
 
128
-	* Added addMolecule and subMolecule 
132
+	* Added addMolecule and subMolecule
129 133
 	to do arithmetics with adducts / fragments
130 134
 
131 135
 2007-05-04  Steffen Neumann  <sneumann@ipb-halle.de>
... ...
@@ -134,7 +138,7 @@
134 138
 
135 139
 2007-04-30  Steffen Neumann  <sneumann@ipb-halle.de>
136 140
 
137
-	* R/zzz.R: removed unnecesary library.dynam() 
141
+	* R/zzz.R: removed unnecesary library.dynam()
138 142
 	which also caused problems on systems without libR.so
139 143
 	* src/ims/*: moved imslib sources into src/ subfolder
140 144
 	instead of using unchanged libims.a
... ...
@@ -1,7 +1,7 @@
1 1
 Package: Rdisop
2 2
 Title: Decomposition of Isotopic Patterns
3
-Version: 1.19.0
4
-Date: 2012-04-21
3
+Version: 1.19.1
4
+Date: 2012-10-26
5 5
 Author: Anton Pervukhin <apervukh@minet.uni-jena.de>, Steffen Neumann <sneumann@ipb-halle.de> 
6 6
 Maintainer: Steffen Neumann <sneumann@ipb-halle.de>
7 7
 Description: Identification of metabolites using high precision mass
... ...
@@ -130,13 +130,16 @@ subMolecules <- function(formula1, formula2,
130 130
 
131 131
 
132 132
 decomposeMass <- function(mass, ppm=2.0, mzabs=0.0001,
133
-                          elements=NULL, filter=NULL, z=0, maxisotopes=10) {
133
+                          elements=NULL, filter=NULL, z=0, maxisotopes=10,
134
+                          minElements="C0", maxElements="C999999") {
134 135
     decomposeIsotopes(c(mass), c(1), ppm=ppm, mzabs=mzabs,
135
-                      elements=elements, filter=filter, z=z, maxisotopes=maxisotopes)
136
+                      elements=elements, filter=filter, z=z, maxisotopes=maxisotopes,
137
+                      minElements=minElements, maxElements=maxElements)
136 138
 }
137 139
 
138 140
 decomposeIsotopes <- function(masses, intensities, ppm=2.0, mzabs=0.0001,
139
-                              elements=NULL, filter=NULL, z=0, maxisotopes=10)
141
+                              elements=NULL, filter=NULL, z=0, maxisotopes=10,
142
+                              minElements="C0", maxElements="C999999")
140 143
 {
141 144
     # Use limited limited CHNOPS unless stated otherwise
142 145
     if (!is.list(elements) || length(elements)==0 ) {
... ...
@@ -174,6 +177,7 @@ decomposeIsotopes <- function(masses, intensities, ppm=2.0, mzabs=0.0001,
174 177
     molecules <- .Call("decomposeIsotopes",
175 178
                        masses, intensities, ppm, elements, element_order, z,
176 179
                        maxisotopes,
180
+                       minElements, maxElements,
177 181
                        PACKAGE="Rdisop")
178 182
 
179 183
     molecules
180 184
new file mode 100644
... ...
@@ -0,0 +1,15 @@
1
+test.empty <- function() {
2
+   checkEquals(length(decomposeMass(12)$formula), 1)
3
+}
4
+
5
+test.exact <- function() {
6
+   checkEquals(length(decomposeMass(12, minElements="C", maxElements="C")$formula), 1)
7
+}
8
+
9
+test.remove <- function() {
10
+  checkEquals(length(decomposeMass(12, minElements="C2", maxElements="C4")$formula), 0)
11
+}
12
+
13
+## Martin: Add more unit tests
14
+
15
+# decomposeMass(1.00785, minElements="C0", maxElements="C99999")
... ...
@@ -8,9 +8,10 @@
8 8
   obtained e.g.\ by FTICR or TOF mass spectrometers 
9 9
 }
10 10
 \usage{
11
-decomposeMass(mass, ppm=2.0, mzabs=0.0001, elements=NULL, filter=NULL, z=0, maxisotopes = 10)
11
+decomposeMass(mass, ppm=2.0, mzabs=0.0001, elements=NULL, filter=NULL,
12
+z=0, maxisotopes = 10, minElements="C0", maxElements="C999999")
12 13
 decomposeIsotopes(masses, intensities, ppm=2.0, mzabs=0.0001,
13
-elements=NULL, filter=NULL,  z=0, maxisotopes = 10)
14
+elements=NULL, filter=NULL,  z=0, maxisotopes = 10, minElements="C0", maxElements="C999999")
14 15
 isotopeScore(molecule, masses, intensities, elements = NULL, filter = NULL, z = 0)
15 16
 }
16 17
 \arguments{
... ...
@@ -24,6 +25,8 @@ isotopeScore(molecule, masses, intensities, elements = NULL, filter = NULL, z =
24 25
   \item{maxisotopes}{maximum number of isotopes shown in the resulting
25 26
     molecules}
26 27
   \item{elements}{list of allowed chemical elements, defaults to CHNOPS}
28
+  \item{minElements, maxElements}{Molecular formulas, which contain
29
+    lower and upper boundaries of allowed formula respectively}
27 30
   \item{filter}{NYI, will be a selection of DU, DBE and Nitrogen rules}
28 31
   \item{molecule}{a molecule as obtained from getMolecule() or
29 32
     decomposeMass / decomposeIsotopes}
... ...
@@ -79,7 +79,6 @@ float getDBE(const ComposedElement& molecule, int z) {
79 79
 
80 80
 // }}}
81 81
 
82
-
83 82
 char getParity(const ComposedElement& molecule, int charge=0) {
84 83
   // {{{ 
85 84
 
... ...
@@ -95,7 +94,7 @@ char getParity(const ComposedElement& molecule, int charge=0) {
95 94
 // }}}
96 95
 
97 96
 bool isValidMyNitrogenRule(const ComposedElement& molecule, int z) {
98
-  // {{{ 
97
+  // {{{
99 98
 
100 99
   bool massodd =  static_cast<int>(molecule.getNominalMass()) % 2 == 1 ? true : false;
101 100
   bool masseven = !massodd;
... ...
@@ -117,6 +116,30 @@ bool isValidMyNitrogenRule(const ComposedElement& molecule, int z) {
117 116
 }
118 117
 // }}}
119 118
 
119
+bool isWithinElementRange(const ComposedElement& molecule, const ComposedElement& minElements, const ComposedElement& maxElements) {
120
+// {{{
121
+
122
+  for (ComposedElement::container::const_iterator it = molecule.getElements().begin(); 
123
+       it != molecule.getElements().end(); ++it) {
124
+    
125
+    int mincount = static_cast<int>(minElements.getElementAbundance((it->first).getName()));
126
+    int maxcount = static_cast<int>(maxElements.getElementAbundance((it->first).getName()));
127
+
128
+    int count = static_cast<int>(molecule.getElementAbundance((it->first).getName()));
129
+      
130
+      if (count < mincount) {
131
+	return false;
132
+      }
133
+
134
+      // TODO: Fails e.g. for "C2N0" 
135
+      if (maxcount>0 && count > maxcount) {
136
+	return false;
137
+      }
138
+  }
139
+
140
+  return true;
141
+}
142
+// }}}
120 143
 
121 144
 //
122 145
 // Decomposition of Mass / Isotope Pattern
... ...
@@ -124,7 +147,8 @@ bool isValidMyNitrogenRule(const ComposedElement& molecule, int z) {
124 147
 
125 148
 RcppExport SEXP decomposeIsotopes(SEXP v_masses, SEXP v_abundances, SEXP s_error, 
126 149
 				  SEXP l_alphabet, SEXP v_element_order, 
127
-				  SEXP z, SEXP i_maxisotopes) {
150
+				  SEXP z, SEXP i_maxisotopes,
151
+				  SEXP s_minElements, SEXP s_maxElements) {
128 152
 // {{{ 
129 153
 
130 154
     typedef DistributionProbabilityScorer scorer_type;
... ...
@@ -229,15 +253,29 @@ RcppExport SEXP decomposeIsotopes(SEXP v_masses, SEXP v_abundances, SEXP s_error
229 253
 	// - chemical filter is applied
230 254
 	// - isotopic pattern is calculated
231 255
 	// - isotopic pattern is matched against input spectrum
256
+
257
+	// Initialize minimum/maximum element count "molecules"
258
+	ComposedElement minElements(CHAR(Rf_asChar(s_minElements)), alphabet);
259
+	ComposedElement maxElements(CHAR(Rf_asChar(s_maxElements)), alphabet);
260
+
232 261
 	for (decompositions_t::iterator decomps_it = decompositions.begin(); 
233 262
 		decomps_it != decompositions.end(); ++decomps_it) {
234 263
 
235 264
 		// creates a candidate molecule out of elemental composition and a set of elements
236 265
 		ComposedElement candidate_molecule(*decomps_it, alphabet);
266
+
267
+		// Check minimum/maximum element counts
268
+		if (!isWithinElementRange(candidate_molecule, minElements, maxElements)) {
269
+			continue;
270
+		} 
271
+
272
+
237 273
 		// checks on chemical filter
238 274
 // 		if (!isValidMyNitrogenRule(candidate_molecule, z)) {
239 275
 // 			continue;
240 276
 // 		} 
277
+
278
+
241 279
 		// updates molecules isotope distribution (since its not calculated upon creation: 
242 280
 		// it would be time consuming before applying chemical filter)
243 281
 		candidate_molecule.updateIsotopeDistribution();
... ...
@@ -383,7 +421,7 @@ RcppExport SEXP calculateScore(SEXP v_predictMasses, SEXP v_predictAbundances, S
383 421
 
384 422
 RcppExport SEXP getMolecule(SEXP s_formula, SEXP l_alphabet, 
385 423
 			    SEXP v_element_order, SEXP z, SEXP i_maxisotopes) {
386
-  // {{{ 
424
+// {{{ 
387 425
 
388 426
   SEXP  rl=R_NilValue; // Use this when there is nothing to be returned.
389 427
 
... ...
@@ -819,7 +857,7 @@ extern "C" {
819 857
       {"getMolecule", (void* (*)())&getMolecule, 4},
820 858
       {"addMolecules", (void* (*)())&addMolecules, 4},
821 859
       {"subMolecules", (void* (*)())&subMolecules, 4},
822
-      {"decomposeIsotopes", (void* (*)())&decomposeIsotopes, 7},
860
+      {"decomposeIsotopes", (void* (*)())&decomposeIsotopes, 9},
823 861
       {"calculateScore", (void* (*)())&calculateScore, 7},
824 862
       {NULL, NULL, 0}
825 863
     };