Merge pull request #3 from targeted-lipidomics/dev

Release v0.1.0

DESCRIPTION

@@ -11,13 +11,16 @@ License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a
     license
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.0
 Collate: 
     'setClasses.R'
     'createMSIDatamatrix.R'
     'create_cal_curve.R'
+    'extdata.R'
     'int2conc.R'
     'int2response.R'
+    'read_mrm.R'
     'remove_blank_mzs.R'
+    'setCommonAxis.R'
     'summarise_cal_levels.R'
     'zero2na.R'

NAMESPACE

@@ -2,6 +2,7 @@
 
 export(calibrationInfo)
 export(quant_MSImagingExperiment)
+export(read_mrm)
 export(tissueInfo)
 exportClasses(calibrationInfo)
 exportClasses(quant_MSImagingExperiment)
@@ -11,6 +12,7 @@ exportMethods(create_cal_curve)
 exportMethods(int2conc)
 exportMethods(int2response)
 exportMethods(remove_blank_mzs)
+exportMethods(setCommonAxis)
 exportMethods(summarise_cal_levels)
 exportMethods(zero2na)
 import(Cardinal)

R/createMSIDatamatrix.R

@@ -8,43 +8,67 @@ setGeneric("createMSIDatamatrix", function(MSIobject, ...) standardGeneric("crea
 #' @import dplyr
 #' @include setClasses.R
 #'
-#' @param MSIobject MSI object from Cardinal
+#' @param MSIobject MSI object from Cardinal, pData to include..... sample_ID..
 #' @param inputNA Whether to convert 0's in matrix to NA (default = T)
+#' @param roi_header Header in pData pertaining to ROIs to average. Set to NA to skip generating average df
 #' @return MSIobject with slots updated for i) matrix of average ng/pixel of m/z (rows = m/z and cols = cal level) in tissue ROIs ii) sample/ROI metadata
 #'
 #' @export
 setMethod("createMSIDatamatrix", "quant_MSImagingExperiment",
-          function(MSIobject, inputNA = T){
+          function(MSIobject, inputNA = T, roi_header = NA){
 
+            # Subset pixels in ROIs only (based on roi_header)
+            if(is.na(roi_header)){
+              pData(MSIobject)$ROI = 1:nrow(pData(MSIobject))
+            } else{
+              pData(MSIobject)$ROI = pData(MSIobject)[[roi_header]]
+              MSIobject = MSIobject[, - which(is.na(pData(MSIobject)$ROI))]
+            }
 
-            MSIobject = MSIobject[, - which(is.na(pData(MSIobject)$ROI))]
-
+            # Update pixel data
             pixel_df = data.frame(pData(MSIobject)) %>%
               subset(!is.na(ROI)) %>%
               tibble::rownames_to_column("pixel_ind")
 
-            df = data.frame(matrix(ncol = length(unique(pData(MSIobject)$sample_ID)),
-                                   nrow = nrow(fData(MSIobject))))
 
-            colnames(df) = unique(pData(MSIobject)$sample_ID)
-            rownames(df) = fData(MSIobject)@mz
+            # All pixel df
+            all_pixel_df = data.frame(spectra(MSIobject)[,])
+            colnames(all_pixel_df ) = sprintf("pixel_%s", pixel_df$pixel_ind)
+            rownames(all_pixel_df ) = fData(MSIobject)$name
 
-            for(roi in unique(pixel_df$sample_ID)){
+            if(inputNA){
+              all_pixel_df <- replace(all_pixel_df, all_pixel_df==0, NA)
+            }
 
-              pixel_inds = as.numeric(pixel_df$pixel_ind[which(pixel_df$sample_ID == roi)])
 
-              temp_df = data.frame(spectra(MSIobject)[, pixel_inds])
-              rowMeans(temp_df, na.rm=T)
+            # Create empty average df
+            if(!is.na(roi_header)){
+              ave_df = data.frame(matrix(ncol = length(unique(pData(MSIobject)$sample_ID)),
+                                         nrow = nrow(fData(MSIobject))))
 
-              df[[roi]] = rowMeans(temp_df, na.rm=T)
+              colnames(ave_df) = unique(pData(MSIobject)$sample_ID)
+              rownames(ave_df) = fData(MSIobject)$name
 
-            }
 
-            if(inputNA){
-              df <- replace(df, df==0, NA)
+              for(roi in unique(pixel_df$sample_ID)){
+
+                pixel_inds = as.numeric(pixel_df$pixel_ind[which(pixel_df$sample_ID == roi)])
+
+                temp_df = data.frame(spectra(MSIobject)[, pixel_inds])
+                rowMeans(temp_df, na.rm=T)
+
+                ave_df[[roi]] = rowMeans(temp_df, na.rm=T)
+
+              }
+
+              if(inputNA){
+                ave_df <- replace(ave_df, ave_df==0, NA)
+              }
+
+              MSIobject@tissueInfo@roi_average_matrix = ave_df
             }
 
-            MSIobject@tissueInfo@conc_matrix = df
+            MSIobject@tissueInfo@all_pixel_matrix = all_pixel_df
             MSIobject@tissueInfo@sample_metadata = pixel_df
 
             return(MSIobject)

R/create_cal_curve.R

@@ -12,7 +12,7 @@ setGeneric("create_cal_curve", function(MSIobject, ...) standardGeneric("create_
 #' @include setClasses.R
 #'
 #' @param response_matrix matrix of average ng/pixel of m/z (rows = m/z and cols = cal level)
-#' @param cal_type string of approach to generate claibration curve - 'std_addition' is default
+#' @param cal_type string of approach to generate claibration curve - 'std_addition' if standards are on tissue and 'cal' if direct onto glass slide.
 #' @return MSIobject with slots updated for i) cal_list - List of linear models for each m/z (response v concentration, where concentration is ng/pixel) and ii) r2 values for each calibration iii) calibration metadata
 #'
 #' @export
@@ -23,8 +23,8 @@ setMethod("create_cal_curve", "quant_MSImagingExperiment",
 
             cal_metadata = MSIobject@calibrationInfo@cal_metadata %>%
               mutate(pixel_count = sapply(sample,
-                                          function(sample_name) pixel_count[[sample_name]]),
-                     ng_per_pixel = amount_ng / pixel_count)
+                                          function(sample_name) pixel_count[[sample_name]])) %>%
+              mutate(ng_per_pixel = amount_ng / pixel_count)
 
             background_sample = cal_metadata$sample[which(cal_metadata$ng_per_pixel == 0)]
 

R/extdata.R

@@ -0,0 +1,82 @@
+#' Dataset tissue_MRM_data
+#'
+#' This is a dataset containing Waters MRM data of SphingoMyelin 16:00 in Guinea pig lung tissue.
+#'
+#' @name tissue_MRM_data
+#'
+#' @section tissue_MRM_data.raw/imaging/Analyte 1.txt:
+#'
+#' This data is used in read_mrm().
+NULL
+
+
+#' Dataset cal_MRM_data
+#'
+#' This is a dataset containing Waters MRM data of SphingoMyelin 16:00 in Guinea pig lung tissue.
+#'
+#' @name cal_MRM_data
+#'
+#' @section cal_MRM_data.raw/imaging/Analyte 1.txt:
+#'
+#' This data is used in read_mrm().
+NULL
+
+
+#' Dataset ion_library
+#'
+#' This is a table containing a library of MRM transitions included in targeted MSI experiments.
+#'
+#' @name ion_library
+#'
+#' @section ion_library.txt:
+#'
+#' This data is used in read_mrm().
+NULL
+
+
+#' Dataset cal_rois
+#'
+#' This is a table containing a logical data about caliobration levels pixels are associated with.
+#'
+#' @name cal_rois
+#'
+#' @section cal_rois.csv:
+#'
+#' This data is used to identify calibration spots.
+NULL
+
+
+#' Dataset calibration_metadata
+#'
+#' This is a table containing information about concentrations in the calibration spots.
+#'
+#' @name calibration_metadata
+#'
+#' @section calibration_metadata.csv:
+#'
+#' This data is used to create calibration curves in create_cal_curve().
+NULL
+
+
+#' Dataset tissue_pixels
+#'
+#' This is a table containing a logical data about whether pixel is part of a calibration level.
+#'
+#' @name tissue_pixels
+#'
+#' @section tissue_pixels.csv:
+#'
+#' This data is used to identify tissue pixels.
+NULL
+
+
+#' Dataset tissue_rois
+#'
+#' This is a table containing a logical data about which tissue types pixels are associated with.
+#'
+#' @name tissue_rois
+#'
+#' @section tissue_rois.csv:
+#'
+#' This data is used to identify and label tissue types.
+NULL

R/read_mrm.R

@@ -1,11 +1,27 @@
-read_mrm = function(name,
-                    folder){
+require(Cardinal)
+
+#' Function to create data matrix from MSI object
+#' @import Cardinal
+#' @import dplyr
+#' @include setClasses.R
+#'
+#' @param name Name of Waters raw imaging file (including the imaging/Analyte 1.txt file after processing)
+#' @param folder Location of folder with raw imaging data in
+#' @param lib_ion_path Full path to library of transitions for targeted MSI experiments. Must include headers: 'transition_id',	'precursor_mz',	'product_mz',	'collision_eV',	'cone_V',	'Polarity',	'Type'. Where 'Type' is "Analyte" for analytes in tissue and "IS" for internal standards.
+#' @param polarity String indicating whether data was collected in 'Positive' or 'Negative' polarity.
+
+#' @return MSIobject with slots updated for i) matrix of average ng/pixel of m/z (rows = m/z and cols = cal level) in tissue ROIs ii) sample/ROI metadata
+#'
+#' @export read_mrm
+read_mrm = function(name, folder, lib_ion_path, polarity){
+
+  ion_lib = read.table(file = lib_ion_path, sep = "\t", header = T)
 
   imaging_folder = sprintf("%s/%s.raw/imaging", folder, name)
 
   analyte_fn = list.files(imaging_folder, full.names = T)
 
-  analyte_df = read.table(analyte_fn, skip=2, check.names = F, fill = TRUE, sep="\t", header=F)
+  analyte_df = read.table(analyte_fn, fill = TRUE, sep="\t", header=F, blank.lines.skip = T)[-1,]
 
   transitions = t(analyte_df[1:3,]) %>%
     `colnames<-`(c("transition_id", "precursor_mz", "product_mz")) %>%
@@ -36,11 +52,21 @@ read_mrm = function(name,
   run <- factor(rep(name, nrow(coord)))
   pdata <- PositionDataFrame(run=run, coord=coord)
 
+  # Gather info of MRM transitions from the ion library
+  ion_lib = ion_lib %>%
+    subset(Polarity == polarity) %>%
+    dplyr::right_join(y=transitions,  by = c('precursor_mz', 'product_mz'), suffix = c("_name", "_int")) %>%
+    mutate(transition_id_name = ifelse(is.na(transition_id_name), transition_id_int, transition_id_name),
+           Polarity = ifelse(is.na(Polarity), polarity, Polarity),
+           Type = ifelse(is.na(Type), "Unknown", Type)) %>%
+    arrange(transition_id_int)
+
   # feature metadata
-  fdata <- MassDataFrame(mz=transitions$transition_id,
-                         analyte = "analyte",
-                         precursor_mz = transitions$precursor_mz,
-                         product_mz = transitions$product_mz)
+  fdata <- MassDataFrame(mz=ion_lib$transition_id_int,
+                         analyte = ion_lib$Type,
+                         precursor_mz = ion_lib$precursor_mz,
+                         product_mz = transitions$product_mz,
+                         name = ion_lib$transition_id_name)
 
   # intensity data
   idata = t(analyte_df[, grep("transition", colnames(analyte_df))])

R/setClasses.R

@@ -24,7 +24,8 @@ calibrationInfo = setClass("calibrationInfo",
 #' @export
 tissueInfo = setClass("tissueInfo",
                            slots = c(
-                             conc_matrix = "data.frame",
+                             roi_average_matrix = "data.frame",
+                             all_pixel_matrix = "data.frame",
                              sample_metadata = "data.frame",
                              feature_metadata = "data.frame"
                            )

R/setCommonAxis.R

@@ -0,0 +1,60 @@
+library(Cardinal)
+
+setGeneric("setCommonAxis", function(MSIobjects, ...) standardGeneric("setCommonAxis"))
+
+#' Function to update intensity with concentration values
+#' @import Cardinal
+#' @include setClasses.R
+#'
+#' @param MSIobjects List of MSIobjects
+#' @param ref_object MSIobject with all transitions of interets included
+#' @return List of MSIobjects with common fData and matching spectra
+#'
+#' @export
+setMethod("setCommonAxis", "list",
+          function(MSIobjects, ref_object){
+
+            features = lapply( 1:length(MSIobjects), FUN=function(x){
+              data.frame(fData(MSIobjects[[x]])) })
+
+            ref_features = data.frame(fData(ref_object))
+
+            # Common axis
+            for(i in 1:length(MSIobjects)){
+
+              if(all(dim(ref_features) == dim(features[[i]]))){
+                if(all(ref_features == features[[i]])) next
+              }
+
+              feat = features[[i]]
+
+              # Correct fData
+              feat = merge(x=ref_features, y=feat, by = c("precursor_mz", "product_mz"),
+                           all.x = T, all.y = F, suffixes = c("","_old")) %>%
+                mutate(name = name_old, analyte = analyte_old) %>%
+                arrange(mz)
+
+              # Set empty iData
+              idata = matrix(nrow=nrow(feat), ncol=ncol(MSIobjects[[i]]))
+
+              # Update iData to spectral info corresponding to fData channels
+              for(mz in (feat %>% subset(!is.na(analyte)))$mz){
+
+                old_mz = feat$mz_old[which(feat$mz == mz)]
+
+                idata[mz, ] = spectra(MSIobjects[[i]][old_mz,])
+
+              }
+
+              MSIobjects[[i]] <- MSImagingExperiment(imageData= idata,
+                                         featureData= MassDataFrame(mz=as.numeric(feat$mz),
+                                                                   analyte = feat$analyte,
+                                                                   precursor_mz = as.numeric(feat$precursor_mz),
+                                                                   product_mz = as.numeric(feat$product_mz),
+                                                                   name = feat$name),
+                                         pixelData=pData(MSIobjects[[i]]))
+            }
+
+            return(MSIobjects)
+
+          })

README.rst

@@ -2,7 +2,8 @@
 quantMSImageR
 ==============================================
 
-Software tools (largely extended from the R package Cardinal) for processing and quantification of targeted (multiple reaction monitoring) mass spectrometry imaging (MSI) data.
+=======
+Software tools for processing and quantifying targeted multiple reaction monitoring (MRM) mass spectrometry imaging (MSI) data.
 
 ------------
 Install