Skip to content

Commit

Permalink
updated filepath, removed local imviewer_lsm fig from repo tracking.
Browse files Browse the repository at this point in the history
  • Loading branch information
@AlistairBoettiger
AlistairBoettiger committed Nov 22, 2011
2 parents ea59dc1 + 7cfa0ce commit d5d4806
Show file tree
Hide file tree
Showing 65 changed files with 1,896 additions and 932 deletions.
375 changes: 204 additions & 171 deletions CheckDotUpDown.m 100644 → 100755
View file

Large diffs are not rendered by default.

409 changes: 238 additions & 171 deletions CheckDotUpDown.m~ 100644 → 100755
View file

Large diffs are not rendered by default.

Empty file modified CheckGaussFit.m 100644 → 100755
View file
Empty file.
Empty file modified CheckGaussFit.m~ 100644 → 100755
View file
Empty file.
Empty file modified DepthDots.m 100644 → 100755
View file
Empty file.
Empty file modified DuplicateDots.m 100644 → 100755
View file
Empty file.
Empty file modified DuplicateDots.m~ 100644 → 100755
View file
Empty file.
Empty file modified Readme.txt 100644 → 100755
View file
Empty file.
186 changes: 94 additions & 92 deletions Unsupervised_DotFinding.m 100644 → 100755
View file
@@ -1,90 +1,101 @@
%% Unsupervised_DotFinding.m
%
% Alistair Boettiger Date Begun: 03/10/11
% Levine Lab Last Modified: 04/18/11
% Levine Lab Last Modified: 07/07/11
%

clear all;

tot_time = tic;
% Input options
old_lab = 0;
folder = '/Users/alistair/Documents/Berkeley/Levine_Lab/Projects/mRNA_counting/Data/'; % '/Users/alistair/Documents/Berkeley/Levine_Lab/Projects/Enhancer_Modeling/Data/';
rawfolder = '/Volumes/Data/Lab Data/Raw_Data/2011-05-22/'; % '/Volumes/Data/Lab Data/Raw_Data/02-17-11/'; %% %
old_lab = 0; Es = 0; ver = '';
folder = '/home/alistair/Documents/Research/Projects/mRNA_counting/Data/2011-05-22/'; %2011-06-20/'; % 2011-04_and_earlier/'; % % 2011-05-22/'; % 2011-06-20/'; % '/Users/alistair/Documents/Berkeley/Levine_Lab/Projects/Enhancer_Modeling/Data/';
rawfolder = '/home/alistair/Documents/Research/Raw_Data_Temp/2011-05-22/'; %2011-06-20/'; % 2011-04_and_earlier/'; %'; % 2011-06-20/'; % '/Volumes/Data/Lab Data/Raw_Data/02-17-11/'; %% %
rawfolder = '/media/GRAID/Raw_Data/2011-05-22/'
stackfolder = 's05_MP06/';% 's07_MP08/'; % 's07_MP05Hz/';% 's04_MP10/';% 'MP07Hz/'; % 'MP12Hz/'; % 's02_MP01/';% 's01_MP09/';% 'sna2.8Hz/' ;%'s06_MP10_sna18/'; %'s21_MP07/';% 'MP07Hz/';% 's11_G4B/' % 's06_MP10_sna18/'; % %'s10_bcd1x/';% 's11_bcd6x/'; %'s14_comp_cntrl/'; % 's12_cntrl_2label/'; %'MP02_22C/'; %'MP01_22C/'; % 'MGa1x/'; % 'MP10_22C/'; %'MP05_22C/'; %'YW_ths_sog/'; % 'MP10_22C/'; % % 'MP09_22C/'; % 'MGa2x/'; % 'MGa1x/'; % 'MGa2x/'; % 'MP10_22C_sna_y_c/'; %
fname = 's05_MP06Hz_b', ver = '_v3'% 'MP08Hz_snaD_22C_b'; % 's07_MP05Hz_22C'; ver = '_v2'; % 'MP10Hz_c'; %'MP07Hz_snaD_22C_b' ; ver = '_v3';% 'MP12Hz_snaD_22C_b'; % 's04_MP10Hz'; % 's02_MP01_Hz_22C_b'; % 's01_MP09_Hz_22C_c'; %'sna2.8Hz_snaD_22C'; % 's06_MP10_sna18_b'; % 'MP07het_snaD_22C'; % 'MP07Hz_snaD_22C';%'s11_G4B_LacZ';% 's06_MP10_sna18_b'; % 's05_MP06Hz'; % %'s10_bcd1x';% 's11_bcd6x'; % 's14_comp_cntrl'; Es =1; % 's12_cntrl_2label'; Es = 1; % 'MP09_22C_hb_y_f'; Es = 7; % 'MP02_22C_hb_y'; Es = 9; % 'MP02_22C_hb_y_b'; Es = 10; % % 'MP01_22C_hb_y_f'; Es = 12; % 'MP01_22C_hb_y_c'; Es = 10; % 'MP01_22C_hb_y'; Es = 13; % 'MGa1x_LacZ_b'; Es = 12; % 'MP10_22C_sna_y_e'; Es = 12; % 'MP05_22C_sna_y_c'; Es =7; % 'MP10_22C_sna_y_d3'; Es = 1; %'YW_ths_sog'; Es = 12; % % 'MP09_22C_hb_y_e'; Es = 10; % 'MP09_22C_hb_y_d'; Es=11; % 'MGa2x_LacZ_sna_b'; Es = 10; % 'MP10_22C_sna_y_d'; % 'MGa_LacZ'; %'MGa2x_LacZ_sna'; %'MP10_22C_sna_y_c'; old_lab = 1; % 'MP05_22C_sna_y'; old_lab = 1; %
mRNA_channels = 2;% 2; % 3; % 1; % total mRNA channels

