library(MSnbase)
Load 6 of the CDF files from the faahKO
cdf_files <- dir(system.file("cdf", package = "faahKO"), recursive = TRUE,
full.names = TRUE)[c(1:3, 7:9)]
Define the sample grouping.
s_groups <- rep("KO", length(cdf_files))
s_groups[grep(cdf_files, pattern = "WT")] <- "WT"
s_groups
## [1] "KO" "KO" "KO" "WT" "WT" "WT"
Define a data.frame that will be used as phenodata
pheno <- data.frame(sample_name = sub(basename(cdf_files), pattern = ".CDF",
replacement = "", fixed = TRUE),
sample_group = s_groups, stringsAsFactors = FALSE)
pheno
## sample_name sample_group
## 1 ko15 KO
## 2 ko16 KO
## 3 ko18 KO
## 4 wt15 WT
## 5 wt16 WT
## 6 wt18 WT
Read the data and extract chromatograms
raw_data <- readMSData2(cdf_files, pdata = new("NAnnotatedDataFrame", pheno))
## Polarity can not be extracted from netCDF files, please set manually the polarity with the 'polarity' method.
## Warning in nrow(fullhd): bytecode version mismatch; using eval
# here we do some arbitrary mz and rt intervals to show how that would work
# This does not need to be supplied to get TICs and BPIs
rtr <- rbind(c(45*60, 55*60), c(65*60, 70*60))
mzr <- rbind(c(200, 220), c(300, 320))
bpis <- chromatogram(raw_data, rt = rtr, mz = mzr, aggregationFun = "max")
bpis
## Chromatograms with 2 rows and 6 columns
## ko15.CDF ko16.CDF ko18.CDF wt15.CDF
## <Chromatogram> <Chromatogram> <Chromatogram> <Chromatogram>
## [1,] length: 384 length: 384 length: 384 length: 384
## [2,] length: 192 length: 192 length: 192 length: 192
## wt16.CDF wt18.CDF
## <Chromatogram> <Chromatogram>
## [1,] length: 384 length: 384
## [2,] length: 192 length: 192
library(RColorBrewer)
sample_colors <- brewer.pal(3, "Set1")[1:2]
names(sample_colors) <- c("KO", "WT")
sample_colors
## KO WT
## "#E41A1C" "#377EB8"
col <- paste0(sample_colors[raw_data$sample_group], "80")
col
## [1] "#E41A1C80" "#E41A1C80" "#E41A1C80" "#377EB880" "#377EB880" "#377EB880"
The “base” plot using the plotting feature from MSnbase.
plot(bpis, col = col)
.ChromatogramLongFormat <- function(x){
data.frame(rtime = x@rtime,
intensity = x@intensity,
spectrum_id = names(x@rtime),
fromFile = x@fromFile,
mz_interval = paste0(x@filterMz[1]," - ",paste0(x@filterMz[2])),
row.names = NULL,
stringsAsFactors = FALSE
)
}
.ChromatogramsLongFormat <- function(x, pdata = NULL){
table <- lapply(x@.Data, .ChromatogramLongFormat)
table <- do.call(rbind.data.frame,table)
table$mz_interval <- as.factor(table$mz_interval)
# add filename
filename <- cbind.data.frame(fromFile = 1:attributes(x)$dim[2], filename = attributes(x)$dimnames[[2]], stringsAsFactors = FALSE)
table <- merge(table, filename, all.y = TRUE, by = "fromFile")
# if we supplied pdata add the corresponding annotation to each row
if(!is.null(pdata)){
# we merge the pdata columns to be able to supply a text column that can be shown in plotly tooltips
text <- lapply(names(pdata@data), function(x) paste0(x,": " , as.character(as.matrix(pdata@data[x])),"<br>"))
text <- apply(as.data.frame(text),1,paste, collapse="")
pdata <- cbind.data.frame(fromFile = 1:nrow(pdata), pdata@data, text = text, stringsAsFactors = FALSE)
table <- merge(table, pdata, all.y = TRUE, by = "fromFile")
}
return(table)
}
ggplotChromatograms <- function(bpis, rt_unit = "sec", pdata = NULL) {
require(ggplot2)
bpis_long <- .ChromatogramsLongFormat(bpis, pdata = pdata)
# show in sec or min
if(rt_unit=="min"){
bpis_long$rtime <- bpis_long$rtime/60
x_lab <- "Retention time (min)"
}
if(rt_unit=="sec"){
x_lab <- "Retention time (sec)"
}
# if pdate is supplied add it otherwise only the spectrum_id and filename goes to the text field
if(!is.null(pdata)){
bpis_long <- transform(bpis_long, text=paste0("filename: ", filename, "<br>","spectrum_id: ",spectrum_id,"<br>",text))
}else{
bpis_long <- transform(bpis_long, text=paste0("filename: ", filename, "<br>","spectrum_id: ", spectrum_id))
}
# Prettier name for RT
names(bpis_long)[names(bpis_long) == 'rtime'] <- 'Retention time'
# actual plot
p <- ggplot(data = bpis_long, aes(x = `Retention time`, y = intensity, text = text, group = interaction(fromFile,mz_interval))) +
geom_line() +
theme_classic() +
labs(x = x_lab, y = "Intensity")
return(p)
}
ggplotChromatograms(bpis)
## Loading required package: ggplot2
p <- ggplotChromatograms(bpis, rt_unit = "min")
p
… assigning a new color to each sample using default ggplot2 color
p + aes(color=as.factor(fromFile))
This does the same as the original base plot
names(col) <- 1:length(col)
p + aes(color=as.factor(fromFile)) +
scale_colour_manual(name = "File", values = col)
We can instead add some info about each sample.
We can then use the columns in this table to modify the plot.
p <- ggplotChromatograms(bpis, rt_unit = "min", pdata = new("NAnnotatedDataFrame", pheno)) +
aes(color=sample_group)
p
p + scale_colour_manual(name = "Group", values = sample_colors)
# Interactivity with plotly
We can easily make an interactive plot now.
Notice the hover over text.
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
p <- p +
scale_colour_manual(name = "Group", values = sample_colors)
ggplotly(p, dynamicTicks = TRUE)
You can also split by the mz interval group.
p <- p + facet_wrap(~mz_interval, scales = "free")
ggplotly(p, dynamicTicks = TRUE)