Skip to content

Commit 9f11e1f

Browse files
YohkoYohko
committed
load and process Agilent GCMS data
1 parent fcff63b commit 9f11e1f

File tree

5 files changed

+439
-0
lines changed

5 files changed

+439
-0
lines changed

‎gc data analysis/GC_AgilentloadDATAMS.m‎

Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
%Licence: GNU General Public License version 2 (GPLv2)
2+
function data = GC_AgilentloadDATAMS(fid)
3+
data = [];
4+
fseek(fid,4,'bof');
5+
strlen = fread(fid,1,'uint8');
6+
datatype = convertCharsToStrings(char(fread(fid,strlen,'char*1')));
7+
% sample description
8+
fseek(fid,24,'bof');
9+
strlen = fread(fid,1,'uint8');
10+
filename = convertCharsToStrings(char(fread(fid,strlen,'char*1')));
11+
12+
fseek(fid,178,'bof');
13+
strlen = fread(fid,1,'uint8');
14+
timestr = convertCharsToStrings(char(fread(fid,strlen,'char*1')));
15+
16+
fseek(fid,208,'bof');
17+
strlen = fread(fid,1,'uint8');
18+
convertCharsToStrings(char(fread(fid,strlen,'char*1'))); % GCMS
19+
20+
fseek(fid,228,'bof');
21+
strlen = fread(fid,1,'uint8');
22+
methodstr = convertCharsToStrings(char(fread(fid,strlen,'char*1')));
23+
24+
fseek(fid,448,'bof'); % GCMS
25+
strrep(convertCharsToStrings(char(fread(fid,168/2,'char*1',1))),char(00),'');
26+
% method directory
27+
fseek(fid,616,'bof');
28+
strrep(convertCharsToStrings(char(fread(fid,510/2,'char*1',1))),char(00),'');
29+
% name of method
30+
fseek(fid,1126,'bof'); % 616 + 510
31+
strrep(convertCharsToStrings(char(fread(fid,510/2,'char*1',1))),char(00),'');
32+
% tune file directory
33+
fseek(fid,1636,'bof'); % 616 + 510 + 510
34+
strrep(convertCharsToStrings(char(fread(fid,510/2,'char*1',1))),char(00),'');
35+
% MSD tune file
36+
fseek(fid,2146,'bof'); % 616 + 510 + 510 + 510
37+
strrep(convertCharsToStrings(char(fread(fid,510/2,'char*1',1))),char(00),'');
38+
% sample description
39+
fseek(fid,4236,'bof'); % 616 + 510 + 510 + 510
40+
strrep(convertCharsToStrings(char(fread(fid,510/2,'char*1',1))),char(00),'');
41+
42+
% get offset to data
43+
fseek(fid,266,'bof');
44+
dataoffset = 2*fread(fid,1,'uint16','ieee-be')-2;
45+
% get number of scans
46+
switch datatype
47+
case 'GC / MS Data File'
48+
fseek(fid,322,'bof');
49+
otherwise
50+
fseek(fid,280,'bof');
51+
end
52+
nscans = fread(fid,1,'uint16','ieee-be');
53+
54+
% MSD data
55+
fseek(fid,dataoffset,'bof');
56+
points = 0;
57+
pointvector = zeros(nscans+1,1); % points for each scan (for each time point)
58+
TIC = zeros(nscans,2);
59+
datams = zeros(0,3);
60+
fseek(fid,dataoffset,'bof');
61+
for ii = 1:nscans
62+
% get the position of the next scan
63+
npos = ftell(fid) + 2 * fread(fid,1,'uint16','ieee-be');
64+
% keep a running total of how many measurements
65+
points = points + (npos - ftell(fid) - 26) / 4;
66+
pointvector(ii+1) = points;
67+
% get time, 60000 --> ms to min
68+
TIC(ii,1) = fread(fid,1,'uint32','ieee-be') / 60000;
69+
% what are these 12 bytes??
70+
%disp(ftell(fid))
71+
fseek(fid,ftell(fid)+12,'bof');
72+
npts = pointvector(ii+1)-pointvector(ii);
73+
lastscan = ii+1;
74+
if(npts<0)
75+
lastscan = ii;
76+
break
77+
else
78+
% ions and counts
79+
datapoints = fread(fid,2*npts,'uint16','ieee-be');
80+
datams((pointvector(ii)+1):(pointvector(ii+1)),2) = datapoints(1:2:(end)); % m/z
81+
datams((pointvector(ii)+1):pointvector(ii+1),3) = datapoints(2:2:end); % counts
82+
datams((pointvector(ii)+1):pointvector(ii+1),1) = TIC(ii,1);
83+
% what are the next 6bytes?
84+
%disp(ftell(fid))
85+
%disp('6bytes at end')
86+
%fread(fid,3,'uint16','ieee-be')
87+
fseek(fid,npos-4,'bof');
88+
TIC(ii,2) = fread(fid,1,'uint32','ieee-be');
89+
end
90+
fseek(fid,npos,'bof'); % goto next scan
91+
end
92+
93+
datams(:,2) = datams(:,2)/20;
94+
datams(:,3) = bitand(datams(:,3), 16383,'uint16').*8.^ bitshift(datams(:,3),-14,'uint16');
95+
TIC = TIC(1:(lastscan-1),:);
96+
%% return data
97+
timecodes = textscan(timestr,'%s');
98+
day = str2double(timecodes{1}{1});
99+
year = 2000+str2double(timecodes{1}{3});
100+
tmpstr = timecodes{1}{4};
101+
hr = str2double(tmpstr(1:2));
102+
minute = str2double(tmpstr(4:5));
103+
seconds = 0;
104+
month = 0;
105+
switch timecodes{1}{5}
106+
case 'pm'
107+
hr = hr+12;
108+
case 'am'
109+
if(hr == 12)
110+
hr = 0;
111+
end
112+
end
113+
switch timecodes{1}{2}
114+
case 'Jan'
115+
month = 1;
116+
case 'Feb'
117+
month = 2;
118+
case 'Mar'
119+
month = 3;
120+
case 'Apr'
121+
month = 4;
122+
case 'May'
123+
month = 5;
124+
case 'Jun'
125+
month = 6;
126+
case 'Jul'
127+
month = 7;
128+
case 'Aug'
129+
month = 8;
130+
case 'Sep'
131+
month = 9;
132+
case 'Oct'
133+
month = 10;
134+
case 'Nov'
135+
month = 11;
136+
case 'Dec'
137+
month = 12;
138+
otherwise
139+
disp('Unknown Month');
140+
disp(timecodes{1}{3});
141+
end
142+
timecode = posixtime(datetime(year,month,day,hr,minute,seconds));
143+
data = {timecode; datams; TIC};
144+
end

