function [output1] =  mismatchbars(varargin)

tic % initialization start
indir = '';
savedir='';
Sname = '';

screen_size=get(0,'ScreenSize');
scrW =screen_size(3);
scrH =screen_size(4);
fh = figure(1);
set(fh, 'Position', [0.01*scrW 0.7*scrH 0.3*scrW 0.5*scrH]);

% length of the read you want to analyze
myL = 70;   % for SxCxxxC library
QUA=0;  % make it 0 if you don't want to do any quality selection

Xnucleotide = '';
SAVE = '';
mismatches = 2;
chunk = 500000;
newline = [char(13) char (10)]; % enforse the Windows type newline!
START = 0;

if exist('varargin','var')
    L = length(varargin);
    if rem(L,2)~=0, error('Parameters/Values must come in pairs.'); end

    % read input variables
    for ni = 1:2:L
        switch lower(varargin{ni})
            case 'name', Sname = varargin{ni+1};
            case 'nuc', Xnucleotide = varargin{ni+1};
            case 'dir', indir=varargin{ni+1};
            case 'qua', QUA=varargin{ni+1};
            case 'length', myL=varargin{ni+1};
            case 'type', type=varargin{ni+1};
            case 'save', SAVE=varargin{ni+1};
            case 'savedir', savedir=varargin{ni+1};
            case 'mismatches', mismatches=varargin{ni+1}; 
            case 'chunk',      chunk = 500000;
        end
    end
end

FID = fopen(fullfile(indir,Sname),'r');

 % find the reading parameters for quick reading
[TOT, lines]=findPAR( FID, chunk );

if TOT>1
    disp([Sname ' has ' num2str(lines) ' lines reading in ~'...
          num2str(TOT) ' chunks of ' num2str(chunk) ]);
else
    fprintf([Sname '; ']);
end


for i=1:TOT
    
    fprintf(['chunk ' num2str(i) '/' num2str(TOT) '; ']);

    [tags, ReadF, ReadR, quaF, quaR, START, isEOF] = ...
                                loadNparse2(indir,Sname,chunk,START,newline);
    if ~isempty(tags)
        Tags{i} = cellstr(tags);
        NucF{i} = cellstr(ReadF);
        NucR{i} = cellstr(ReadR);
        scoreF{i} = cellstr(quaF);
        scoreR{i} = cellstr(quaR);
    end
end

TAG ={};
SeqF = {};
SeqR = {};
QuaF = {};
QuaR = {};

for i=1:size(Tags,2)
    
    TAG = [TAG;Tags{i}];
    SeqF = [SeqF; NucF{i}];
    SeqR = [SeqR; NucR{i}];
    QuaF = [QuaF; scoreF{i}];
    QuaR = [QuaR; scoreR{i}];
end
  
 TAG = char(TAG);                           
 SeqF = char(SeqF);
 SeqR = char(SeqR);
 QuaF = char(QuaF);
 QuaR = char(QuaR);

