function pik = readlcmstxtfile(ithreshperc, chrglw, chrghg, deamidatn, oxidatn)

% function 'readlcmstxtfile': this function reads text files from experimental 
% LC-MS1 data and converts them into a [4,npeaks] peak table required by MicrobeMS
% parseuniprot - create in silico databases for mass spectrometry-
% based microbial identification, last revision: January 17, 2020
% Copyright Peter Lasch, RKI/ZBS6 Berlin, 2017 - 2020
%
% Usage:
%
% >> pik = readlcmstxtfile(ithreshperc, chrglw, chrghg, deamidatn, oxidatn);
%
% Input argument: ithreshperc  - a scalar between 0 and 150 with an impact on the
%                                number of MS1 peaks to be removed
%                                the higher ithreshperc, the lower the resulting 
%                                number of peaks
%                                
%                 oxidatn      - logical value, indicates to perform / not perform 
%                                the test for oxidation peaks (default: 1)
%                 deamidatn    - logical value, indicates to perform / not perform 
%                                the test for deamidation peaks (default: 1)
%                 chrglw       - minimum allowed charge state (2 - default)
%                 chrghg       - maximum allowed charge state (4 - default)

%
% Output argument: pik - a matrix of size [4, npeaks]
%                                npeaks - the number of MS1 peaks 
%                                pik(1,1:npeaks) - the vector of peptide mass values (M)
%                                pik(2,1:npeaks) - the respective intensity values
%                                pik(3,1:npeaks) - not used but required by MicrobeMS
%                                pik(4,1:npeaks) - not used but required by MicrobeMS

%% Definitions 
mwproton  = 1.007276466; % proton mass 
mwoxid    = 15.99491;    % monoisotopic mass shift of oxidized peptides
mwdeamid  = 0.98402;     % monoisotopic mass shift due to deamidation of Asn and Gln 
tol       = 0.00004;     % allowed mass tolerance (in mass units)

% Basic checks
if nargin == 0
    ithreshperc = 50;
    deamidatn = 1;       % 1- perform / 0 - do not perform check for deamidation peaks
    oxidatn   = 1;       % 1- perform / 0 - do not perform check for oxidation peaks    
    chrglw    = 2;       % minimum allowed charge state (2 - default)
    chrghg    = 4;       % maximum allowed charge state (4 - default)
elseif nargin == 1
    deamidatn = 1;
    oxidatn   = 1;
    chrglw    = 2;
    chrghg    = 4;
elseif nargin == 3
    deamidatn = 1;
    oxidatn   = 1;
elseif nargin ~= 5
    error('inadequate number of arguments for function ''readlcmstxtfile''');
end

if isnumeric(ithreshperc) == 0 || isinf(ithreshperc) == 1 || isnan(ithreshperc) == 1 || ithreshperc > 150 || ithreshperc <= 0
    warning('invalid input for ''ithreshperc'' - the default value of ithreshperc is used (50)')
    ithreshperc = 50;
end
if isnumeric(chrglw) == 0 || isinf(chrglw) == 1 || isnan(chrglw) == 1 || chrglw > 7 || chrglw < 0
    warning('invalid input for ''chrglw'' - the default value of ''chrglw'' is used (2)')
    chrglw = 2;
end
if isnumeric(chrghg) == 0 || isinf(chrghg) == 1 || isnan(chrghg) == 1 || chrghg > 10 || chrghg <= 1
    warning('invalid input for ''chrghg'' - the default value of ''chrghg'' is used (4)')
    chrghg = 4;
end

if chrglw >= chrghg
    error('''chrglw'' must not equal or larger than ''chrghg''');
end

 
% PART 1. open the text file
[fname, pname] = uigetfile('*.txt','Select a LC-MS1 text file');
if fname == 0
   pik = [];
   disp('file not loaded')
   return
else
    fullfilename = [pname fname];
end
cd(pname);
disp(['loading file: ' fullfilename])

lcmsfid         = fopen(fullfilename);
if lcmsfid == -1
   error('error opening LC-MS1 file')
end

% PART 2. display settings
disp('pik = readlcmstxtfile(ithreshperc, deamidatn, oxidatn)')
disp(['   ithreshperc : ' num2str(ithreshperc)]);
disp(['   deamidatn   : ' num2str(deamidatn)]);
disp(['   oxidatn     : ' num2str(oxidatn)]);
disp(['   chrglw      : ' num2str(chrglw)]);
disp(['   chrghg      : ' num2str(chrghg)]);
disp(['   pname       : ' pname]);
disp(['   fname       : ' fname]);

