diff MassbankEntry.R @ 2:20d69a062da3 draft

planemo upload for repository https://github.com/workflow4metabolomics/lcmsmatching.git commit d4048accde6bdfd5b3e14f5394902d38991854f8
author prog
date Thu, 02 Mar 2017 08:55:00 -0500
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MassbankEntry.R	Thu Mar 02 08:55:00 2017 -0500
@@ -0,0 +1,129 @@
+###########################
+# MASSBANK SPECTRUM CLASS #
+###########################
+
+MassbankEntry <- methods::setRefClass("MassbankEntry", contains = "BiodbEntry")
+
+###########
+# FACTORY #
+###########
+
+createMassbankEntryFromTxt <- function(contents, drop = TRUE) {
+
+	entries <- list()
+
+	# Define fields regex
+	regex <- character()
+	regex[[BIODB.ACCESSION]] <- "^ACCESSION: (.+)$"
+	regex[[BIODB.MSDEV]] <- "^AC\\$INSTRUMENT: (.+)$"
+	regex[[BIODB.MSDEVTYPE]] <- "^AC\\$INSTRUMENT_TYPE: (.+)$"
+	regex[[BIODB.MSTYPE]] <- "^AC\\$MASS_SPECTROMETRY: MS_TYPE (.+)$"
+	regex[[BIODB.MSPRECMZ]] <- "^MS\\$FOCUSED_ION: PRECURSOR_M/Z (.+)$"
+	regex[[BIODB.NB.PEAKS]] <- "^PK\\$NUM_PEAK: ([0-9]+)$"
+	regex[[BIODB.MSPRECANNOT]] <- "^MS\\$FOCUSED_ION: PRECURSOR_TYPE (.+)$"
+	regex[[BIODB.CHEBI.ID]] <- "^CH\\$LINK: CHEBI\\s+(.+)$"
+	regex[[BIODB.KEGG.ID]] <- "^CH\\$LINK: KEGG\\s+(.+)$"
+	regex[[BIODB.INCHI]] <- "^CH\\$IUPAC:\\s+(.+)$"
+	regex[[BIODB.INCHIKEY]] <- "^CH\\$LINK: INCHIKEY\\s+(.+)$"
+	regex[[BIODB.CHEMSPIDER.ID]] <- "^CH\\$LINK: CHEMSPIDER\\s+(.+)$"
+	regex[[BIODB.CAS.ID]] <- "^CH\\$LINK: CAS\\s+(.+)$"
+	regex[[BIODB.FORMULA]] <- "^CH\\$FORMULA:\\s+(.+)$"
+	regex[[BIODB.SMILES]] <- "^CH\\$SMILES:\\s+(.+)$"
+	regex[[BIODB.MASS]] <- "^CH\\$EXACT_MASS:\\s+(.+)$"
+	regex[[BIODB.PUBCHEMCOMP.ID]] <- "^CH\\$LINK: PUBCHEM\\s+.*CID:([0-9]+)"
+	regex[[BIODB.PUBCHEMSUB.ID]] <- "^CH\\$LINK: PUBCHEM\\s+.*SID:([0-9]+)"
+
+	for (text in contents) {
+
+		# Create instance
+		entry <- MassbankEntry$new()
+
+		if ( ! is.null(text) && ! is.na(text)) {
+
+			# Read text
+			lines <- strsplit(text, "\n")
+			for (s in lines[[1]]) {
+
+				# Test generic regex
+				parsed <- FALSE
+				for (field in names(regex)) {
+					g <- stringr::str_match(s, regex[[field]])
+					if ( ! is.na(g[1,1])) {
+						entry$setField(field, g[1,2])
+						parsed <- TRUE
+						break
+					}
+				}
+				if (parsed)
+					next
+
+				# Name
+				if (is.na(entry$getField(BIODB.NAME))) {
+					g <- stringr::str_match(s, "^CH\\$NAME:\\s+(.+)$")
+					if ( ! is.na(g[1,1]))
+						entry$setField(BIODB.NAME, g[1,2])
+				}
+		
+				# PubChem
+				g <- stringr::str_match(s, "^CH\\$LINK: PUBCHEM\\s+([0-9]+)$")
+				if ( ! is.na(g[1,1]))
+					entry$setField(BIODB.PUBCHEMSUB.ID, g[1,2])
+
+				# MS MODE
+				g <- stringr::str_match(s, "^AC\\$MASS_SPECTROMETRY: ION_MODE (.+)$")
+				if ( ! is.na(g[1,1])) {
+					entry$setField(BIODB.MSMODE, if (g[1,2] == 'POSITIVE') BIODB.MSMODE.POS else BIODB.MSMODE.NEG)
+					next
+				}
+
+				# PEAKS
+				if (.parse.peak.line(entry, s))
+					next
+			}
+		}
+
+		entries <- c(entries, entry)
+	}
+
+	# Replace elements with no accession id by NULL
+	entries <- lapply(entries, function(x) if (is.na(x$getField(BIODB.ACCESSION))) NULL else x)
+
+	# If the input was a single element, then output a single object
+	if (drop && length(contents) == 1)
+		entries <- entries[[1]]
+
+	return(entries)
+}
+
+###################
+# PARSE PEAK LINE #
+###################
+
+.parse.peak.line <- function(entry, line) {
+
+	peaks <- BIODB.PEAK.DF.EXAMPLE
+	
+	# Annotation
+	g <- stringr::str_match(line, "^\\s+([0-9][0-9.]*) ([A-Z0-9+-]+) ([0-9]+) ([0-9][0-9.]*) ([0-9][0-9.]*)$")
+	if ( ! is.na(g[1,1]))
+		peaks[1, c(BIODB.PEAK.MZ, BIODB.PEAK.FORMULA, BIODB.PEAK.FORMULA.COUNT, BIODB.PEAK.MASS, BIODB.PEAK.ERROR.PPM)] <- list(as.double(g[1,2]), g[1,3], as.integer(g[1,4]), as.double(g[1,5]), as.double(g[1,6]))
+
+	# Peak
+	g <- stringr::str_match(line, "^\\s+([0-9][0-9.]*) ([0-9][0-9.]*) ([0-9]+)$")
+	if ( ! is.na(g[1,1]))
+		peaks[1, c(BIODB.PEAK.MZ, BIODB.PEAK.INTENSITY, BIODB.PEAK.RELATIVE.INTENSITY)] <- list(as.double(g[1,2]), as.double(g[1,3]), as.integer(g[1,4]))
+
+	if (nrow(peaks) > 0) {
+
+		# Get curent peaks and merge with new peaks
+		current.peaks <- entry$getField(BIODB.PEAKS)
+		if ( ! is.null(current.peaks))
+			peaks <- rbind(current.peaks, peaks)
+
+		entry$setField(BIODB.PEAKS, peaks)
+
+		return(TRUE)
+	}
+
+	return(FALSE)
+}