‎gc data analysis/GC_AgilentloadTCDbin.m‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
%Licence: GNU General Public License version 2 (GPLv2)
2+
function data = GC_AgilentloadTCDbin(fid)
3+
fseek(fid,282,'bof');
4+
starttime = fread(fid,1,'single','ieee-be')/ 60000;
5+
endtime = fread(fid,1,'single','ieee-be')/ 60000;
6+
%
7+
fseek(fid,326,'bof');
8+
strlen = fread(fid,1,'uint8');
9+
version = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
10+
11+
% filetype
12+
fseek(fid,347,'bof');
13+
strlen = fread(fid,1,'uint8');
14+
filetype = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
15+
16+
% filename
17+
fseek(fid,858,'bof');
18+
strlen = fread(fid,1,'uint8');
19+
filename = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
20+
21+
% date and time
22+
fseek(fid,2391,'bof');
23+
strlen = fread(fid,1,'uint8');
24+
timecodes = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
25+
26+
%model number
27+
fseek(fid,2492,'bof');
28+
strlen = fread(fid,1,'uint8');
29+
instrumentmodel = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
30+
31+
%instrument
32+
fseek(fid,2533,'bof');
33+
strlen = fread(fid,1,'uint8');
34+
instrumentname = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
35+
36+
%method file
37+
fseek(fid,2574,'bof');
38+
strlen = fread(fid,1,'uint8');
39+
methodname = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
40+
41+
fseek(fid,3089,'bof');
42+
strlen = fread(fid,1,'uint8');
43+
programname = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
44+
fseek(fid,4125,'bof');
45+
datarate = fread(fid,1,'uint8'); % datarate?, yes in Hz
46+
%y-unit??
47+
fseek(fid,4172,'bof');
48+
strlen = fread(fid,1,'uint8');
49+
instrumentmodel = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
50+
51+
fseek(fid,4213,'bof');
52+
strlen = fread(fid,1,'uint8');
53+
signalname = convertCharsToStrings(char(fread(fid,strlen,'char*1',1)));
54+
% data
55+
fseek(fid,6144,'bof'); % jump to beginning of data
56+
yaxis = fread(fid,'double');
57+
ftell(fid);
58+
xaxis = linspace(starttime,endtime,length(yaxis));
59+
y(:,1) = xaxis;
60+
y(:,2) = yaxis;
61+
% data
62+
% fseek(fid,6152,'bof'); % jump to beginning of data
63+
% yaxis = fread(fid,'double');
64+
% ftell(fid);
65+
% xaxis = linspace(0,1/datarate/60*(length(yaxis)-1),length(yaxis));
66+
% y(:,1) = xaxis;
67+
% y(:,2) = yaxis;
68+
%
69+
timecodes = textscan(timecodes,'%s');
70+
day = str2double(timecodes{1}{1});
71+
year = 2000+str2double(timecodes{1}{3});
72+
tmpstr = timecodes{1}{4};
73+
hr = str2double(tmpstr(1:2));
74+
minute = str2double(tmpstr(4:5));
75+
seconds = 0;
76+
month = 0;
77+
switch timecodes{1}{5}
78+
case 'pm'
79+
hr = hr+12;
80+
case 'am'
81+
if(hr == 12)
82+
hr = 0;
83+
end
84+
end
85+
switch timecodes{1}{2}
86+
case 'Jan'
87+
month = 1;
88+
case 'Feb'
89+
month = 2;
90+
case 'Mar'
91+
month = 3;
92+
case 'Apr'
93+
month = 4;
94+
case 'May'
95+
month = 5;
96+
case 'Jun'
97+
month = 6;
98+
case 'Jul'
99+
month = 7;
100+
case 'Aug'
101+
month = 8;
102+
case 'Sep'
103+
month = 9;
104+
case 'Oct'
105+
month = 10;
106+
case 'Nov'
107+
month = 11;
108+
case 'Dec'
109+
month = 12;
110+
otherwise
111+
disp('Unknown Month');
112+
disp(timecodes{1}{3});
113+
end
114+
timecode = posixtime(datetime(year,month,day,hr,minute,seconds));
115+
data = {timecode; y};
116+
end