% filter TAG to retain only TAG with at most four mismatches
  [r1,~]=find(TAG=='*');
  
  r2 = setdiff(1:size(TAG,1), r1);
  
  TAG = TAG(r2',:);
  SeqF = SeqF(r2',:);
  SeqR = SeqR(r2',:);
  QuaF = QuaF(r2',:);
  QuaR = QuaR(r2',:);

  % transfer TAG to make the position of mismatches count from the very beginning
  % of the reads
 TAG = tag_trans(TAG);

if QUA
    disp('selecting an error-free subset');
    disp(  unique(QuaF(1:1000,:))   );
    [rF,~]=find( QuaF>' ' & QuaF<=QUA);
    [rR,~]=find( QuaR>' ' & QuaR<=QUA);
    R=unique([rF; rR]);
    
    keep = setdiff(1:1:size(QuaF,1),R); 
    
    disp(['selected ' num2str(numel(keep)) ' out of ' ...
          num2str(size(QuaF,1)) ' total reads ']);

    SeqF = SeqF(keep,:);
    SeqR = SeqR(keep,:);
    QuaF = QuaF(keep,:);
    QuaR = QuaR(keep,:);
    TAG = TAG(keep,:);
    disp('done');
end

[uTAG, Fr]=uniquef(cellstr(TAG));

cTAG=cellstr(TAG);

% find the location of each unique tag in the original TAG file
[~,IX1,IX2]=intersect(uTAG,cTAG);

counter1 = 0.1*ones(1,myL);
counter2 = 0.1*ones(1,myL);

for i=1:numel(IX2)
    
    tag_now = uTAG{IX1(i)};

    [N,L]=tag2num(tag_now);%N, the position of mismatch
    
    if (L+numel(N))==myL & N~=0 & (numel(N)==1 | numel(N)==2)
    
        READ=SeqF(IX2(i),:); % extract one read
        
        % find all reads that have this tag
        IX3 = find ( strcmp(tag_now,cTAG) );
        READS = cellstr( SeqF(IX3, : ));
        
        matchs = cell2mat( regexp(READS,'[ATCG]{1,100}','start') )-1;
        read   =           regexp(READ, '[ATCG]{1,100}','match');
        
        % unique shifts and number of times they are found
        [M, Mn] = uniquef(matchs);  
        
        if numel(N)==1
            for j=1:numel(M)
                errors = N-M(j);
                
                if errors <= size(counter1,2)
                
                    for jj=1:numel(errors)
                    
                        counter1(errors(jj)) = counter1(errors(jj)) + Mn(j);

                    end
                end
            end  
            
            bar(log10(counter1)); 
            
        else
            continue 
        end
            

        Xaxis={};  % tick labels for X axis
        if isempty(Xnucleotide)
            Xnucleotide = read{1};
        end
        
        for j=1:numel(Xnucleotide)
            Xaxis{j}=Xnucleotide(j);
        end

        disp([ read{1} ' ' tag_now ' [' num2str(M') '] ' ...
               num2str(errors(1)) ' ' ...
               num2str( Fr(IX1(i)) ) ]);

        set(gca,'Yscale','linear','Xtick',1:myL,'XTickLabel',Xaxis,...
            'Xlim',[0 myL+1], 'TickDir', 'out', 'Ylim', [-1 6] );
        set(gca,'yTick', 0:2:6);
        set(gca, 'YTickLabel', {'10^0','10^2','10^4','10^6'});
        drawnow;
    end

    
    
end


IX = find ( strcmp(['L{' num2str(myL)  '}M{0}'],uTAG)  );


for i=1
    
    subplot(1,1,i);
    title([Sname ...
           ' reads total: ' num2str(sum(Fr)+size(r1,1))...
           'error-free: '   num2str(Fr(IX)) ]);
           
    line([0 myL+1],[log10(Fr(IX)) log10(Fr(IX))], 'Color', 'k');
    line([0 myL+1],[log10(sum(Fr)) log10(sum(Fr))], 'Color', 'k');
end

output1 = [sum(Fr) Fr(IX) counter1 ];
output2 = [sum(Fr) Fr(IX) counter2 ];

for i=1:numel(uTAG)
    
    if (L+numel(N))==64 & N~=0
        
        IX=strcmp(uTAG{i},cTAG);
        A = min(find(IX)); % location of the first read with this TAG

        M = numel(find(SeqF(A,:)=='~')); % number of '~' preceeding the read
        
        errors = N-M;
        
        for j=1:numel(errors)
            counter(errors(j)) = counter(errors(j)) + Fr(i);
        end
    end
    
end

if ~isempty(SAVE)
    set(gcf,'PaperPositionMode','auto')
    print(fh, '-dtiff', fullfile(savedir, ['A-' SAVE '.tiff' ]) );
    print(fh, '-depsc',  fullfile(savedir, ['A-' SAVE '.eps'  ]) );    
end

end
