R/centroid.R

Defines functions centroid .determine.peakgroups .determine.peakgroups.orbitrap .is.peak .trapez .getProfileMS2

Documented in centroid

#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)
}
protViz/protViz documentation built on Jan. 19, 2024, 8:10 a.m.