stackfolder = 's14_comp_cntrl/'; % 's12_cntrl_2label/'; %'MP02_22C/'; %'MP01_22C/'; % 'MGa1x/'; % 'MP10_22C/'; %'MP05_22C/'; %'YW_ths_sog/'; % 'MP10_22C/'; % % 'MP09_22C/'; % 'MGa2x/'; % 'MGa1x/'; % 'MGa2x/'; % 'MP10_22C_sna_y_c/'; %
fname ='s14_comp_cntrl'; Es =1; % 's12_cntrl_2label'; Es = 1; % 'MP09_22C_hb_y_f'; Es = 7; % 'MP02_22C_hb_y'; Es = 9; % 'MP02_22C_hb_y_b'; Es = 10; % % 'MP01_22C_hb_y_f'; Es = 12; % 'MP01_22C_hb_y_c'; Es = 10; % 'MP01_22C_hb_y'; Es = 13; % 'MGa1x_LacZ_b'; Es = 12; % 'MP10_22C_sna_y_e'; Es = 12; % 'MP05_22C_sna_y_c'; Es =7; % 'MP10_22C_sna_y_d3'; Es = 1; %'YW_ths_sog'; Es = 12; % % 'MP09_22C_hb_y_e'; Es = 10; % 'MP09_22C_hb_y_d'; Es=11; % 'MGa2x_LacZ_sna_b'; Es = 10; % 'MP10_22C_sna_y_d'; % 'MGa_LacZ'; %'MGa2x_LacZ_sna'; %'MP10_22C_sna_y_c'; old_lab = 1; % 'MP05_22C_sna_y'; old_lab = 1; %
mRNA_channels = 3; %2; % 1; % total mRNA channels
sname = fname; % 'MP07het_snaD_22C_1';% '_1'; % additional label on slide.


ver = '';% '_v2';

% MP10_22C_sna_y_c and MP05_22C all done at 3.5, 4, 0.03, 30, 30
% MGa2x and MGa1x all done at 2.5, 3, 0.03, 30, 30

% Focus on subset of image:
m = 1/2048; % .9;% .7; % .5; .7; % 1/2048; %
Zs = 50; % Upper limit on number of Z sections length(Im);
h = 2048; w=2048;

% xp1= floor(h/2*m)+1; xp2 = floor(h/2*(2-m))+1; yp1 = floor(w/2*m)+1; yp2 = floor(w/2*(2-m))+1;
% hs = yp2-yp1+1; ws = xp2-xp1+1;

hs = 200; ws = 200;
xp1 = 1801; yp1 = 1301; xp2 = xp1 + ws -1; yp2 = yp1 + hs - 1;

emb = '01';
Imax = imread([rawfolder,stackfolder,'max_',fname,'_',emb,'.tif']);
Imax_dots = Imax(xp1:xp2,yp1:yp2,1:3);
figure(2); clf; imagesc(Imax_dots);


%%

disp(['Coordinates: ', num2str(xp1), ' : ', num2str(xp2), ', ' num2str(yp1), ' : ', num2str(yp2) ] );



filename = [rawfolder,'/',fname];
load([rawfolder,stackfolder,sname,'.mat'])

getpreciseZ = 0;
consec_layers = 2;
ovlap = 3;
w = Datas.Stack1.Image1.IMG.width;
h = Datas.Stack1.Image1.IMG.height;
if Es==0
Zs = Datas.LSM_info.DimensionZ;
Es = length(fields(Datas)) - 3; % Number of Stacks
end
% ------- Option: Focus on subset of image: ------------------- %
m = 1/2048; % .7; % .5; .7; % 1/2048; %
xp1= floor(h/2*m)+1; xp2 = floor(h/2*(2-m))+1; yp1 = floor(w/2*m)+1; yp2 = floor(w/2*(2-m))+1;
hs = yp2-yp1+1; ws = xp2-xp1+1;