% PART 3. load the data
% get the number of lines in the file
% skip first line 
tmp = fgetl(lcmsfid);
k = 0;
while ~feof(lcmsfid)
   tmp = fgetl(lcmsfid);
   k=k+1;
end
frewind(lcmsfid);

% define vectors:
mdivz = zeros(1,k);
z     = zeros(1,k);
intns = zeros(1,k);

% read the data from file lcmsfid
% skip first line 
tmp = fgetl(lcmsfid);
t   = 0;
for i=1:k
   txtln        = fgetl(lcmsfid);
   txtln        = strrep(txtln,'"','');
   txtln        = strrep(txtln,'False','');
   flt32        = str2num(txtln);
   chrg         = flt32(2);   % charge state
   if chrg >= chrglw  && chrg <= chrghg
       t            = t+1;   
       mdivz(t)     = flt32(1);
       z(t)         = chrg;
       intns(t)     = flt32(10);         
   end
end
fclose(lcmsfid);

mdivz          = mdivz(1:t);
z              = z(1:t);
intns          = intns(1:t);         

% PART 4. create the peak table
% in this approach subtract first the proton mass from the m/z value
% of the first column and multiply then the difference mdivz by the charge
mdivz     = mdivz - mwproton;
mw        = mdivz.*z;
% resort vectors
[mw, idx] = sort(mw);
intns     = intns(idx);
t         = size(mw,2);
disp(['   total number of peaks in spectrum            : ' num2str(t)])

% PART 5. prepare ui elements
scrsz  = get(0,'ScreenSize');
 posw  = [(scrsz(3)-1000)/2 (scrsz(4)-800)/2 1000 800]; % size of ther figure
posax1 = [ 70 450 400 300];
posax2 = [570 450 400 300];
posax3 = [ 70  70 400 300];
posax4 = [570  70 400 300];

h.fig   = figure('NumberTitle','off','Units','Pixel','Name',['LC-MS spectrum ' fname],'Position',posw,'MenuBar','figure','Resize','off','PaperPositionMode','auto');
h.ax01  = axes('Parent',h.fig,'Units','Pixel','Box','on','CameraUpVector',[0 1 0],'CameraUpVectorMode','manual', ...
               'Position',	posax1, 'Interruptible','off','SelectionHighlight','off','XMinorTick','on','TickLength',[0.003 0], ...
               'HandleVisibility','on','NextPlot','ReplaceChildren','XGrid','off','YGrid','off','ZGrid','off');
set(get(h.ax01,'Title'),'String','Original data: Number of peaks with log10 intensity')
set(get(h.ax01,'XLabel'),'String','log(10) intensity')
set(get(h.ax01,'YLabel'),'String','# of peaks')
h.ax02  = axes('Parent',h.fig,'Units','Pixel','Box','on','CameraUpVector',[0 1 0],'CameraUpVectorMode','manual', ...
               'Position',	posax2, 'Interruptible','off','SelectionHighlight','off','XMinorTick','on','TickLength',[0.003 0], ...
               'HandleVisibility','on','NextPlot','ReplaceChildren','XGrid','off','YGrid','off','XLim',[500 7000]);
set(get(h.ax02,'Title'),'String','Original data: Number of peaks per mass region')
set(get(h.ax02,'XLabel'),'String','mass')
set(get(h.ax02,'YLabel'),'String','# of peaks')
h.ax03  = axes('Parent',h.fig,'Units','Pixel','Box','on','CameraUpVector',[0 1 0],'CameraUpVectorMode','manual', ...
               'Position',	posax3, 'Interruptible','off','SelectionHighlight','off','XMinorTick','on','TickLength',[0.003 0], ...
               'HandleVisibility','on','NextPlot','ReplaceChildren','XGrid','off','YGrid','off','ZGrid','off');
