Spectra

In R, we have two ways to visualize your spectra, one is quick and efficient and the other is interactive. Here is an example of quick and efficient plotting.

plot(scratch_OpenSpecy) # quick and efficient

This is an example of an interactive plot. You can plot two different datasets simultaneously to compare.

# This will min-max normalize your data even if it isn't already but are not
# changing your underlying data
plotly_spec(scratch_OpenSpecy, json_example)

Maps

Spectral maps can also be visualized as overlaid spectra but in addition the spatial information can be plotted as a heatmap. It is important to note that when multiple spectra are uploaded in batch they are prescribed x and y coordinates, this can be helpful for visualizing summary statistics and navigating vast amounts of data but shouldn’t be confused with data which actually has spatial coordinates.

In R the user can control what values form the colors of the heatmap by specifying z, it is handy to put the information you want in the metadata for this reason. This example just shows the x values of the spectra.

heatmap_spec(spectral_map,
             z = spectral_map$metadata$x)

An interactive plot of the heatmap and spectra overlayed can be generated with the interactive_plot() function. A user can hover over the heatmap to identify the next row id they are interested in and update the value of select to see that spectrum.

interactive_plot(spectral_map, select = 100, z = spectral_map$metadata$x)

Threshold Signal-Noise

Considering whether you have enough signal to analyze spectra is important. Classical spectroscopy would recommend your highest peak to be at least 10 times the baseline of your processed spectra before you begin analysis. If your spectra is below that threshold even after processing, you may want to consider recollecting it. In practice, we are rarely able to collect spectra of that good quality and more often use 5. The “run_sig_over_noise” metric searches your spectra for high and low regions and conducts division on them to derive the signal to noise ratio and step specifies the number of intensities to search. This is a good way to automatically calculate the signal to noise ratio. In the example below you can see that our signal to noise ratio is increased by the processing, the goal of processing is generally to maximize this. “sig_times_noise” metric multiplies the mean signal by the standard deviation of the signal and “tot_sig” metric sums the intensities. The latter can be really useful for thresholding spectral maps to identify particles which we will discuss later. If you know where your signal region and noise regions are you can specify them with sig_min, sig_max, noise_min, and noise_max.

sig_noise(processed, metric = "run_sig_over_noise") > 
  sig_noise(raman_hdpe, metric = "run_sig_over_noise")

If analyzing spectra in batch, we recommend looking at the heatmap and optimizing the percent of spectra that are above your signal to noise threshold to determine the correct settings instead of looking through spectra individually. Setting the min_sn will threshold the heatmap image to only color spectra which have a sn value over the threshold.

spectral_map_p <- spectral_map |>
  process_spec(flatten_range = T)

spectral_map_p$metadata$sig_noise <- sig_noise(spectral_map_p)

heatmap_spec(spectral_map_p, sn = spectral_map_p$metadata$sig_noise, min_sn = 5)

Intensity Adjustment

Open Specy assumes that intensity units are in absorbance units but Open Specy can adjust reflectance or transmittance spectra to absorbance units. The transmittance adjustment uses the \(\log_{10} 1/T\) calculation which does not correct for system or particle characteristics. The reflectance adjustment use the Kubelka-Munk equation \(\frac{(1-R)^2}{2R}\).

This is the respective R code for a scenario where the spectra doesn’t need intensity adjustment:

trans_raman_hdpe <- raman_hdpe
trans_raman_hdpe$spectra <- 2 - trans_raman_hdpe$spectra^2
    
rev_trans_raman_hdpe <- trans_raman_hdpe |> 
  adj_intens(type = "transmittance")

plotly_spec(trans_raman_hdpe, rev_trans_raman_hdpe)

Conforming

Conforming spectra is essential before comparing to a reference library and can be useful for summarizing data when you don’t need it to be highly resolved spectrally. We set the default spectral resolution to 8 because this tends to be pretty good for a lot of applications and is in between 4 and 8 which are commonly used wavenumber resolutions. There are several ways that this function can be specified.

conform_spec(raman_hdpe, res = 8) |> # convert res to 8 wavenumbers.
  summary()
#> $wavenumber
#>  Length Min. Max.     Res.
#>     362  304 3192 7.977901
#> 
#> $spectra
#>  Number Min. Intensity Max. Intensity
#>       1       42.85152       602.4639
#> 
#> $metadata
#>   Min. Max.
#> x    1    1
#> y    1    1
#> [1] "x"                 "y"                 "user_name"        
#> [4] "spectrum_type"     "spectrum_identity" "organization"     
#> [7] "license"           "session_id"        "file_id"

