/
peakShape.jl
78 lines (64 loc) · 2.68 KB
/
peakShape.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
push!(LOAD_PATH, pwd())
import PyPlot
import HDF5
import InterpolationFunctions
import BaselineFunctions
import PeakshapeFunctions
function baselineAndPeakshape(
filepath;
outputfilename="results/_result.hdf5",
peakshapeRegions=8,
peakshapeRegionStretch=0.5,
peakshapeQuantileValue = 0.1,
peakfindingNoiseThresholdValue = 25,
peakWindowWidth = 200,
peakfindingSignalLimit = 0.2
)
file = joinpath(filepath, outputfilename)
massAxis = HDF5.h5read(file, "MassAxis")
avgSpectrum = HDF5.h5read(file, "AvgSpectrum")
baselinePoints, baselineValues, baselineNoise = BaselineFunctions.calculateBaseline(massAxis,avgSpectrum,baselinePointWidth = 0.8, threshold=0.2)
baselineNoiseInterpolated = InterpolationFunctions.interpolate(massAxis, baselinePoints, baselineNoise)
baselineInterpolated = InterpolationFunctions.interpolate(massAxis, baselinePoints, baselineValues)
baselineCorrectedAvgSpec = avgSpectrum[:,1] - baselineInterpolated;
peakIndices = PeakshapeFunctions.findPeakIndices(massAxis, avgSpectrum, baselineInterpolated, baselineNoiseInterpolated, noiseThreshold = peakfindingNoiseThresholdValue, signalLimit = peakfindingSignalLimit)
peakShapesCenterMass, peakShapesY = PeakshapeFunctions.calculatePeakshapes(massAxis, baselineCorrectedAvgSpec, peakIndices, nbrMassRegions = peakshapeRegions,regionStretch=peakshapeRegionStretch, peakWindowWidth = peakWindowWidth, quantileValue = peakshapeQuantileValue)
PyPlot.figure()
PyPlot.title("Baseline Correction")
PyPlot.semilogy(baselinePoints,baselineValues,".-")
PyPlot.semilogy(baselinePoints,baselineValues + baselineNoise,".")
#semilogy(peakMasses, peakValues, "x")
PyPlot.semilogy(massAxis,avgSpectrum)
############ delete h5 data that will be overwritten ###########
fh = HDF5.h5open(file,"r+")
if HDF5.exists(fh, "AvgBaseline")
HDF5.o_delete(fh,"AvgBaseline")
end
if HDF5.exists(fh, "MassDepPeakshape")
HDF5.o_delete(fh,"MassDepPeakshape")
end
if HDF5.exists(fh, "MassDepPeakshapeCenterMasses")
HDF5.o_delete(fh,"MassDepPeakshapeCenterMasses")
end
if HDF5.exists(fh, "BaseLines")
HDF5.o_delete(fh,"BaseLines")
end
HDF5.close(fh)
HDF5.h5write(file, "AvgBaseline", baselineInterpolated)
HDF5.h5write(file, "MassDepPeakshape", peakShapesY)
HDF5.h5write(file, "MassDepPeakshapeCenterMasses", peakShapesCenterMass)
fh = HDF5.h5open(file,"r")
ds = fh["SumSpecs"]
spectra = size(ds)[2]
HDF5.close(fh)
if (spectra > 0)
SumSpecs = HDF5.h5read(file, "SumSpecs")
fill!(SumSpecs,0)
HDF5.h5write(file, "BaseLines", SumSpecs)
else
HDF5.h5write(file, "SumSpecs", avgSpectrum)
fill!(avgSpectrum,0)
HDF5.h5write(file, "BaseLines", avgSpectrum)
end
println("DONE")
end