‎gc data analysis/GC_AgilentloadTICASCII.m‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
%Licence: GNU General Public License version 2 (GPLv2)
2+
function data = GC_AgilentloadTICASCII(fid)
3+
timecode = 0;
4+
h = fgets(fid);
5+
timecodes = textscan(h,'%s');
6+
day = str2double(timecodes{1}{2});
7+
year = 2000+str2double(timecodes{1}{4});
8+
tmpstr = timecodes{1}{5};
9+
hr = str2double(tmpstr(1:2));
10+
minute = str2double(tmpstr(4:5));
11+
seconds = 0;
12+
month = 0;
13+
switch timecodes{1}{6}
14+
case 'pm'
15+
hr = hr+12;
16+
case 'am'
17+
if(hr == 12)
18+
hr = 0;
19+
end
20+
end
21+
switch timecodes{1}{3}
22+
case 'Jan'
23+
month = 1;
24+
case 'Feb'
25+
month = 2;
26+
case 'Mar'
27+
month = 3;
28+
case 'Apr'
29+
month = 4;
30+
case 'May'
31+
month = 5;
32+
case 'Jun'
33+
month = 6;
34+
case 'Jul'
35+
month = 7;
36+
case 'Aug'
37+
month = 8;
38+
case 'Sep'
39+
month = 9;
40+
case 'Oct'
41+
month = 10;
42+
case 'Nov'
43+
month = 11;
44+
case 'Dec'
45+
month = 12;
46+
otherwise
47+
disp('Unknown Month');
48+
disp(timecodes{1}{3});
49+
end
50+
timecode = posixtime(datetime(year,month,day,hr,minute,seconds));
51+
fgets(fid); % 'Start of data points'
52+
y = textscan(fid,'%f','delimiter',',','emptyvalue', NaN);
53+
rows = size(y{1},1)/2;
54+
y =[reshape(y{1},[2,rows])]';
55+
data = {timecode; y};
56+
end