# Force one spectrum to have the exact same wavenumbers as another
conform_spec(asp_example, range = ps_example$wavenumber, res = NULL) |> 
  summary()
#> $wavenumber
#>  Length     Min.     Max.     Res.
#>    1736 651.8728 3998.025 1.927507
#> 
#> $spectra
#>  Number Min. Intensity Max. Intensity
#>       1   0.0009781419      0.5156938
#> 
#> $metadata
#>   Min. Max.
#> x    1    1
#> y    1    1
#> [1] "x"          "y"          "file_name"  "license"    "col_id"    
#> [6] "session_id" "file_id"

# Specify the wavenumber resolution and use a rolling join instead of linear
# approximation (faster for large datasets). 
conform_spec(spectral_map, range = ps_example$wavenumber, res = 10,
             type = "roll") |> 
  summary()
#> $wavenumber
#>  Length Min. Max.     Res.
#>     329  720 4000 9.969605
#> 
#> $spectra
#>  Number Min. Intensity Max. Intensity
#>     208      -1.317072       1.168225
#> 
#> $metadata
#>   Min. Max.
#> x    0   15
#> y    0   12
#>  [1] "y"             "x"             "file_name"     "license"      
#>  [5] "description"   "samples"       "lines"         "bands"        
#>  [9] "header offset" "data type"     "interleave"    "z plot titles"
#> [13] "pixel size"    "col_id"        "session_id"    "file_id"

Smoothing

The Savitzky-Golay filter is used for smoothing. Higher polynomial numbers lead to more wiggly fits and thus less smoothing, lower numbers lead to more smooth fits. The SG filter is fit to a moving window of 11 data points by default where the center point in the window is replaced with the polynomial estimate. Larger windows will produce smoother fits. The derivative order is set to 1 by default which transforms the spectra to their first derivative. A zero order derivative will have no derivative transformation. When smoothing is done well, peak shapes and relative heights should not change. The absolute value is primarily useful for first derivative spectra where the absolute value results in an absorbance-like spectrum which is why we set it as the default. You’ll notice a new function we are using c_spec() which is used to combine spectral objects into one OpenSpecy object.

Examples of smoothing:

none <- make_rel(raman_hdpe)
p1 <- smooth_intens(raman_hdpe, polynomial = 1, derivative = 0, abs = F)
p4 <- smooth_intens(raman_hdpe, polynomial = 4, derivative = 0, abs = F)

c_spec(list(none, p1, p4)) |> 
  plot()

Sample raman_hdpe spectrum with different smoothing polynomials.

Derivative transformation can happen with the same function.

none <- make_rel(raman_hdpe)
d1 <- smooth_intens(raman_hdpe, derivative = 1, abs = T)
d2 <- smooth_intens(raman_hdpe, derivative = 2, abs = T)

c_spec(list(none, d1, d2)) |> 
  plot()

Sample raman_hdpe spectrum with different derivatives.

Baseline Correction

The goal of baseline correction is to get all non-peak regions of the spectra to zero absorbance. The higher the polynomial order, the more wiggly the fit to the baseline. If the baseline is not very wiggly, a more wiggly fit could remove peaks which is not desired. The baseline correction algorithm used in Open Specy is called “iModPolyfit” (Zhao et al. 2007). This algorithm iteratively fits polynomial equations of the specified order to the whole spectrum. During the first fit iteration, peak regions will often be above the baseline fit. The data in the peak region is removed from the fit to make sure that the baseline is less likely to fit to the peaks. The iterative fitting terminates once the difference between the new and previous fit is small. An example of a good baseline fit below. Manual baseline correction can also be specified by providing a baseline OpenSpecy object.

alternative_baseline <- smooth_intens(raman_hdpe, polynomial = 1, window = 51,
                                      derivative = 0, abs = F, make_rel = F) |>
  flatten_range(min = 2700, max = 3200, make_rel = F)

none <- make_rel(raman_hdpe)
d2 <- subtr_baseline(raman_hdpe, type = "manual",
                     baseline = alternative_baseline)
d8 <- subtr_baseline(raman_hdpe, degree = 8)

c_spec(list(none, d2, d8)) |> 
  plot()

Sample raman_hdpe spectrum with different degrees of background subtraction (Cowger et al., 2020).