% ws = 2048; hs = 2048; xp1 = 1; yp1 = 1;
% xp2 = xp1 + ws -1; yp2 = yp1 + hs - 1;
disp(['Coordinates: ', num2str(xp1), ' : ', num2str(xp2), ', ' num2str(yp1), ' : ', num2str(yp2) ] );
% ------------------------------------------------------------- %


show_projected = 1; % show max-project with all dots and linked dots.
% -------------- Graphing and Display Options ------------------ %
show_projected = 1; % show max-project with all dots and linked dots.
plotdata = 0; % CheckDotUpDown display parameter
plotZdata = 0 ;% show z-map of data
showhist = 1; % show histogram of mRNA counts per cell.
showim = 1; % show colorcoded mRNA counts per cell
bins = 40; % bins for histograms of mRNA
t = 0; %.45; % threshold for region definition plotting
spread = 1.3; % over/under
% ------------------------------------------------------------- %


%---- Dot Finding Parameters ----- %
% dotfinder's parameters
sigmaE = 3;% IMPORTANT
sigmaI = 4; % IMPORTANT
% min_int = 0.04; % 5 ;% .05 % not necessary Fix at Zero
min_int = 0.04; % 5 ;% .05 % not necessary Fix at Zero
FiltSize = 30;%
min_size = 30;%
min_peak = 1000; %

% Build the Gaussian Filter
Ex = fspecial('gaussian',FiltSize,sigmaE); % excitatory gaussian
Ix = fspecial('gaussian',FiltSize,sigmaI); % inhibitory gaussian
Filt = Ex -Ix;
% sphere finding parameters
getpreciseZ = 0;
consec_layers = 3;
ovlap = 2;
% large ovlap yields confusing dots and then watershed splits these up
% in weird dot-distructive ways
%---------------------------------%


%Data = cell(10,mRNA_channels);

% Build the Gaussian Filter
Ex = fspecial('gaussian',FiltSize,sigmaE); % excitatory gaussian
Ix = fspecial('gaussian',FiltSize,sigmaI); % inhibitory gaussian

disp('Running DotFinder1');

%%
for e= 1:Es
%%
disp('loading data...');
tic
disp('loading data...');

if e<10
emb = ['0',num2str(e)];
else
emb = num2str(e);
end



try load([rawfolder,stackfolder,fname,'_',emb,'_nucdata.mat']);

Expand All @@ -104,55 +115,46 @@
continue
end
end

% filename = [rawfolder,'/',fname];
% Im = lsm_read_mod([filename,'.mat'],str2double(emb),1.5E4);





toc

% thresh = .1;
disp(['analyzing embryo, ',emb,'...']);

for mRNAchn = 1:mRNA_channels % mRNAchn =2

if mRNAchn == 1;
min_int = 0.05; % just for speed
else
min_int = 0.05; %
end


DotLabels= cell(1,Zs);
DotData = cell(1,Zs);
Inds = cell(1,Zs);
Ints = cell(1,Zs);
im_folder = cell(1,Zs);

tic; disp('finding dots...');
for z = 1:Zs % z = 11
try
im_folder{z} = [rawfolder,stackfolder,fname,'_',emb,'_z',num2str(z),'.tif'];
Iin_z = imread(im_folder{z});
catch meZ
Zs = z-1;
disp(meZ.message);
disp(['stack depth = ',num2str(Zs)]);
for z = 1:Zs % z = 11
im_folder{z} = [rawfolder,stackfolder,fname,'_',emb,'_z',num2str(z),'.tif'];
try
Iin_z = imreadfast(im_folder{z});
[DotLabels{z},DotData{z},Inds{z},Ints{z}] = dotfinder(Iin_z(xp1:xp2,yp1:yp2,mRNAchn),Ex,Ix,min_int,min_size,min_peak);
catch err
disp(err.message);
Zs = z-1;
break
end
[DotLabels{z},DotData{z},Inds{z},Ints{z}] = dotfinder(Iin_z(xp1:xp2,yp1:yp2,mRNAchn),Ex,Ix,min_int,min_size);
end
end
end
toc;

