#R
################################################################################
# Christian Panse <cp@fgcz.ethz.ch>
# 2020-04-14
# see also:
# https://github.com/lgatto/MSnbase/blob/2b6d6f5162e7464c39b48d0daa4d981c6ec30bbd/R/functions-Spectrum.R#L649
# Orbitrap Fusion Lumos FSN20242
# p2722/.../.../stds_pos_neg_MS_highconc_UVPD_50_300.raw
# scan 1959
.getProfileMS2 <- function(){
intensity <- c(0, 0, 0, 0, 0, 0, 0, 0, 4592.90478515625,
7091.9931640625, 8611.9755859375, 8369.5830078125, 5866.6279296875,
0, 0, 0, 0, 0, 0, 0, 0, 3156.986328125, 7038.2119140625,
9386.5810546875, 8286.3486328125, 4776.27978515625, 0, 0, 0, 0, 0, 0,
0, 0, 0, 5108.396484375, 5117.8349609375, 84061.6640625, 118886.5625,
320283.5, 1089958.375, 3251054.75, 5548869, 7183389.5, 5562785,
3252974.5, 1067703, 279433.28125, 88928.6015625, 78262.015625,
18332.697265625, 13843.185546875, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1030.65051269531, 4999.087890625, 13393.734375,
21808.88671875, 34574.78515625, 18202.544921875, 9088.7490234375,
2902.12182617188, 0, 0, 0, 0, 0, 0, 0, 0, 3872.65942382812,
7224.23388671875, 8925.8095703125, 7839.21875, 4900.6533203125, 0, 0,
0, 0, 0, 0, 0, 0, 2410.578125, 6061.23779296875, 19732.857421875,
35713.9296875, 46617.46875, 35048.6953125, 19386.44921875,
6379.92822265625, 2889.75024414062, 0, 0, 0, 0, 0, 0, 0, 0,
2862.50463867188, 6281.21435546875, 17308.380859375, 33795.66796875,
46966.4921875, 37150.4921875, 21012.11328125, 7166.1201171875,
2516.72485351562, 0, 0, 0, 0, 0, 0, 0, 0, 4881.4619140625,
5934.54150390625, 25561.189453125, 52422.23046875, 80021.7578125,
76278.0078125, 47473.97265625, 20932.923828125, 4777.568359375, 0, 0,
0, 0, 0, 0, 0, 0, 7265.79638671875, 14160.1884765625,
17088.591796875, 13269.5, 6262.5625, 0, 6208.72900390625,
12367.2783203125, 15205.4873046875, 11901.078125, 5150.2880859375, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9042.5947265625, 21935.7734375,
28896.01171875, 15938.279296875, 6433.23876953125, 4853.4267578125,
11348.7666015625, 20853.142578125, 23384.671875, 16632.787109375,
6872.1298828125, 0, 0, 0, 0, 0, 0, 0, 5160.78466796875,
11936.0849609375, 14292.1318359375, 10230.9365234375,
4040.53955078125, 0, 0, 0, 0, 0, 0, 0, 0, 4704.51513671875,
12345.115234375, 28885.31640625, 40221.6640625, 32290.67578125,
18112.1015625, 6414.52392578125, 0, 0, 0, 0, 0, 0, 0, 0,
2931.01000976562, 7218.38330078125, 10585.7294921875,
10754.474609375, 24779.36328125, 97241.171875, 316750.1875,
563389.625, 762858.0625, 620968.4375, 369190.3125, 136320.453125,
40065.2890625, 15975.5263671875, 12490.22265625, 3184.92211914062, 0,
0, 0, 0, 0, 0, 0, 0, 4765.2080078125, 10852.9716796875,
18578.505859375, 24274.90625, 23566.21484375, 16251.158203125,
7466.9970703125, 0, 0, 0, 0, 0, 0, 0, 0, 2977.11328125,
9619.7529296875, 26667.041015625, 43905.8515625, 45744.2109375,
29386.36328125, 14009.4697265625, 4155.806640625, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13210.6142578125,
8176.30517578125, 0, 0, 0, 0, 0, 0, 0, 0, 8506.875, 11196.6259765625,
12399.216796875, 11524.4951171875, 8141.75537109375, 0, 0, 0, 0, 0,
0, 0, 0, 0, 7342.37109375, 33556.92578125, 70206.71875, 108179.5625,
102746.5, 63951.015625, 29043.640625, 6719.19384765625, 0, 0, 0, 0,
0, 0, 0, 0, 4147.06005859375, 17146.6640625, 56676.96875,
97892.1640625, 126585.9609375, 97967.5234375, 57022.9453125,
17630.9765625, 3792.12329101562, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 80385.0703125, 306936.90625, 342259.71875,
864612.4375, 3378294.75, 10811007, 19051660, 25577258, 20697266,
12314792, 4545463, 1350920, 530626.875, 403116.25, 115720.609375,
20698.126953125, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
2686.4072265625, 18135.4375, 83746.125, 245605.015625, 416970.96875,
526961.125, 400583.1875, 231464.40625, 76335.2421875, 21159.19921875,
6607.9921875, 0, 0, 0, 0, 0, 0, 0, 0)
mZ <- c(99.0035825901103,
99.0046403943294, 99.0056982155017, 99.0067560536276,
104.620208956094, 104.621358041487, 104.622507145812,
104.623656269068, 104.624805411256, 104.625954572377,
104.627103752431, 104.628252951418, 104.62940216934,
104.630555142021, 104.631704397811, 104.632853672536,
104.634002966197, 109.365155510419, 109.366383649002,
109.367611808273, 109.368839988231, 109.370068188878,
109.371296410214, 109.372524652239, 109.373752914954,
109.37498119836, 109.376202821526, 109.377431146314,
109.378659491793, 109.379887857964, 112.033842424927,
112.035115789622, 112.036389176027, 112.037662584142,
112.038936013967, 112.040209465503, 112.04148293875, 112.04275643371,
112.044029950382, 112.045303488767, 112.046577048865,
112.047850630677, 112.049124234203, 112.050397859445,
112.051671506402, 112.052945175074, 112.054218865463,
112.055492577569, 112.056766311393, 112.058040066934,
112.059313844193, 112.060587643171, 112.061861463869,
112.06313527078, 112.064409134918, 112.065683020776,
112.066956917471, 112.068230846773, 112.069504797797,
112.070778770544, 112.07205030673, 112.073324322924,
112.074598360842, 112.075872420485, 113.04402324308,
113.045313869001, 113.046604517024, 113.04789518715, 113.04918587938,
113.050476593714, 113.051767330153, 113.053058088697,
113.054348869347, 113.055639672103, 113.056930496966,
113.058221343936, 113.059527481029, 113.060818372215,
113.06210928551, 113.063400220915, 113.711647730784,
113.712949806734, 113.71425190505, 113.71555402573, 113.716856168776,
113.718158334188, 113.719460521966, 113.720762732112,
113.722064964625, 113.723369748746, 113.724672025996,
113.725974325615, 113.727276647604, 128.997844460787,
128.99941772846, 129.000991024914, 129.002564350151, 129.00413770417,
129.005711086974, 129.007284498562, 129.008857938935,
129.010431408094, 129.012004906039, 129.013578432772,
129.015151988293, 129.016725572602, 129.018306923798,
129.019880565687, 129.021454236367, 129.023027935838,
135.030417436874, 135.032102345484, 135.033787285631,
135.035472257316, 135.037157260538, 135.0388422953, 135.040527361601,
135.042212459442, 135.043897588825, 135.04558274975,
135.047267942218, 135.04895316623, 135.050638421785,
135.052333359397, 135.054018678044, 135.055704028238,
135.057389409979, 140.995162029413, 140.996959804026,
140.998757613022, 141.000555456404, 141.002353334172,
141.004151246327, 141.005949192869, 141.007747173801,
141.009545189121, 141.011343238832, 141.013141322935,
141.014939441429, 141.016737594317, 141.018546394251,
141.020344615928, 141.022142871999, 141.023941162468,
147.471669798511, 147.473592854024, 147.475515947153,
147.477439077899, 147.479362246262, 147.481285452244,
147.483208695845, 147.485131977067, 147.48705529591,
147.488978633779, 147.490902027868, 147.492825459582,
147.494748928922, 147.496672435887, 147.498595980481,
147.500519543494, 147.502443163345, 147.504366820826,
147.506290515939, 147.508214248684, 147.510138019062,
147.512061794274, 147.513985639922, 147.515909523206,
147.517833444128, 147.519757402687, 147.521681398886,
147.523605432725, 147.525529504205, 147.527453613327,
147.529377760092, 147.531301944502, 147.533226148294,
147.535150407994, 147.537074705342, 147.538999040338,
147.540923412982, 147.542847823277, 147.544772233941,
147.546696719538, 147.548621242789, 147.550545803693,
147.552470402253, 147.554395038468, 147.55631971234, 147.55824442387,
147.560169173059, 147.562093959907, 147.564018784417,
147.565943646588, 147.567881589027, 147.569806526525,
147.571731501687, 147.573656514515, 150.088643255271,
150.090617725709, 150.09259223511, 150.094566783475,
150.096541370804, 150.098515997098, 150.100490662359,
150.102465366588, 150.104440109785, 150.106414891952,
150.10838971309, 150.110376097007, 150.112350996089,
150.114325934145, 150.116300911176, 168.084792409383,
168.087132440456, 168.089472520395, 168.091812649201,
168.094152826876, 168.096493053422, 168.098833328839,
168.10117365313, 168.103514026294, 168.105854448335,
168.108194919252, 168.110535439048, 168.112876007724,
168.11521662528, 168.11755729172, 168.119898007043, 168.122238771252,
168.124579584347, 168.12692044633, 168.129261357202,
168.131601146235, 168.13394215489, 168.136283212439,
168.138624318882, 169.024717645059, 169.027077331693,
169.029437067741, 169.031796853205, 169.034156688086,
169.036516572385, 169.038876506104, 169.041236489244,
169.043596521807, 169.045956603793, 169.048316735205,
169.050676817338, 169.053037047604, 169.0553973273, 169.057757656427,
169.102613303524, 169.104974621557, 169.10733598905,
169.109697406004, 169.11205887242, 169.1144203883, 169.116781953646,
169.119143568458, 169.121505232739, 169.123866946489,
169.12622870971, 169.128590522403, 169.130961226008,
169.133323137649, 169.135685098767, 169.138047109363,
178.200216622696, 178.202771037331, 178.205325506891,
178.207880031377, 178.210434610791, 178.212989245135,
178.21554393441, 178.218098678617, 178.220653477759,
178.223208331837, 178.225763240852, 178.228318204806, 178.2308732237,
178.233428297537, 178.235983410106, 178.238538593832,
178.241176891941, 178.243732185561, 178.246287534132,
178.248842937655, 226.861659683684, 226.865328879715,
226.868998164762, 226.87266753883, 226.876337001922,
226.880006554039, 226.883676195185, 226.887345925364,
226.891015744577, 226.894685546307, 226.898355543598,
226.902025629933, 226.905695805314, 227.052575871578,
227.056249700338, 227.059923618265, 227.063597625362,
227.067271721633, 227.070945907079, 227.074620181705,
227.078294545512, 227.081968998504, 227.085643540684,
227.089318172054, 227.092992892617, 227.096667702377,
227.100474751961, 227.104149740121, 227.107824817487,
227.11149998406, 261.11116858176, 261.11569930404, 261.120230144245,
261.124761102379, 261.129292178445, 261.133823372448,
261.138354684393, 261.142886114282, 261.147417662121,
261.151949327913, 261.156481111662, 261.161013013373,
261.16554503305, 261.17006052747, 261.174592783089, 261.179125156687,
261.183657648266, 296.006212936471, 296.011681596413,
296.017150407905, 296.022619370952, 296.028088485561,
296.033557751737, 296.039027169485, 296.044496738812,
296.049966459722, 296.055436330667, 296.060906354762,
296.066376530456, 296.071846857758, 296.07731733667,
296.082787967201, 296.088258749354, 296.093729683135,
296.099200768551, 296.104672005607, 296.110143394307,
296.115614934659, 296.121086626668, 296.126558470338,
296.132030465677, 296.137502612688, 296.142974911379,
296.148447361754, 296.15391996382, 296.159392717581,
296.164865623044, 296.170338679011, 296.175811887894,
296.181285248494, 296.186758760819, 296.192232424512,
296.197706240301, 296.20318020783, 296.208654308923,
296.214128579951, 296.219603002736, 296.225077577284,
297.069976837928, 297.075475003739, 297.080973322191,
297.086471793289, 297.091970417041, 297.097469193451,
297.102968122524, 297.108467204268, 297.113966438686,
297.119465825786, 297.124965365572, 297.13046505805,
297.135964903226, 297.141464901106, 297.146965051695,
297.152465354998, 297.157965811023, 297.163466419773,
297.168967181255, 310.095889034636, 310.101752762872,
310.107616657429, 310.113480718314)
idx <- order(mZ)
list(mZ=mZ[idx], intensity=intensity[idx])
}
# https://de.wikipedia.org/wiki/Trapez-Methode
.trapez <- function(x, y){
idx <- 2:length(x)
((x[idx] - x[idx-1]) %*% (y[idx] + y[idx-1])) / 2
}
# https://github.com/cpanse/p2135R/blob/55eab2376d52e35da127dc3eb74924a2acd9e400/vignettes/7_RSintensityCmp_20181221.Rmd#L160
.is.peak <- function(idx, x){
# criteria for not being a peak
if (idx == 1){return(FALSE)}
if (idx == length(x)){return(FALSE)}
if (x[idx - 1] > x[idx]){return(FALSE)}
if (x[idx + 1] > x[idx]){return(FALSE)}
TRUE
}
# splits profile peakgroups by 0 intensities and ignores peak shapes
.determine.peakgroups.orbitrap <- function(x){
V <- x != 0.0;
n <- length(V);
peakgroup <- rep(0, n)
count <- 0;
for (i in 1:n){
if (V[i]){
if(!V[i-1]){count <- count + 1};
peakgroup[i] <- count}
}
peakgroup
}
# simple heuristic
#
# Note: this method will not detect overlapping peaks as it
# can be done by IMSTOF http://www.tofwerk.com/ libraries.
#
# TODO(cp): if the while loops are to slow replace it by some Rcpp constructs
.determine.peakgroups <- function(mZ, x, tolppm){
peak.idx <- which(sapply(1:length(x), .is.peak, x=x))
n <- length(x)
peakgrps <- rep(0, n)
count <- 1
for (i in peak.idx){
# mid
idx <- i
peakgrps[i] <- count
eps <- 1e-06 * tolppm * mZ[idx]
#lower
while (x[idx - 1] < x[idx] & idx > 2 & abs(mZ[idx - 1] - mZ[idx]) < eps){
peakgrps[idx] <- count
idx <- idx - 1
}
#upper
idx <- i
while (x[idx + 1] < x[idx] & idx < n & abs(mZ[idx + 1] - mZ[idx]) < eps){
peakgrps[idx] <- count
idx <- idx + 1
}
count <- count + 1
}
peakgrps
}
centroid <- function(mZ, intensity, tolppm=100, debug=FALSE){
stopifnot(length(mZ) == length(intensity))
if (debug){
plot((diff(mZ)) ~ mZ[2:length(mZ)], log='y');
abline(h = 0.1, col='red')
points(mZ , tolppm * 1e-06 * mZ, col='green', type='l')
}
n <- length(mZ)
# peakgrps <- c(0, cumsum(diff(mZ) > eps))
peakgrps <- split(1:n, .determine.peakgroups(mZ, intensity, tolppm))
peakgrps <- peakgrps[names(peakgrps) != '0']
rv <- lapply(peakgrps , FUN=function(i){
if(length(i) > 2){
intensity.auc <- .trapez(mZ[i], intensity[i])
mZ.centroid <- weighted.mean(x=mZ[i], w=intensity[i])
# TODO(cp): replace by WEW's f: x -> ax^2+bx+c
# interpolation using the three highes peaks to derive AUC
if (debug){
plot(diff(mZ[i]))
plot(mZ[i], intensity[i],
type='h',
xlab='mZ',
ylab='intensity',
ylim = c(0, max(intensity[i])))
abline(v=mZ.centroid, col='red')
points(mZ.centroid, intensity.auc,
type='h', col='red', lwd=4)
legend("topright", NULL, round(intensity.auc, 3), title="AUC")
axis(3, mZ.centroid, mZ.centroid)
}
return (data.frame(mZ=mZ.centroid, intensity=intensity.auc))
}
NULL
})
do.call('rbind', rv)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.