set(get(h.ax03,'Title'),'String','Pre-processed data: Number of peaks with log10 intensity')
set(get(h.ax03,'XLabel'),'String','log(10) intensity')
set(get(h.ax03,'YLabel'),'String','# of peaks')
h.ax04  = axes('Parent',h.fig,'Units','Pixel','Box','on','CameraUpVector',[0 1 0],'CameraUpVectorMode','manual', ...
               'Position',	posax4, 'Interruptible','off','SelectionHighlight','off','XMinorTick','on','TickLength',[0.003 0], ...
               'HandleVisibility','on','NextPlot','ReplaceChildren','XGrid','off','YGrid','off','XLim',[500 7000]');
set(get(h.ax04,'Title'),'String','Pre-processed data: Number of peaks per mass region')
set(get(h.ax04,'XLabel'),'String','mass')
set(get(h.ax04,'YLabel'),'String','# of peaks')

% PART 6. apply an arbitrary intensity threshold
% ithreshperc: defines the number of peaks to be accepted
% sort peaks according to the intensity:
% axis 01: Number of peaks with log10 intensity (original data)
[~,ctrs1]   = hist(log10(intns),100);
ctrs1       = ctrs1';
hist(h.ax01,log10(intns),ctrs1);

% axis 02: Histogram: Number of peaks per mass region (original data)
[~,ctrs2]   = hist(mw,100);
ctrs2       = ctrs2';
hist(h.ax02,mw,ctrs2);
drawnow;

% remove low intensity peaks: low mass peaks are preferentially 
% removed, each mass bin is treated separately, see function redpeaknmb

[mw,intns] = redpeaknmb(mw,intns,ithreshperc);
% [mw,intns] = redpeaknmb_old(mw,intns,ithreshperc); % outcommended older version of redpeaknmb
t          = size(mw,2);
disp(['   # peaks after applying an intensity threshold: ' num2str(t)])

% PART 7. find & remove peaks due to oxidation
if isnumeric(oxidatn) == 1 && oxidatn == 1
    for i=1:t
        diff = abs(mw+mwoxid-mw(i));
        didx = find(diff<tol,1);
        if isempty(didx) ~= 1
            intns(i) = 0;
        end
    end
    idxnrm    = intns>0;
    intns     = intns(idxnrm);
    mw        = mw(idxnrm);
    t         = size(mw,2);
    disp(['   # peaks after removing oxidation peaks       : ' num2str(t)])
end

% PART 8. find & remove peaks due to deamidation
if isnumeric(deamidatn) == 1 && deamidatn == 1
    for i=1:t
        diff = abs(mw+mwdeamid-mw(i));
        didx = find(diff<tol,1);
        if isempty(didx) ~= 1
            intns(i) = 0;
        end
    end
    idxnrm    = intns>0;
    intns     = intns(idxnrm);
    mw        = mw(idxnrm);
    t         = size(mw,2);
    disp(['   # peaks after removing deamidation peaks     : ' num2str(t)])
end

% axis 03: Number of peaks with log10 intensity (corrected)
xlim   = get(h.ax01,'XLim');
ylim   = get(h.ax01,'YLim');
hist(h.ax03,log10(intns),ctrs1);
set(h.ax03,'XLim',xlim,'YLim',ylim);
h.ptch3 = findobj(h.ax03,'Type','patch');
set(h.ptch3,'FaceColor',[1 0 0],'EdgeColor','w')

% axis 04: Histogram: Number of peaks per mass region (corrected)
xlim   = get(h.ax02,'XLim');
ylim   = get(h.ax02,'YLim');
hist(h.ax04,mw,ctrs2);
set(h.ax04,'XLim',xlim,'YLim',ylim);
h.ptc4 = findobj(h.ax04,'Type','patch');
set(h.ptc4,'FaceColor',[1 0 0],'EdgeColor','w')

% create the array 'pik':
pik      = ones(4,t);
pik(1,:) = mw;
pik(2,:) = intns;
%%

%% ************************************************************************
function [nmw,nintns] = redpeaknmb(mw,intns,ithreshperc)
% *************************************************************************
% function to reduce peak number by removing low intensity peaks preferentially 
% from low mass bins 
% bin(i)   - bin number i
% pcrm(i)  - percentage of data to be removed in the i-th bin
ithresh    = ithreshperc/100; 
bin(1,:)   = [  0   700]; 
bin(2,:)   = [700   900]; 
bin(3,:)   = [900  1100]; 
bin(4,:)   = [1100 1300]; 
bin(5,:)   = [1300 1500]; 
bin(6,:)   = [1500 1700]; 
bin(7,:)   = [1700 1900]; 
bin(8,:)   = [1900 2100]; 
bin(9,:)   = [2100 2300]; 
bin(10,:)  = [2300 2500];
bin(11,:)  = [2500 2800]; 
bin(12,:)  = [2800 3300]; 
bin(13,:)  = [3300 3800]; 
bin(14,:)  = [3800 4400]; 
bin(15,:)  = [4400 5000]; 
bin(16,:)  = [5000 10000];

nbin       = size(bin,1);                         % number of bins (--> 16)
rmf1       = 0.50;
rmf2       = 0.02;
pcrmstp    = (rmf1-rmf2)/(nbin-1);
pcrm       = [rmf1:-pcrmstp:rmf2];
pcrm       = pcrm +(ithresh-0.5);
tidx       = pcrm>1;
pcrm(tidx) = 1;
tidx       = pcrm<0;
pcrm(tidx) = 0;
n          = zeros(nbin,1);
t          = zeros(nbin,1);
thresh     = zeros(nbin,1);

for i=1:nbin
   idx       = find(mw>=bin(i,1) & mw<bin(i,2));  % indices of signals of the i-th bin
   if numel(idx) == 0
       continue
   end
   lgthidx   = size(idx,2);
   binint    = intns(idx);                        % intensities of respective signals of the i-th bin
   binmw     = mw(idx);                           % mw values of signals of the i-th bin
   n(i)      = numel(binint);                     % total number of peaks present in the i-th bin
   srtint    = sort(binint);                      % sorted signal intensity vector of the i-th bin
   srtindidx = 1+floor(pcrm(i)*n(i));             % actual index in srtind
   if srtindidx <= lgthidx
        thresh(i) = srtint(srtindidx);            % intensity threshold in the given bin
   else
        thresh(i) = srtint(srtindidx-1);          % intensity threshold in the given bin
   end
   idxint    = binint>thresh(i);                  % indices of intensity values above the bin-specific threshold
   rbinint   = binint(idxint);                    % reduced intensity vector of bin i
   rbinmw    = binmw(idxint);                     % reduced mass values vector of bin i
   t(i)      = numel(rbinint);                    % reduced number of signals in bin i
   % create new vectors intns and mw
   if i==1
       nintns = rbinint;
       nmw    = rbinmw;
   else
       nintns = [nintns rbinint]; 
       nmw    = [nmw rbinmw];
   end
end
%%

%% *************************************************************************
function [nmw,nintns] = redpeaknmb_old(mw,intns,ithreshperc)
% *************************************************************************
% this is the outdated function from Mar 2019 to reduce peak number by 
% removing low intensity peaks
% bin(i)   - bin number i
% pcrm(i)  - percentage of data to be removed in the i-th bin
ithresh    = ithreshperc/100; 
bin(1,:)   = [  0   700]; pcrm(1)  = 0.90;
bin(2,:)   = [700   900]; pcrm(2)  = 0.74;
bin(3,:)   = [900  1100]; pcrm(3)  = 0.56;
bin(4,:)   = [1100 1300]; pcrm(4)  = 0.42;
bin(5,:)   = [1300 1500]; pcrm(5)  = 0.32;
bin(6,:)   = [1500 1700]; pcrm(6)  = 0.24;
bin(7,:)   = [1700 1900]; pcrm(7)  = 0.18;
bin(8,:)   = [1900 2100]; pcrm(8)  = 0.12;
bin(9,:)   = [2100 2300]; pcrm(9)  = 0.10;
bin(10,:)   = [2300 2500];pcrm(10) = 0.08;
bin(11,:)  = [2500 2800]; pcrm(11) = 0.06;
bin(12,:)  = [2800 3300]; pcrm(12) = 0.05;
bin(13,:)  = [3300 3800]; pcrm(13) = 0.04;
bin(14,:)  = [3800 4400]; pcrm(14) = 0.03;
bin(15,:)  = [4400 5000]; pcrm(15) = 0.02;
bin(16,:)  = [5000 10000];pcrm(16) = 0.01;

for i=1:size(bin,1)
   idx       = find(mw>=bin(i,1) & mw<bin(i,2));  % indices of peaks with a mw in the i-th bin
   binint    = intns(idx);                        % intensities of respective peaks in the i-th bin
   binmw     = mw(idx);                           % mw values of respective peaks in the i-th bin
   n(i)      = numel(binint);                     % total number of peaks in the i-th bin
   srtint    = sort(binint);                      % sorted intensity vector of the i-th bin
   thresh(i) = srtint(round(pcrm(i)*n(i)));       % intensity threshold in the given bin
   idxint    = binint>thresh(i);                  % indices of intensity values above the bin-specific threshold
   rbinint   = binint(idxint);                    % reduced vector of intensities in the given bin
   rbinmw    = binmw(idxint);                     % reduced vector of intensities in the given bin
   t(i)      = numel(rbinint);                    % reduced number of peaks in the i-th bin
   % create new vectors intns and mw
   if i==1
       nintns = rbinint;
       nmw    = rbinmw;
   else
       nintns = [nintns rbinint]; 
       nmw    = [nmw rbinmw];
   end
end
%%