Analysis of olefins TICs from GC-MS with R

Angel Angelov | aangeloo@gmail.com | 2017

This script can be used to perform analysis of quantitative data from GC-MS runs from the Shimadzu instrument in R (including figure generation and statistics). The starting file is an “ASCIIData.txt” file, generated by the Shimadzu software, here is an example file you can use. It takes the peak areas for compounds generated by the Shimadzu software (GCMC Postrun) from the ASCIIData file and calculates the amounts of compounds based on an internal standard, in this case this is the C30 alkane peak area. Presumably several GCMS runs were done for each sample.

Load required libraries:

library(ggplot2)
library(dplyr)
library(reshape2)

Read in the text file. This file should be processed (with Excel may be) so that you have the samples (observations) in rows and the values (area of each olefin) in columns, tab separated. The command below opens a window so you can select your text file:

df <- read.csv(file.choose(), header = TRUE, sep = "\t", dec = ",")

A more programmatic approach to directly read the ASCIIData.txt file generated by the Shimadzu software is below. First off, it is a horrible format to clean and feed into R. Anyway, here is a dirty solution. It gets the number of runs in the file (using grep). The 8 rows skipped at the beginning and afterwards and the nrows = have to be changed for each ASCIIData file, depending on the particular GCMS run. Implemented in a responsive context (the ASCIIData file has to be in your working directory):

  # n <- as.integer(system("grep -c Header ASCIIData.txt", intern = TRUE)) #count samples (runs) in file
  n <- as.integer(grep("Header", readLines("ASCIIData.txt")) %>% length()) # without the system call above
  
  ole <- readline(prompt="Enter number of olefins detected for these runs(including C30): ")
11
  ole <- as.integer(ole)
  
  df <- data.frame() # make an empty data.frame
  skip <- seq(8, length.out = n, by = 20) # tell on which lines is the data, define what to skip in the next step
  for (i in skip) {
    df <- rbind(df,read.table("ASCIIData.txt", header = F, skip = i, sep="\t", nrows = ole))
    }
  df <- df %>% select(olefin = V2, area=V10) # take columns 2 and 10 (peak area)
  

Get the names of the samples in the ASCIIData file

names <- system("grep \"Data File Name\" ASCIIData.txt | cut -f 11 -d \"-\"", intern = TRUE)
tt <- list() # make an empty list
for (i in names) {tt <- cbind(tt,print(rep(i,ole)))} # repeat each sample name `ole` number of times
df$names <- paste(tt) # add sample names to df
#df$names <- substr(df$names,1,3) # in this case, take only sample name, removing replicate info
#df$names <- as.factor(df$names)
rm(tt)
df <- dcast(df,names~olefin, value.var = "area")

Now the ASCIIData.txt file is finally read and in the form I need it. Calculations and plotting:

C30conc <- readline(prompt="Enter the concentration of internal standard (µg/ml): ")
10
C30conc <- as.integer(C30conc)
df.norm <- df[,2:ncol(df)]/df$C30*C30conc # normalize peak areas to C30 peak
df.norm[,ncol(df.norm)] <- NULL # delete standard (last column)
df.norm$total <- rowSums(df.norm) # get total olefins per sample
df.norm$names <- df$names
df.norm$samples <- substr(df.norm$names, 6,8) # in this case the sample names are located in positions 6 to 8 in the names string
# first melt data, then get statistics and plot
df.norm.melt <- df.norm %>% select(-names) %>% melt(id.vars = "samples") # get rid of the "names" column here
df.norm.melt %>% group_by(samples, variable) %>% summarise(mean=mean(value),SD= sd(value),N= n()) %>% write.csv("out-quant-olefins.csv") # mean for each olefin and for total, grouped by sample
ggplot(df.norm.melt, aes(x= variable, y= value)) + geom_boxplot(aes(color=samples)) # and so on...