Range Selection

Sometimes an instrument operates with high noise at the ends of the spectrum and, a baseline fit produces distortions, or there are regions of interest for analysis. Range selection accomplishes those goals. You should look into the signal to noise ratio of your specific instrument by wavelength to determine what wavelength ranges to use. Distortions due to baseline fit can be assessed from looking at the process spectra. Additionally, you can restrict the range to examine a single peak or a subset of peaks of interests. Multiple ranges can be specified simultaneously.

none <- make_rel(raman_hdpe)
r1 <- restrict_range(raman_hdpe, min = 1000, max = 2000)
r2 <- restrict_range(raman_hdpe, min = c(1000, 1800), max = c(1200, 2000))

compare_ranges <- c_spec(list(none, r1, r2), range = "common")
# Common argument crops the ranges to the most common range between the spectra
# when joining. 

plot(compare_ranges)

Sample raman_hdpe spectrum with different degrees of range restriction.

Flattening Ranges

Sometimes there are peaks that really shouldn’t be in your spectra and can distort your interpretation of the spectra but you don’t necessarily want to remove the regions from the analysis because you believe those regions should be flat instead of having a peak. One way to deal with this is to replace the peak values with the mean of the values around the peak. This is the purpose of the flatten_range function. By default it is set to flatten the CO2 region for FTIR spectra because that region often needs to be flattened when atmospheric artifacts occur in spectra. Like restrict_range, the R function can accept multiple ranges.

single <- filter_spec(spectral_map, 120) # Function to filter spectra by index
                                         # number or name or a logical vector. 
none <- make_rel(single)
f1 <- flatten_range(single) #default flattening the CO2 region. 
f2 <- flatten_range(single, min = c(1000, 2500), max = c(1200, 3000))

compare_flats <- c_spec(list(none, f1, f2))

plot(compare_flats)

Sample raman_hdpe spectrum with different degrees of background subtraction (Cowger et al., 2020).

Min-Max Normalization

Often we regard spectral intensities as arbitrary and min-max normalization allows us to view spectra on the same scale without drastically distorting their shapes or relative peak intensities. In the package, most of the processing functions will min-max transform your spectra by default if you do not specify otherwise. This plot shows two plots that are nearly identical except for the min-max transformed spectrum has a y axis that ranges from 0-1.

relative <- make_rel(raman_hdpe)

Reading Libraries

Reference libraries are spectra with known identities. The Open Specy library now has over 10,000 spectra in it an is getting so large that we cannot fit it within the R package size limit of 5 MB. We host the reference libraries on OSF and have a function to pull the libraries down automatically. Running get_lib by itself will download all libraries to your package director or you can specify which libraries you want and where you want them.

get_lib(type = "derivative")

After download the libraries you can load them into your active environment. You should load them one at a time. Libraries are also OpenSpecy objects so you can use any Open Specy object as a library.

lib <- load_lib(type = "derivative")

Matches

Before attempting to use a reference library to identify spectra it is really important to understand what format the reference library is in. All the OpenSpecy reference libraries are in Absorbance units. derivative has been absolute first derivative transformed, nobaseline has been baseline corrected, raw is the rawest form of the reference spectra (not recommended except for advanced uses). The previously mentioned libraries all have Raman and FTIR spectra in them. mediod is the mediod compressed library version of the absolute first derivative for FTIR, model is an exception because it is a multinomial regression approach for FTIR to identification of absolute derivative spectra. In this example we use the data("test_lib") which is a subsampled version of the absolute derivative library and data("raman_hdpe") which is an unprocessed Raman spectrum in absorbance units of HDPE plastic. With single spectra it is easy to look at the spectra but when doing in batch, also refer to the Thresholding Signal-Noise section for guidance on making sure your batch spectra are processed to quality specs.

data("test_lib")
data("raman_hdpe")

processed <- process_spec(x = raman_hdpe, 
                          conform_spec = F, #We will conform during matching. 
                          smooth_intens = T #Conducts the default derivative transformation. 
                          )

# Check to make sure that the signal to noise ratio of the processed spectra is
# greater than 10.
print(sig_noise(processed) > 10)
plotly_spec(raman_hdpe, processed)

