%sort sequencing images generated by seq-o-matic. Sequence files(folders)
%should be named by sample name followed by 'seq' followed by an index
%(usually automaticaly generated by mm.
folders=dir('*seq*');
folders=sort_nat({folders.name});
for i=1:length(folders)
    rootdir=cd(folders{i});
    poslist=dir('*Pos*');
    poslist=sort_nat({poslist.name});
    
    parfor n=1:length(poslist)
        cd(poslist{n});
        files=dir('*.tif');
        files={files.name};
        info=imfinfo(files{1});
        %parse file names
        tcz=cell(length(files),3);
        for m=1:length(files)
            tcz(m,:)=textscan(files{m},'%*s %u %s %u','Delimiter',{'_','.'});
            %find sequencing channels
            
        end
        tcz(:,2)=cellfun(@cell2mat,tcz(:,2),'UniformOutput',false);
        
        %sort channels into GTAC + other channels. Needs work to
        %generalize.
        ch=unique(tcz(:,2)); 
        %assuming ch is ACGT, sort them.
        if length(ch)==4
            ch=ch([3 4 1 2]);
        else
            ch=ch([3 4 1 2 5:end]);
        end
        
        
        %make max projectios
        for m=1:length(ch)
            idx=find(strcmp(tcz(:,2),ch(m)));
            im=zeros(info.Height,info.Width,length(idx));
            for q=1:length(idx)
                im(:,:,q)=imread(files{idx(q)});
            end
            if m==1
                imwrite(uint16(max(im,[],3)),['../../MAX_',folders{i},'_',poslist{n},'_seq',num2str(i),'.tif']);
            else
                imwrite(uint16(max(im,[],3)),['../../MAX_',folders{i},'_',poslist{n},'_seq',num2str(i),'.tif'],'WriteMode','Append');
            end
        end
        cd ..
    end
    cd ..
end
%% separate by Pos
for i=1:length(poslist)
    mkdir(poslist{i});
    files=dir(['*',poslist{i},'*.tif']);
    files={files.name};
    for n=1:length(files)
        movefile(files{n},[poslist{i},'/',files{n}]);
    end
end
    

%% Fix and align max proj files and make rgb files
%load chprofile %if not present, use chprofile =[];



for i=5:9:length(poslist)
    cd(poslist{i});
    %xlightsubtractseqbgn('MAX',0.65)
    %xlightseqalignment1('bgnsub');
    [chprofile,tforms]=mmseqalignment('MAX',[]);
    cd aligned
    %make rgb. Needs work
    rgbout1('aligned',3);
    mkdir('../RGB');
    movefile('RGB*.tif','../RGB/');
    cd ..
    cd ..
    
end

%% manually select cells and basecall. Not needed for generalization.
for m=1:length(poslist)
    rootdir=cd(poslist{m});
    cd RGB;
    files=dir('RGB*.tif');
    files=sort_nat({files.name});
    I=imread(files{1});
    for i=2:length(files)
        I(:,:,:,i)=imread(files{i});
    end
    implay(I,2);
    im=I(:,:,:,1);
    
    cd ..
    
    
    xall=[];yall=[];
    figure;imshow(im);hold on;title('Press "Enter" when done.');
    while 1
        [x,y]=myginput(1,'crosshair');
        if numel(x)==1
            xall=[xall;x];
            yall=[yall;y];
            scatter(x,y,10,'r','filled');
        else
            break
        end
    end
    close all;
    
    
    cd aligned;
    files=dir('aligned*.tif');
    files=sort_nat({files.name});
    im=conv2(imread(files{1}),eye(3));
    for n=2:4
        im(:,:,n)=conv2(imread(files{1},n),eye(3));
    end
    for i=2:length(files)
        for n=1:4
            im(:,:,n,i)=conv2(imread(files{i},n),eye(3));
        end
    end
    
    sig=zeros(length(xall),1,size(im,3),size(im,4));
    %basecall all selected points
    for i=1:length(xall)
        sig(i,1,:,:)=im(uint16(yall(i)),uint16(xall(i)),:,:);
    end
    sig=permute(squeeze(sig),[1 3 2]);
    [intensity,seq]=max(sig,[],3);
    qual=intensity./sqrt(sum(sig.^2,3));
    seqC=char(seq);
    seqC(seq==1)='G';
    seqC(seq==2)='T';
    seqC(seq==3)='A';
    seqC(seq==4)='C';
    
    save('seq.mat','seq','seqC','qual','intensity','sig');
    cd(rootdir);
end

%% plot quality and intensity over cycles
for i=1:length(poslist)
    rootdir=cd(poslist{i});
    cd aligned;
    load seq.mat
    figure
    subplot(1,2,1);
    errorbar(log2(1+mean(intensity)),((log2(1+mean(intensity))-log2(1+std(intensity)))/sqrt(size(intensity,1))),'LineWidth',2,'Color','k');
    set(gca,'ylim',[0 14]);
    %plot(log2(1+mean(intensity)),'LineWidth',2,'Color','k');
    xlabel('Cycles');
    ylabel('Intensity(bits)');
    title(poslist{i});
    subplot(1,2,2);
    errorbar(mean(qual),std(qual)/sqrt(size(qual,1)),'Linewidth',2,'Color','k');
    set(gca,'ylim',[0.5 1.05],'ytick',0.5:0.1:1);
    xlabel('Cycles');
    ylabel('Quality core');
    cd(rootdir);
    title(poslist{i});
end

%% remake RGB
for i=1:length(poslist)
    cd(poslist{i});
    cd aligned
    rgbout1('aligned',1);
    mkdir('../RGB');
    movefile('RGB*.tif','../RGB/');
    cd ..
    cd ..
    
end


