% This is code for work published here: https://www.biorxiv.org/content/10.1101/2020.10.26.355560v1 % Three-dimensional biofilm growth supports a mutualism involving matrix and nutrient sharing % Heidi A. Arjes, Lisa Willis, Haiwen Gui, Yangbo Xiao, Jason Peters, Carol Gross, Kerwyn Casey Huang % doi: https://doi.org/10.1101/2020.10.26.355560 % This is a custom code designed by Heidi A. Arjes that imports a matrix % including measured GFP and RFP fluorescent data from colony images (from % the "typhoon_image_analysis.m" code. Note that this example from a % 96-well plate with a titration row. The "huge_matrix.mat" file included % contains the raw measurements from the 16 h timepoint in the sacA::GFP % library. % To run properly, run this script from a folder containing % "huge_matrix.mat,", "controls.mat", and "genes.mat." %% %Goal - normalize GFP to RFP signal per plate % start with huge matrix clear all; close all; %% %huge_matrix is a 96x32 matrix generated from typhoon_image_analysis code, %where the columns are organized as follows: %1-4: LB_GFP (Plate 1-4) %5-8: LB_RFP (Plate 1-4) %9-12: LBxyl_GFP (Plate 1-4) %13-16: LBxyl_RFP (Plate 1-4) %17-20: MSgg_GFP (Plate 1-4) %21-24: MSgg_RFP (Plate 1-4) %25-28: MSggxyl_GFP (Plate 1-4) %29-32: MSggxyl_RFP (Plate 1-4) load('huge_matrix.mat'); load('genes.mat'); genesidx = find(strcmp(genes,'NA') == 0); genes_302 = genes(genesidx); %% need to subtract blanks on each plate GFP_matrix = nan(96,16); RFP_matrix = nan(96,16); GFP_matrix = huge_matrix(:,[1:4,9:12,17:20,25:28]); RFP_matrix = huge_matrix(:,[5:8,13:16,21:24,29:32]); for i = 1:16 %go plate by plate, subtract blank from each plate if (mod(i,4)==1) | (mod(i,4)==2)%mod is remainder after division (i = 1,5,9,13 all have mod of 1), plates 1 or 2 blank = GFP_matrix(12,i); GFP_matrix_blanked(:,i) = GFP_matrix(:,i)-blank; blank = RFP_matrix(2,i); RFP_matrix_blanked(:,i) = RFP_matrix(:,i)-blank; elseif (mod(i,4)==3) | (mod(i,4)==0) blank = GFP_matrix(96,i); GFP_matrix_blanked(:,i) = GFP_matrix(:,i)-blank; blank = RFP_matrix(86,i); RFP_matrix_blanked(:,i) = RFP_matrix(:,i)-blank; end end %should I make anything less than 0 blank here? or later? (I don't think it %matters, do now GFP_matrix_blanked(GFP_matrix_blanked<0) = 0; RFP_matrix_blanked(RFP_matrix_blanked<0) = 0; figure subplot(121) imagesc(GFP_matrix_blanked); colorbar title('Blanked GFP values'); subplot(122) imagesc(RFP_matrix_blanked); colorbar title('Blanked RFP values'); %% % for figure, just plot MSgg plate 1 figure; subplot(3,1,1) imagesc(reshape(GFP_matrix_blanked(:,9),12,8)'); colorbar title('Blanked GFP values - MSgg P1') subplot(3,1,2) imagesc(reshape(RFP_matrix_blanked(:,9),12,8)'); colorbar title('Blanked RFP values - MSgg P1') % %divide GFP by RFP % % GRRatios = cat(2,(huge_matrix(:,1:4)./ huge_matrix(:,5:8)),(huge_matrix(:,9:12)./ huge_matrix(:,13:16)),... % (huge_matrix(:,17:20)./ huge_matrix(:,21:24)),(huge_matrix(:,25:28)./ huge_matrix(:,29:32))); GRRatios = GFP_matrix_blanked ./RFP_matrix_blanked; % % for figure, just plot MSgg plate 1; subplot(3,1,3); imagesc(reshape(GRRatios(:,9),12,8)'); colorbar; title('GFP ./ RFP'); %% plot all 16 plates figure; imagesc(GRRatios); colorbar; %% Grabbed the titration row data and compiled it % titrations are 100% GFP to 100% RFP in 10% steps - they are in wells % A2:A12 (2:12) in P1 and P2 and H2:H12 (86:96) GRRatio_titrations = nan(11,16); % 11 rows of titrations, 16 plates for col = 1:16 %specifies which rows are titration if(mod(col,4) == 1 || mod(col,4) == 2) GRRatio_titrations(:,col) = GRRatios(2:12,col); %rowA else GRRatio_titrations(:,col) = GRRatios(86:96,col); %rowH end end figure; imagesc(GRRatio_titrations); colorbar; %% %look at titration row and fit a curve %GFPtoRFPtitration = data_matrix(1:11,1:4) ./ data_matrix(1:11,5:8); ratiofun = @(fit_con,x)(fit_con(1)*x./(1-fit_con(2)*x));%forcing through zero, use with blanked data % linear function that goes through 0 % linearfun = @(fit_con,x)(fit_con(1)*x+fit_con(2)); % fit: 1. make a guess and use that in the nlinfit function % 2. call nlinfit, which will return the fit par5ameters xdata = [0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0]; conditions = {'LB' 'LB-xyl' 'MSgg' 'MSgg-xyl'}; mapped_raw = nan(96,16); guess = [1 1 0]; %going through each condition (1:4 = LB P1:P4; 5:8 = LBX P1:P4; 9:12 = MSgg %P1:P4; 13:16 = MSggxyl P1:P4), normalizing each plate to the titration row %on that plate for i = 1:16 % switch to 9 for just MSgg Plate 1 for figure 1 if (mod(i,4)==1)%mod is remainder after division (i = 1,5,9,13 all have mod of 1) figure; %make new figure when we get to a new condition end if mod(i,4)==0 %mod is remainder after division (i = 4,8,12,16 all have mod of 0) subplot(1,4,4); % this is Plate4 from the condition and should be plotted on 4th place title('Plate 4'); else subplot(1,4,mod(i,4)); %remainder is 2 or 3, sublot in 1,2 or 3 (plate 1,2,3) title(['Plate',num2str(mod(i,4))]); end real_fit = nlinfit(xdata,GRRatio_titrations(2:11,i)',ratiofun,guess);%leave out 100% GFP % ask KC how to get r2 value scatter(xdata,GRRatio_titrations(2:11,i)'); hold on; plot(xdata,ratiofun(real_fit,xdata)); % Reconcile with summer code if mod(i,4)==0 title('Plate 4'); else title(['Plate ',num2str(mod(i,4))]); end xlabel('ratio GFP'); ylabel('GFP/RFP intensity (AU)'); ylim([0 10]); %map so that every data value is normalized to the curve on its plate mapped_raw(:,i) = (GRRatios(:,i)-real_fit(3))./(real_fit(1)+real_fit(2).*(GRRatios(:,i)-real_fit(3))); % mapped_raw(:,i)= end mapped_raw(mapped_raw<0) = 0; %% hmm, should I normalize each plate so that average control is 1? % start small by plotting all the controls for LB P1_control = [13,25,37,49,61,73,85]; P2_control = [24,36,48,60,72,84,96]; P3_control = [1,13,25,37,49,61,73]; P4_control = [12,24,36,48,60,72,84]; mapped_scaled = nan(96,16); figure; counter = 1; for i = 1:16 %cycle through mapped_raw if (mod(i,4)==1)%mod is remainder after division (i = 1,5,9,13 all have mod of 1) control = P1_control; subplot(1,4,counter) counter = counter +1; elseif (mod(i,4) ==2) control = P2_control; elseif (mod(i,4) ==3) control = P3_control; elseif (mod(i,4) ==0) control = P4_control; end % take mean of control, find normalization value, multiply entire plate % by scalar scalar = 1/mean(mapped_raw(control,i)); mapped_scaled(:,i) = mapped_raw(:,i) * scalar; plot(mapped_raw(control,i)) hold on xlim([0 8]); ylim([0.5 0.9]); title(conditions{counter-1}); end % double checking controls are now average of 1 figure; counter = 1; for i = 1:16 %cycle through mapped_raw if (mod(i,4)==1)%mod is remainder after division (i = 1,5,9,13 all have mod of 1) control = P1_control; subplot(1,4,counter) counter = counter +1; elseif (mod(i,4) ==2) control = P2_control; elseif (mod(i,4) ==3) control = P3_control; elseif (mod(i,4) ==0) control = P4_control; end plot(mapped_scaled(control,i)) hold on xlim([0 8]); %ylim([0.5 0.9]); title(conditions{counter-1}); end %% plot mapped scaled - 96 by 16 heatmap for figure for paper figure imagesc(mapped_scaled); colorbar xticks([]); yticks([]); title('GFP/RFP mapped to titration row and normalized to controls'); %% reshape mapped_raw to 384 mapped_384 = reshape(mapped_raw,384,4); mapped_scaled_384 = reshape(mapped_scaled,384,4); figure; subplot(1,2,1) imagesc(mapped_384); ax = gca; ax.XTick = [1:4]; ax.XTickLabel = conditions; ax.YTick = [1:size(genes)]; ax.YTickLabel = genes; colorbar; title('mapped to each plate'); subplot(1,2,2) imagesc(mapped_scaled_384); ax = gca; ax.XTick = [1:4]; ax.XTickLabel = conditions; ax.YTick = [1:size(genes)]; ax.YTickLabel = genes; colorbar; title('mapped to each plate & scaled so control = 1'); % conclusion - need to do a bit more comparison (scatter plot, etc) to see % what is better; but my initial instinct is that the scaled one is better, % so using that for now %% % pull out controls and real data controls = [13,25,37,49,61,73,85,120,132,144,156,168,180,192,193,205,217,229,241,253,265,300,312,324,336,348,360,372]; %plot control values for each plate: reshaped_controls = reshape(controls,7,4); %second dimension is controls per plate figure; imagesc(mapped_384(controls,:)) ax = gca; ax.XTick = [1:4]; ax.XTickLabel = conditions; colorbar %% figure; % not scaled to controls for i = 1:4 subplot(2,2,i) histogram(mapped_384(genesidx,i),'Normalization','probability','BinWidth',0.05) hold on; histogram(mapped_384(controls,i),'Normalization','probability','BinWidth',0.05) title(conditions{i}); end %% figure; %scaled so controls are 1 for i = 1:4 subplot(2,2,i) histogram(mapped_scaled_384(genesidx,i),'Normalization','probability','BinWidth',0.085) hold on; histogram(mapped_scaled_384(controls,i),'Normalization','probability','BinWidth',0.085) title(conditions{i}); end