Skip to content

Instantly share code, notes, and snippets.

@shackett

shackett/rNMR_shackett_custom.R

Last active March 18, 2016 20:12
Show Gist options
  • Save shackett/7950798 to your computer and use it in GitHub Desktop.
Save shackett/7950798 to your computer and use it in GitHub Desktop.
Download ZIP
Using rNMR, call this custom quantification metric which will subtract baseline and create informative plots of each specie before and after fit * First run shack_peakHeight or shack_peakArea as a custom ROI summary type* Then run shack_fitSummary(outList) in the R console to visualize fitting summary
Raw
rNMR_shackett_custom.R
shack_peakHeight <- function( inList ){
#### For baseline subtraction and quantification using peak height ####
# Save spectra baseline file in a list so that original and baseline
# subtracted spectra can be compared side-by-side
# Quantification by the maximum value once the baseline has been removed
require(baseline)
if(!exists("outList")){
outList <- list()
}
fName <- inList$file.par$file.name
w2 <- inList$w2
spectra <- baseline(t(as.matrix(inList$data)))
subOut <- list(name = fName, w2 = w2, spectra = spectra)
outList[[length(outList) + 1 ]] <- subOut
assign("outList", outList, inherits=FALSE, envir=.GlobalEnv)
return(max(getCorrected(spectra)))
}
shack_peakArea <- function( inList ){
#### For baseline subtraction and quantification using peak area ####
# Save spectra baseline file in a list so that original and baseline
# subtracted spectra can be compared side-by-side
# Quantification by sum of corrected spectra
require(baseline)
if(!exists("outList")){
outList <- list()
}
fName <- inList$file.par$file.name
w2 <- inList$w2
spectra <- baseline(t(as.matrix(inList$data)))
subOut <- list(name = fName, w2 = w2, spectra = spectra)
outList[[length(outList) + 1 ]] <- subOut
assign("outList", outList, inherits=FALSE, envir=.GlobalEnv)
return(sum(getCorrected(spectra)))
}
shack_fitSummary <- function( outList = outList ){
#### Generate summary plots comparing raw spectra to baseline subtracted spectra ####
# creates a "Fit_Summary" folder in the same directory as the sourced files split up by ROI #
require(reshape2)
require(ggplot2)
require(data.table)
scatter_facet_theme <- theme(text = element_text(size = 23, face = "bold"), title = element_text(size = 25, face = "bold"),
axis.text = element_text(color = 'white'), axis.text.x = element_text(angle = 90), strip.text = element_text(color = "coral1"),
legend.title = element_text(color = 'white'), legend.text = element_text(color = 'white'),
legend.position = "right", panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
panel.background = element_rect(fill = 'black'), plot.background = element_rect(fill = 'black'),
strip.background = element_rect(fill = "cadetblue2"), legend.background = element_rect(fill = 'gray50'),
legend.key = element_rect(fill = 'black')
)
# summarize each quantified peak - file name, directory, ROI name
peakInfo <- data.frame(index = 1:length(outList),
fullPath = sapply(outList, function(x){x$name}),
w2min = sapply(outList, function(x){min(x$w2)}),
w2max = sapply(outList, function(x){max(x$w2)}), stringsAsFactors = F)
peakInfo$ROI <- mapply(function(wmin,wmax){roiTable$Name[which.min(abs(wmin - roiTable$w2_downfield) + abs(wmax - roiTable$w2_upfield))]},
wmin = peakInfo$w2min, wmax = peakInfo$w2max)
peakInfo$dirname = sapply(peakInfo$fullPath, dirname)
peakInfo$filename = sapply(peakInfo$fullPath, function(x){strsplit(basename(x), split = '.ucsf')[[1]][1]})
# analyze each directory containing spectra seperately
for(a_directory in unique(peakInfo$dirname)){
peakInfo_subset <- peakInfo[peakInfo$dirname == a_directory,]
# check to see if a plot folder exists
if(!file.exists(file.path(a_directory, "NMRplots"))){
dir.create(file.path(a_directory, "NMRplots"))
}
# setup and generate a plot for each ROI seperately
for(a_ROI in unique(peakInfo_subset$ROI)){
ROI_subset <- peakInfo_subset[peakInfo_subset$ROI == a_ROI, c('index', 'filename')]
ROI_stack <- do.call(rbind, lapply(c(1:nrow(ROI_subset)), function(a_row){
data.frame(sample = ROI_subset$filename[a_row], w2 = outList[[ROI_subset $index[a_row]]]$w2,
raw = t(getSpectra(outList[[ROI_subset $index[a_row]]]$spectra)),
corrected = t(getCorrected(outList[[ROI_subset $index[a_row]]]$spectra)))
}))
ROI_stack <- data.table(melt(ROI_stack, id.vars = c("sample", "w2"), variable.name = "class", value.name = "intensity"))
ROI_stack$sample <- factor(ROI_stack$sample); ROI_stack$class <- factor(ROI_stack$class)
ROI_peak <- ROI_stack[,list(intensity = max(intensity), w2 = w2[which.max(intensity)]), by = c("sample", "class")]
ROI_summary <- ROI_peak[,list(text = paste(signif((intensity[class == "raw"] - intensity[class == "corrected"])/intensity[class == "raw"]*100, 0), "%", sep = ""), height = mean(intensity), w2 = w2[class == "corrected"]), by = "sample"]
ROIplot <- ggplot() + facet_wrap(~sample, scale = "free_y") +
geom_line(data = ROI_stack, aes(x = w2, y = intensity, color = class), size = 0.5) + geom_hline(yintercept = 0, color = "darkgoldenrod1") +
geom_segment(data = ROI_peak, aes(x = w2, xend = w2, y = 0, yend = intensity), size = 1, color = "white") +
geom_text(data = ROI_summary, aes(x = w2, y = height, label = text), color = "white", hjust = -0.15) +
scatter_facet_theme + expand_limits(y = 0) + scale_color_brewer(palette = "Set1")
ggsave(file = file.path(a_directory, "NMRplots", paste(a_ROI, ".pdf", sep = "")), plot = ROIplot, height = 6 + sqrt(nrow(ROI_subset))*3, width = 8 + sqrt(nrow(ROI_subset))*4, limitsize = F)
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment