Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[BUG FIX] - Can not import basis set form LCModel correctly - shaokun… #597

Merged
merged 1 commit into from
Mar 30, 2023
Merged

[BUG FIX] - Can not import basis set form LCModel correctly - shaokun… #597

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 11 additions & 8 deletions libraries/FID-A/inputOutput/io_LCMBasis.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,20 +34,23 @@
% OUTPUTS:
% BASIS = Simulated basis set in FID-A structure format.

function [BASIS] = io_LCMBasis(LCMBasisSet, addMMFlag, sequence, editTarget)
function [BASIS] = io_LCMBasis(LCMBasisSet, addMMFlag, sequence, editTarget, conjugate)
%Get Osprey folder
[settingsFolder,~,~] = fileparts(which('OspreySettings.m'));
allFolders = strsplit(settingsFolder, filesep);
ospFolder = strjoin(allFolders(1:end-1), filesep); % parent folder (= Osprey folder)


% Parse input arguments
if nargin < 4
editTarget = 'none';
if nargin < 3
addMMFlag = 1;
if nargin < 2
sequence = 'unedited';
if nargin < 5
conjugate = 1;
if nargin < 4
editTarget = 'none';
if nargin < 3
addMMFlag = 1;
if nargin < 2
sequence = 'unedited';
end
end
end
end
Expand All @@ -57,7 +60,7 @@