After your spectra is processed similarly to the library specifications, you can identify the spectra using match_spec(). All of the libraries have wavenumbers at 8 cm^-1 resolution. Whichever library you choose, you need to get your spectra into a similar enough format to use for comparison. That means conforming the wavenumbers to the same values using conform = T. Alternatively, you could have done the conformation during the processing. The add_library_metadata and add_object_metadata options specify the column name in the metadata that you want to add metadata from and top_n specifies how many matches you want. In this example we just identified a single spectrum with the library but you can also send an OpenSpecy object with multiple spectra. The output matches is a data.table with at least 3 columns, object_id tells you the column names of the spectra in x, library_id tells you the column names from the library that it matched to. match_val is the value of the Pearson correlation coefficient (default) or other correlation if specified in ... or if using the model identification option match_val will be the model confidence. The output in this example returned the correct material type, HDPE, as the top match. If using Pearson correlation, 0.7 is a good threshold to use for a positive ID. In this example, only our top match is greater than the threshold so we would disregard the other matches. If no matches were above our threshold, we would proclaim that the spectrum is of an unknown identity. You’ll also notice in this example that we matched to a library with both Raman and FTIR spectra but the Raman spectra had the highest hits, this is the rationale for lazily matching to a library with both. If you want to just match to a library with FTIR or Raman spectra, you can first filter the library using filter_spec() using SpectrumType.

matches <- match_spec(x = processed, library = test_lib, conform = T,
                      add_library_metadata = "sample_name", top_n = 5)
print(matches[,c("object_id", "library_id", "match_val", "SpectrumType",
                 "SpectrumIdentity")])

Library Metadata

The libraries we have created have over 100 variables of metadata in them and this can be onerous to read through especially given that many of the variables are NA values. We created get_metadata() to remedy this by removing columns from the metadata which are all blank values. The function below will return the metadata for the top match in matches. Remember, similar filter_spec(), you can specify logic for more than one thing at a time.

get_metadata(x = test_lib, logic = matches[[1,"library_id"]], rm_empty = T)

Plot Matches

Overlaying unknown spectra and the best matches can be extremely useful to identify peaks that don’t fit to the reference library which may need further investigation. The example below shows great correspondence between the best match and the unknown spectrum. All major peaks are accounted for and the correct relative height. There are two small peaks in the unknown spectrum near 500 that are not accounted for which could be investigated further but we would call this a positive id to HDPE.

plotly_spec(processed, filter_spec(test_lib, logic = matches[[1,"library_id"]]))

Sharing Reference Data

If you have reference data or AI models that you think would be useful for sharing with the spectroscopy community through OpenSpecy please contact the package administrator to discuss options for collaborating.

Brute Force

data("test_lib")
test_map <- read_any(read_extdata("CA_tiny_map.zip"))

test_map_processed <- process_spec(test_map, conform_spec_args = list(
  range = test_lib$wavenumber, res = NULL)
  )

identities <- match_spec(test_map_processed, test_lib, order = test_map,
                         add_library_metadata = "sample_name", top_n = 1)

features <- ifelse(identities$match_val > 0.7,
                   tolower(identities$polymer_class), "unknown")

id_map <- def_features(x = test_map_processed, features = features)

id_map$metadata$identities <- features
# Also should probably be implemented automatically in the function when a
# character value is provided. 

# Collapses spectra to their median for each particle
test_collapsed <- collapse_spec(id_map)

A Priori Particle Thresholding

data("test_lib")
test_map <- read_any(read_extdata("CA_tiny_map.zip"))

# Characterize the total signal as a threshold value. 
snr <- sig_noise(test_map,metric = "log_tot_sig")

# Use this to find your particles and the sig_noise value to use for
# thresholding. 
heatmap_spec(test_map, z = snr)

# Set define the feature regions based on the threshold. 400 appeared to be
# where I suspected my particle to be in the previous heatmap. 
id_map <- def_features(x = test_map, features = snr > 400) 

# Check that the thresholding worked as expected.
heatmap_spec(id_map, z = id_map$metadata$feature_id)

# Collapse the spectra to their medians based on the threshold. Important to
# note here that the particles with id -88 are anything from the FALSE values
# so they should be your background. 
collapsed_id_map <- id_map |>
  collapse_spec()

# Process the collapsed spectra. 
id_map_processed <- process_spec(collapsed_id_map, conform_spec_args = list(
  range = test_lib$wavenumber, res = NULL)
  )

# Check the spectra.
plot(id_map_processed)

# Get the matches of the collapsed spectra for the particles.
matches <- match_spec(id_map_processed, test_lib,
                      add_library_metadata = "sample_name", top_n = 1)