Cents = cell2mat(DotData');
DotData = DotData(1:Zs);
DotMasks = DotMasks(1:Zs);
% resize;
DotLabels= DotLabels(1:Zs);
DotData = DotData(1:Zs);
Inds = Inds(1:Zs);
Ints = Ints(1:Zs);


%%

intype = class(Iin_z);
intype = class(Iin_z);
dotC = CheckDotUpDown(DotLabels,DotData,Inds,Ints,plotdata,getpreciseZ,consec_layers,ovlap,xp1,xp2,yp1,yp2,intype);

Cents = cell2mat(DotData');

% Project all layers

if show_projected == 1
Expand All @@ -165,7 +167,7 @@

Imax_dots = Imax(xp1:xp2,yp1:yp2,mRNAchn);
figure(2);
Iout = figure(2); clf; imagesc(Imax_dots);
Iout = figure(2); clf; imagesc(Imax_dots); colorbar;
colordef black; set(gcf,'color','k');
colormap hot; hold on;
plot( dotC(:,1),dotC(:,2),'w+','MarkerSize',14 );
Expand All @@ -174,7 +176,7 @@
end
%%

clear Imax Cents DotData DotMasks Iin_z
clear Imax Cents DotData DotLabels Inds Ints Iin_z

%%

Expand All @@ -195,7 +197,7 @@
Nmin = single(NucLabel);
Nmin(Nmin==0)=NaN;
Nstart = min(Nmin(:));
Nucs_list = unique(NucLabel);
Nucs_list = nonzeros(unique(NucLabel));
Nnucs = length(Nucs_list);

% M = NucLabel;
Expand All @@ -221,7 +223,7 @@
mRNA_den(i) = mRNA_cnt(i)/imdata.Area(i);
nuc_area(i) = length(imdata.PixelID{i});
if showim == 1
Plot_mRNA(NucLabel==nn) = single(mRNA_den(i));
Plot_mRNA(NucLabeled==nn) = single(mRNA_den(i));
end
end
% normalize density to the average cell area
Expand Down Expand Up @@ -249,13 +251,13 @@
set(gcf,'color','w');
title(['Cell size adjusted mRNA per cell. mean = ',...
num2str(m_den,4),' std=',num2str(s_den,4)]);
saveas(histfig,[folder,fname,'_',emb,'_chn',num2str(mRNAchn),'_hist',ver,'.jpg'],'jpg');
% saveas(histfig,[folder,fname,'_',emb,'_chn',num2str(mRNAchn),'_hist',ver,'.jpg'],'jpg');
% write to disk?
end

if showim == 1
mRNA_map = figure(3); clf; colordef black;
imagesc(Plot_mRNA); colormap('hot'); colorbar;
imagesc(Plot_mRNA*mean(imdata.Area)); colormap('hot'); colorbar;
set(gcf,'color','k');
saveas(mRNA_map,[folder,fname,'_',emb,'_chn',num2str(mRNAchn),'rvar',ver,'.jpg'],'jpg');
end
Expand All @@ -270,18 +272,18 @@
clear imdata M C W mRNA_map Fig_regvar histfig Iout
%
%% Export data
% Data{e,mRNAchn}.nucarea = nuc_area;
% Data{e,mRNAchn}.dotC = dotC;
% Data{e,mRNAchn}.mRNAcnt = mRNA_cnt;
% Data{e,mRNAchn}.Plot_mRNA = Plot_mRNA;
% Data{e,mRNAchn}.mRNAsadj = mRNA_sadj;
% % Data{e,mRNAchn}.DotData = DotData; % break the camel;
% % Data{e,mRNAchn}.DotMasks = DotMasks;
% % Data{e,mRNAchn}.imdata = imdata;
% % Data{e,mRNAchn}.mRNAden = mRNA_den;

Rpars.sigmaE = sigmaE;
Rpars.sigmaI = sigmaI;
Rpars.min_int = min_int;
Rpars.FiltSize = FiltSize;
Rpars.min_size = min_size;
Rpars.getpreciseZ = getpreciseZ;
Rpars.consec_layers = consec_layers;
Rpars.ovlap = ovlap;

save([folder,fname,'_',emb,'_',num2str(mRNAchn),'_data',ver],...
'nuc_area','dotC','mRNA_cnt','Plot_mRNA','mRNA_sadj');
save([folder,fname,'_',emb,'_chn',num2str(mRNAchn),'_data',ver,'.mat'],...
'nuc_area','dotC','mRNA_cnt','Plot_mRNA','mRNA_sadj','Rpars');

clear nuc_area dotC mRNA_cnt mRNA_sadj Plot_mRNA

Expand All @@ -299,10 +301,10 @@
clear Iin_z DotData DotMasks I_max cent1 bw dL Cents ...
Nmin imdata imdata2 NucLabeled Plot_mRNA M C ...
nuc_area dotC mRNA_cnt mRNA_den mRNA_sadj;

%save([folder,fname,'_slidedata',ver], 'Data');

Tout = toc(tot_time)/(60*60);
disp(['elpased time = ',num2str(Tout), ' hours']);

toc(tot_time)
disp('All slide data saved');


Expand Down

0 comments on commit d5d4806

Please sign in to comment.