‎gc data analysis/GC_MSDanalysiswidget.m‎

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
%Licence: GNU General Public License version 2 (GPLv2)
2+
function GC_MSDanalysiswidget()
3+
data = GC_dloadAgilent;
4+
if(~isempty(data))
5+
S.f = figure;
6+
7+
S.dataTIC = data(1).spectrum;
8+
S.dataTCD = data(2).spectrum;
9+
S.dataMSD = data(3).spectrum;
10+
S.dataMSDTIC = data(4).spectrum;
11+
idx = find(~isnan(S.dataMSD(:,3)) & ~isinf(S.dataMSD(:,3)));
12+
S.F = scatteredInterpolant(S.dataMSD(idx,1), S.dataMSD(idx,2),S.dataMSD(idx,3));
13+
min_x = round(min(S.dataMSD(idx,1))*10)/10;
14+
max_x = round(max(S.dataMSD(idx,1))*10)/10;
15+
min_y = round(min(S.dataMSD(idx,2))*1)/1;
16+
max_y = round(max(S.dataMSD(idx,2))*1)/1;
17+
S.xvec = min_x:0.01:max_x; % 0.01 min
18+
S.yvec = min_y:1:max_y; % only full m/z
19+
[xq, yq] = meshgrid(S.xvec,S.yvec);
20+
S.Atest = S.F(xq, yq);
21+
%S.bg = sum(S.Atest(:,1:5),2)/5;
22+
23+
GC_settings_graph;
24+
25+
S.ax1 = subplot(1,2,1);
26+
S.h1=plot(S.dataMSDTIC(:,1),S.dataMSDTIC(:,2),'linewidth', f_line);
27+
ylabel('counts', 'fontsize',f_caption);
28+
xlabel('Retention Time / min', 'fontsize',f_caption);
29+
set(gca, 'linewidth', f_line);
30+
set(gca, 'fontsize', f_caption);
31+
title('MSD TIC', 'fontsize',10);
32+
set(gca, 'units','normalized', 'Position', [0.1, 0.2, 0.35, 0.7]);
33+
34+
S.ax2 = subplot(1,2,2);
35+
imagesc(S.xvec,S.yvec,S.Atest);
36+
ylabel('m/z', 'fontsize',f_caption);
37+
xlabel('Retention Time / min', 'fontsize',f_caption);
38+
set(gca, 'linewidth', f_line);
39+
set(gca, 'fontsize', f_caption);
40+
title('MSD', 'fontsize',10);
41+
xlim([min(S.dataMSD(:,1)) max(S.dataMSD(:,1))]);
42+
ylim([min(S.dataMSD(:,2)) max(S.dataMSD(:,2))]);
43+
set(gca, 'units','normalized', 'Position', [0.6, 0.2, 0.35, 0.7]);
44+
45+
S.edit_mz = uicontrol('Style','edit',...
46+
'units','pix',...
47+
'position',[10 0 50 20],...
48+
'fontsize',14,...
49+
'string','0',...
50+
'callback',{@edit_mz_call,S});
51+
end
52+
53+
54+
function edit_mz_call(varargin)
55+
mzedit = varargin{1};
56+
S = varargin{3};
57+
axes(S.ax1);
58+
tmpval = str2double(mzedit.String);
59+
try
60+
idx = find(S.yvec == tmpval);
61+
catch
62+
idx = 1;
63+
end
64+
if(isempty(idx))
65+
idx = 1;
66+
mzedit.String = 'error';
67+
end
68+
GC_settings_graph;
69+
plot(S.xvec,S.Atest(idx,:),'linewidth', f_line);
70+
legend(sprintf('m/z = %s',mzedit.String));
71+
ylabel('counts', 'fontsize',f_caption);
72+
xlabel('Retention Time / min', 'fontsize',f_caption);
73+
set(gca, 'linewidth', f_line);
74+
set(gca, 'fontsize', f_caption);
75+
%title('MSD TIC', 'fontsize',10);
76+

0 commit comments

Comments
 (0)