% Load LCMBasisSet
[Read]=io_readlcmraw_basis(LCMBasisSet);
[Read]=io_readlcmraw_basis(LCMBasisSet,conjugate);
disp('We are going to rename your basis set accroding to the Osprey naming convention now.\n')
Osp_names = ['Ala ','Asc ','Asp ', 'bHB ','bHG ', 'Cit ', 'Cr ', 'CrCH2 ', 'EtOH ', 'GABA \n',...
'GPC ', 'GSH ', 'Glc ', 'Gln ', 'Glu ', 'Gly ', 'H2O ', 'mI ', 'Lac ', 'NAA ', 'NAAG \n', ...
Expand Down
196 changes: 153 additions & 43 deletions libraries/FID-A/inputOutput/io_readlcmraw_basis.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,96 +9,119 @@
%
% INPUTS:
% filename = filename of LCModel .basis file.
% conjugate = apply complex conjugate
%
% OUTPUTS:
% out = Input basis set saved as a structure in which each field is
% an individual metabolite basis spectrum in FID-A structure
% format.

function [out]=io_readlcmraw_basis(filename)
function [out]=io_readlcmraw_basis(filename,conjugate)

if nargin < 2
conjugate = 1;
end

% Begin to read the data.
fid=fopen(filename);
line=fgets(fid);

%look for FWHMBA
fwhmba_index=findstr(line,'FWHMBA =');
while isempty(fwhmba_index);
fwhmba_index=contains(line,'FWHMBA');
while ~fwhmba_index;
line=fgets(fid);
fwhmba_index=findstr(line,'FWHMBA =');
fwhmba_index=contains(line,'FWHMBA');
line = GetNumFromString(line);
end
linewidth=str2num(line(fwhmba_index+8:end-2));
linewidth=str2num(line);

%look for HZPPM
hzpppm_index=findstr(line,'HZPPPM =');
while isempty(hzpppm_index);
hzpppm_index=contains(line,'HZPPPM');
while ~hzpppm_index;
line=fgets(fid);
hzpppm_index=findstr(line,'HZPPPM =');
hzpppm_index=contains(line,'HZPPPM');
line = GetNumFromString(line);
end
hzpppm=str2num(line(hzpppm_index+8:end-2));
hzpppm=str2num(line);
Bo=hzpppm/42.577;
linewidth=linewidth*hzpppm;

%look for TE
te_index=findstr(line,'ECHOT =');
while isempty(te_index);
te_index=contains(line,'ECHOT');
while ~te_index;
line=fgets(fid);
te_index=findstr(line,'ECHOT =');
te_index=contains(line,'ECHOT');
line = GetNumFromString(line);
end
te=str2num(line(te_index+8:end-2));
te=str2num(line);

%look for spectral width
badelt_index=findstr(line,'BADELT =');
while isempty(badelt_index);
badelt_index=contains(line,'BADELT');
while ~badelt_index;
line=fgets(fid);
badelt_index=findstr(line,'BADELT =');
badelt_index=contains(line,'BADELT');
line = GetNumFromString(line);
end
dwelltime=str2num(line(badelt_index+8:end-2));
dwelltime=str2num(line);
spectralwidth=1/dwelltime;
fileEnd=false;

while ~feof(fid)
%look for a center frequency
ppmsep_index=findstr(line,'PPMSEP =');
while isempty(ppmsep_index);
ppmsep_index=contains(line,'PPMSEP');
while ~ppmsep_index;
line=fgets(fid);
ppmsep_index=findstr(line,'PPMSEP =');

if findstr(line,'METABO =')
break
end

if ischar(line)
ppmsep_index=contains(line,'PPMSEP');
if contains(line,'METABO')
break
end
end
end
if ~isempty(ppmsep_index)
centerFreq=str2num(line(ppmsep_index+8:end-2));
if ppmsep_index
line = GetNumFromString(line);
centerFreq=str2num(line);
else
centerFreq = [];
end

%Look for the metabolite name
metab_index=findstr(line,'METABO =');
while isempty(metab_index);
metab_index=contains(line,'METABO');
while ~metab_index;
line=fgets(fid);
metab_index=findstr(line,'METABO =');
if ischar(line)
metab_index=contains(line,'METABO');
if metab_index
break
end
end
end
if strcmp(line(end-3), '''')
metabName=line(metab_index+10:end-4);
metabName=line(10:end-4);
else
metabName=line(metab_index+10:end-3);
metabName=line(10:end-3);
end


hdrEnd_index=findstr(line,'$END');
while isempty(hdrEnd_index);
hdrEnd_index=contains(line,'$END');
while ~hdrEnd_index;
line=fgets(fid);
hdrEnd_index=findstr(line,'$END');
if ischar(line)
hdrEnd_index=contains(line,'$END');
if hdrEnd_index
break
end
end
end

line=fgets(fid);

nmused_index=findstr(line,'$NMUSED');
nmused_index=contains(line,'$NMUSED');
basis_index=contains(line,'$BASIS');
linenum=1;
RF=[];
% If the line is empty skip it
while ~isempty(line) && isempty(nmused_index) && ~fileEnd
while ~isempty(line) && ~nmused_index && ~basis_index && ~fileEnd
%dataline=line(1:semicol_index-2);
[A,count, errmsg, nextindex] = sscanf(line, '%f', inf);
% If read failed, output the error
Expand All @@ -113,10 +136,29 @@
end
linenum = linenum + 1;
line=fgets(fid);
nmused_index=findstr(line,'$NMUSED');
if ischar(line)
nmused_index=contains(line,'$NMUSED');
basis_index=contains(line,'$BASIS');
end
end
specs=RF(1:2:end) + 1i*RF(2:2:end);
specs=flipud(fftshift(conj(specs),1));

% GO 2022/01/24
% LCModel uses negative BADELT values to encrypt basis sets
% (LCModel.f source code lines 3666 and following)
if dwelltime < 0
dix = 1499;
for rr = 1:length(specs)
[randomresult, dix] = lcmodelrng(dix);
specs(rr) = -specs(rr) .* exp(-20*randomresult + 10);
end
end

if conjugate
specs=flipud(fftshift(conj(specs),1));
else
specs=(fftshift(specs,1));
end
vectorsize=length(specs);
sz=[vectorsize 1];
if mod(vectorsize,2)==0
Expand All @@ -131,24 +173,28 @@
ppm=ppm+4.68;
t=[dwelltime:dwelltime:vectorsize*dwelltime];
txfrq=hzpppm*1e6;
metabName = RemoveWhiteSpaces(metabName);
if strcmp(metabName,'-CrCH2')
metabName = 'CrCH2';
end
if strcmp(metabName,'2HG')
metabName = 'bHG';
end
eval(['out.' metabName '.fids=fids;']);
eval(['out.' metabName '.specs=specs;']);
eval(['out.' metabName '.sz=[vectorsize 1 1 1];']);
eval(['out.' metabName '.n=vectorsize;']);
eval(['out.' metabName '.spectralwidth=spectralwidth;']);
eval(['out.' metabName '.spectralwidth=abs(spectralwidth);']);
eval(['out.' metabName '.sz=sz;']);
eval(['out.' metabName '.Bo=Bo;']);
eval(['out.' metabName '.te=te;']);
eval(['out.' metabName '.tr=[];']);
eval(['out.' metabName '.dwelltime=1/spectralwidth;']);
eval(['out.' metabName '.dwelltime=abs(1/spectralwidth);']);
eval(['out.' metabName '.linewidth=linewidth;']);
eval(['out.' metabName '.ppm=ppm;']);
eval(['out.' metabName '.t=t;']);
eval(['out.' metabName '.txfrq=txfrq;']);
if ~isempty(ppmsep_index)
if ~isempty(centerFreq)
eval(['out.' metabName '.centerFreq=centerFreq;']);
end
eval(['out.' metabName '.date=date;']);
Expand Down Expand Up @@ -177,8 +223,72 @@
end

fclose(fid);
end

% RF=RF';
% rf(:,1)=RF(:,2)*180/pi;
% rf(:,2)=RF(:,1);
% rf(:,3)=ones(length(RF(:,1)),1);
% rf(:,3)=ones(length(RF(:,1)),1);

function str3 = GetNumFromString(str)
str1 = regexprep(str,'[,;=]', ' ');
str2 = regexprep(regexprep(str1,'[^- 0-9.eE(,)/]',''), ' \D* ',' ');
str3 = regexprep(str2, {'\.\s','\E\s','\e\s','\s\E','\s\e'},' ');
end

function str = RemoveWhiteSpaces(str)
pattern = '[ \t\n]'; % Match zero or more spaces, tabs, or newlines, followed by a double quote
replacement = ''; % Replace the matched string with just a double quote
str = regexprep(str, pattern, replacement);
end


% GO 2022/01/24
% LCModel uses negative BADELT values to encrypt basis sets
% (LCModel.f source code lines 3666 and following)
%++++++++++++++++ doubleprecision VERSION 2DP (MAR 1984) ++++++++++++++ 4222
% FUNCTION RANDOM. PRODUCES A PSEUDORANDOM REAL ON THE OPEN INTERVAL 4223
% (0.,1.). 4224
% DIX (IN doubleprecision) MUST BE INITIALIZED TO A WHOLE NUMBER 4225
% BETWEEN 1.0D0 AND 2147483646.0D0 BEFORE THE FIRST CALL TO RANDOM 4226
% AND NOT CHANGED BETWEEN SUCCESSIVE CALLS TO RANDOM. 4227
% BASED ON L. SCHRAGE, ACM TRANS. ON MATH. SOFTWARE 5, 132 (1979). 4228
%----------------------------------------------------------------------- 4229
function [randomresult,dix]=lcmodelrng(dix);

a = 16807.0d0;
b15 = 32768.d0;
b16 = 65536.0d0;
p = 2147483647.0d0;

% 4231
% PORTABLE RANDOM NUMBER GENERATOR 4232
% USING THE RECURSION 4233
% DIX = DIX*A MOD P 4234
% 4235
% 4237
% 7**5, 2**15, 2**16, 2**31-1 4238

% GET 15 HI ORDER BITS OF DIX 4241
xhi = dix./b16;
xhi = xhi - rem(xhi,1.0d0);
% GET 16 LO BITS IF DIX AND FORM LO PRODUCT 4244
xalo =(dix-xhi.*b16).*a;
% GET 15 HI ORDER BITS OF LO PRODUCT 4246
leftlo = xalo./b16;
leftlo = leftlo - rem(leftlo,1.0d0);
% FORM THE 31 HIGHEST BITS OF FULL PRODUCT 4249
fhi = xhi.*a + leftlo;
% GET OVERFLO PAST 31ST BIT OF FULL PRODUCT 4251
k = fix(fhi./b15);
k = fix(k - rem(k,1.0d0));
% ASSEMBLE ALL THE PARTS AND PRESUBTRACT P 4254
% THE PARENTHESES ARE ESSENTIAL 4255
dix =(((xalo-leftlo.*b16)-p)+(fhi-k.*b15).*b16) + k;
% ADD P BACK IN IF NECESSARY 4257
if(dix < 0.0d0)
dix = dix + p;
end
% MULTIPLY BY 1/(2**31-1) 4259
randomresult = dix.*4.656612875d-10;
end %function random