Peakgrouping
Usually one of the first things you will want to do with such data is to produce peak-groups of some kind – grouping peaks with similar mass to identify them with each other. This package supplies two functions for this purpose:
extract_masses, anddbscan.
After introducing these two functions and their uses, I’ll introduce a some heuristics I often use for judging the quality of a peak-group clustering, particularly applicable to peak-grouping by dbscan, but that can also be applied peak-grouping done by extract_masses.
Extract Masses
Most simply, extract_masses can be used to simply extract peaks near known masses of interest, in this dataset for example, m/z = 1570.677 is a mass of interest, and so we could extract it as:
df.cal = data.frame(mass = c(1296.685,
1570.677,
2147.199,
2932.588),
name = c("Angiotensin_I",
"Glu_Fibrinopeptide_B",
"Dynorphin_A",
"ACTH_fragment1_24"),
short = c("AngI",
"GluFib",
"DynA",
"ACTH"))
ppm_err = 20
pl.cal = extract_masses(pl, df.cal$mass, margin = ppm_err, use_ppm = TRUE)This returns pl.cal as a simple subset of pl of peaks within 20 ppm of the m/z values specified in df.cal$mass, with an additional column called group whose values are the values of df.cal$mass to which each peak is matched.
DBSCAN
As mentioned in the documentation, dbscan provides a univariate optimisation of the DBSCAN* algorithm of section 3 of Campello et al. (2013). This is a density based clustering algorithm that, breifly, can be thought of as two steps performed in sequence, first points with insufficient density are removed, and second the remaining points are grouped as per the equivelance classes induced by an equivalence relation. In this context density is measured as the number of points within a certain eps, and the equivalence relation used is “are within eps of each other”. This second step, without the first, is often called tolerance clustering and is often sufficient – this can be acheived by specifing a minimum density of 1:
pl$tol_clus_1 = dbscan(pl$m.z, eps = 0.1, mnpts = 1)
pl$tol_clus_2 = dbscan(pl$m.z, eps = 0.05, mnpts = 1)But more generally the minimum density, mnpts, will need to be adjusted based on the dataset, usually proportionally to the size of the dataset but also depending on other factors. Similarly eps will depend on the dataset, and is directly related to mnpts – the same value of mnpts is a more stringent density cutoff for a smaller value of eps. This test dataset is quite small, so a small value of mnpts is appropriate, such as for example:
pl$dbs_clus_1 = dbscan(pl$m.z, eps = 0.1, mnpts = 100)
pl$dbs_clus_2 = dbscan(pl$m.z, eps = 0.05, mnpts = 100)
pl$dbs_clus_3 = dbscan(pl$m.z, eps = 0.1, mnpts = 250)
pl$dbs_clus_4 = dbscan(pl$m.z, eps = 0.05, mnpts = 250)Note I have generated several different clusterings by varying the parameters, and will compare these in the following section.
Heuristics for diagnosis of peak-grouping.
Here I’ll show how to construct a plot I often use for heuristically diagnosing the quality of the performed peakpicking, and to adjust the tuning parameter selection to get a peak-grouping i am happy with. First, a couple of helper functions for calculating peak-group summaries, and generating the plot to which I alluded respectively:
summarise_peakgroups = function(pl, var) {
return(ddply(pl, var, summarise,
AWM = weighted.mean(m.z, log1p(intensity)),
Centre = (max(m.z) + min(m.z)) / 2,
Range = max(m.z) - min(m.z),
nPeaks = length(m.z),
nDupPeaks = (length(m.z) - length(unique(Acq)))))
}
diag_plot = function(pgs, t) {
return(ggplot(pgs, aes(x=Range, y=nDupPeaks))
+ ggtitle(t)
+ geom_jitter(width = 0, height = 0.4)
+ ylab("# Duplicate Peaks"))
}Given these functions, peak-group summaries can be calculates as:
pgs_tol_1 <- summarise_peakgroups(pl, "tol_clus_1")
pgs_tol_2 <- summarise_peakgroups(pl, "tol_clus_2")
pgs_dbs_1 <- summarise_peakgroups(pl, "dbs_clus_1")
pgs_dbs_2 <- summarise_peakgroups(pl, "dbs_clus_2")
pgs_dbs_3 <- summarise_peakgroups(pl, "dbs_clus_3")
pgs_dbs_4 <- summarise_peakgroups(pl, "dbs_clus_4")and heuristic diagnostic plots can then be generated like so:
print(diag_plot(subset(pgs_tol_1, tol_clus_1 != 0),
"eps: 0.1, mnpts: 1\n0% of peaks unassigned"))
print(diag_plot(subset(pgs_tol_1, tol_clus_1 != 0),
"eps: 0.05, mnpts: 1\n0% of peaks unassigned"))
per.unassigned_1 = round(100 * subset(pgs_dbs_1, dbs_clus_1 == 0)$nPeaks / nrow(pl), 1)
print(diag_plot(subset(pgs_dbs_1, dbs_clus_1 != 0),
paste0("eps: 0.1, mnpts: 100\n", per.unassigned_1,
"% of peaks unassigned")))
per.unassigned_2 = round(100 * subset(pgs_dbs_2, dbs_clus_2 == 0)$nPeaks / nrow(pl), 1)
print(diag_plot(subset(pgs_dbs_2, dbs_clus_2 != 0),
paste0("eps: 0.05, mnpts: 100\n", per.unassigned_2,
"% of peaks unassigned")))
per.unassigned_3 = round(100 * subset(pgs_dbs_3, dbs_clus_3 == 0)$nPeaks / nrow(pl), 1)
print(diag_plot(subset(pgs_dbs_3, dbs_clus_3 != 0),
paste0("eps: 0.1, mnpts: 250\n", per.unassigned_3,
"% of peaks unassigned")))
per.unassigned_4 = round(100 * subset(pgs_dbs_4, dbs_clus_4 == 0)$nPeaks / nrow(pl), 1)
print(diag_plot(subset(pgs_dbs_4, dbs_clus_4 != 0),
paste0("eps: 0.05, mnpts: 250\n", per.unassigned_4,
"% of peaks unassigned")))
Duplicate peaks here refers to multiple peaks from the same spectrum occuring in the same peak-group. Usually, we would not expect this to happen if the peak-picking algorithm is working well (although occasional artefacts can happen). Either way, number of duplicate peaks should be zero or if that is not possible, then it should be as small as possible (and noted, as these duplicates will likely have to be dealt with somehow later in the processing pipeline, either by averaging their intensities or by some other method, depending on the outcome being pursued).
I jitter the points in these plots vertically by +/- 0.2 so that in the case that there are no duplicate peaks, or close to, there will be some visual que for the density of these otherwise overlapping points, giving an idea what kind of range is typical.
Peaks left as unassigned due to insufficient density in the first step of the DBSCAN* algorithm would usually be discarded from this point downstream, so to avoid throwing out important information the number of peaks left unassigned during the peak-grouping step should generally be minimised if possible. On the other hand, if a peak-grouping that acheives zero duplicate peaks in all clusters can be acheived this is preferable as otherwise duplicate peaks would generally have to be dealt with by some exclusion or averaging method in order to avoid ambiguity. On the other hand if there are only a handful of clusters with a number of duplicate peaks greater than zero, these could be manually inspected and a decision be made as to what to do about them. For example, lets consider an example – the clustering performed above with eps: 0.1 and mnpts: 250.
print(subset(pgs_dbs_3, nDupPeaks > 0 & dbs_clus_3 != 0), digits = 6)## dbs_clus_3 AWM Centre Range nPeaks nDupPeaks
## 18 17 1081.52 1081.44 0.481 4421 445
## 25 24 1101.63 1101.72 0.743 4893 481
## 245 244 1592.72 1592.73 0.319 1185 1
## 280 279 1665.96 1665.91 0.336 2233 4
## 372 371 1850.00 1850.00 0.425 1287 10This allows us to identify the mass range in which these peaks occur, and other information about these peakgroups. To go even deeper, lets take a closer look at the m/z distribution of a couple of these peakgroups, specifically the two with over 400 duplicate peaks which seem to be particularly problematic:
pgs = subset(pgs_dbs_3, nDupPeaks > 400 & dbs_clus_3 != 0)$dbs_clus_3
centres = pgs_dbs_3[pgs_dbs_3$dbs_clus_3 %in% pgs, "Centre"]
margins = pgs_dbs_3[pgs_dbs_3$dbs_clus_3 %in% pgs, "Range"] / 2
print(qplot(subset(pl, dbs_clus_3 == pgs[1])$m.z, binwidth = 0.01)
+ ggtitle(paste(round(centres[1], 2), "+/-", round(margins[1], 2)))
+ xlab("m/z"))
print(qplot(subset(pl, dbs_clus_3 == pgs[2])$m.z, binwidth = 0.01)
+ ggtitle(paste(round(centres[2], 2), "+/-", round(margins[2], 2)))
+ xlab("m/z"))
Now lets take a look at how these peaks are clustering by the more stringent clustering with eps: 0.05 and mnpts: 250:
pl.sub = subset(pl, abs(m.z - centres[1]) <= margins[1])
print(qplot(pl.sub$m.z, fill = factor(pl.sub$dbs_clus_4), binwidth = 0.01)
+ ggtitle(paste(round(centres[1], 2), "+/-", round(margins[1], 2)))
+ guides(fill = FALSE) + xlab("m/z"))
pl.sub = subset(pl, abs(m.z - centres[2]) <= margins[2])
print(qplot(pl.sub$m.z, fill = factor(pl.sub$dbs_clus_4), binwidth = 0.01)
+ ggtitle(paste(round(centres[2], 2), "+/-", round(margins[2], 2)))
+ guides(fill = FALSE) + xlab("m/z"))
This looks much better – note that in the figures above colours denote different clusters, with the salmon colour representing peakgroup zero, which is the unassigned peaks due to insufficient density. Of interest is the two seemingly seperate green groups of peaks that are actually clustered together, but as seen previously, do not contain any duplicate peaks. Meaning the two consist of non-overlapping sets of spectra. This is interesting as it could potentially represent a systematic mass error is a portion of the spectra in this dataset, and could be further investigated as such. Particularly as this is peptide data, different resolvable masses are expected to differ by roughly 1Da, and the distance between these groups is much less than that.
