Upload New File

Data Analysis/fitLinearInPhase.m

+function dataOut = fitLinearConcSeries_V03_wBW(dataIn)
+% dataOut = fitLinearConcSeries(dataIn)
+% This function does the quality control and fiting of concentration series
+% data for inPhase using a water-in-oil emulsion approach.
+% Here, "dataIn" contains one specific condition like salt pH and temperature
+
+linFit_Eq = 'a*x + b';
+
+% get protein concentrations (one per experiment)
+proteinConcentration = cell2mat(structCrawler(dataIn, 'proteinConcentration')');
+
+% get data (one per emulsion droplet)
+% careful, the Ind are cells in cells!!!
+volumeFraction = structCrawler(dataIn, 'volumeFraction');
+emulsionVolume = structCrawler(dataIn, 'emulsionVolume');
+dropletVolume = structCrawler(dataIn, 'dropletTotalVolume');
+dropletVolumesInd = structCrawler(dataIn, 'dropletVolumes');
+dropletIntensityInd = structCrawler(dataIn, 'dropletIntensities');
+dropletIntensity = structCrawler(dataIn, 'dropletIntensities_mean');
+soluteIntensity = structCrawler(dataIn, 'soluteIntensity_mean');
+dropletVolumeBF = structCrawler(dataIn, 'dropletTotalVolumeBF');
+
+% combine all the data to make it more easy to work on it
+combinedData = cellfun(@(a,b,c,d,e,f) [a,b,c,d,e,f], volumeFraction,...
+                                                     emulsionVolume,...
+                                                     dropletVolume,...
+                                                     dropletIntensity,...
+                                                     soluteIntensity,...
+                                                     dropletVolumeBF,...
+                       'UniformOutput', false);
+
+% data clean up
+% check for outliers in the volume fraction
+for k = 1 : length(combinedData)
+    currData = combinedData{k};
+    % remove emulsions with nan values as volume fraction
+    indx = ~isnan(currData(:,1));
+    currData = currData(indx, :);
+    dropletVolumesInd{k} = dropletVolumesInd{k}(indx);
+    dropletIntensityInd{k} = dropletIntensityInd{k}(indx);
+    % remove outliers based on a statistical anaylis 
+    % catches empty emulsions to high Vfrac (dust, etc.)
+    indx = ~isoutlier(currData(:,1), 'quartiles');
+    currData = currData(indx, :);
+    dropletVolumesInd{k} = dropletVolumesInd{k}(indx);
+    dropletIntensityInd{k} = dropletIntensityInd{k}(indx);
+    % override old data
+    combinedData{k} = currData;
+end
+
+% get proteinConcentrations in the right format
+currSize = cell2mat(cellfun(@size, combinedData, 'UniformOutput', false)');
+proteinConcentrations = repelem(proteinConcentration, currSize(:,1));
+
+% sort by concentration into cells for averageing
+unqProteinC =unique(proteinConcentrations);
+NConc = length(unqProteinC);
+volFracGroup = cell(1, NConc);
+volFracGroupBF = cell(1, NConc);
+emulsionVolumeGroup = cell(1, NConc);
+dropletVolumeGroup = cell(1, NConc);
+partionFacIntGroup = cell(1, NConc);
+proteinConcentrationsGroup = cell(1, NConc);
+for k = 1 : length(unqProteinC)
+   
+    currIndx = proteinConcentration == unqProteinC(k);
+    currDataGroup = vertcat(combinedData{currIndx});
+    volFracGroup{k} =  currDataGroup(:, 1);
+    sortVec = true(size(volFracGroup{k}));
+    volFracGroup{k} =  volFracGroup{k}(sortVec, 1);
+    volFracGroupBF{k} = currDataGroup(sortVec, 6) ./ currDataGroup(sortVec, 2);
+    % only output variables not used for fit
+    emulsionVolumeGroup{k} = currDataGroup(sortVec, 2);
+    dropletVolumeGroup{k} = currDataGroup(sortVec, 3);
+    partionFacIntGroup{k} = currDataGroup(sortVec, 4) ./ currDataGroup(sortVec, 5);
+    proteinConcentrationsGroup{k} = repmat(unqProteinC(k), size(partionFacIntGroup{k}));
+    dropletVolumesIndGroup{k} = dropletVolumesInd(currIndx);
+    try
+        dropletVolumesIndGroup{k} = cat(2, dropletVolumesIndGroup{k}{:});
+        dropletVolumesIndGroup{k} = dropletVolumesIndGroup{k}(sortVec);
+    catch
+        dropletVolumesIndGroup{k} = dropletVolumesIndGroup{k}{end};
+        dropletVolumesIndGroup{k} = dropletVolumesIndGroup{k}(sortVec);
+    end         
+    dropletIntensityIndGroup{k} = dropletIntensityInd(currIndx);
+    try
+        dropletIntensityIndGroup{k} = cat(2, dropletIntensityIndGroup{k}{:});
+        dropletIntensityIndGroup{k} = dropletIntensityIndGroup{k}(sortVec);
+    catch
+        dropletIntensityIndGroup{k} = dropletIntensityIndGroup{k}{end};
+        dropletIntensityIndGroup{k} = dropletIntensityIndGroup{k}(sortVec);
+    end       
+    partionFacIntIndGroup{k} = cellfun(@(x,y) x ./ y , dropletIntensityIndGroup{k}, num2cell(currDataGroup(sortVec, 5)'), 'UniformOutput', false);
+    
+end
+
+% Brightfield
+% remove undetected condensates from BF detection (they are 0)
+indx = cellfun(@(x) x==0, volFracGroupBF, 'UniformOutput', false);
+volFracGroupBF = cellfun(@(x,y) x(~y), volFracGroupBF, indx, 'UniformOutput', false);
+% remove outliers based on a statistical anaylis 
+    % catches empty emulsions to high Vfrac (dust, etc.)
+volFracGroupBF = cellfun(@(x) rmoutliers(x, 'quartiles'), volFracGroupBF, 'UniformOutput', false);
+
+% only use partiton factors (via intensity) of condensate volumes bigger than
+% 100 µm^3
+indx = cellfun(@(x) x > 100, dropletVolumeGroup, 'UniformOutput', false);
+partionFacIntMean = cellfun(@(x,y) x(y), partionFacIntGroup, indx, 'UniformOutput', false);
+
+% group average
+volFracMean = cellfun(@nanmean, volFracGroup);
+volFracMean_Std = cellfun(@nanstd, volFracGroup);
+volFracBFMean = cellfun(@nanmean, volFracGroupBF);
+volFracBFMean_Std = cellfun(@nanstd, volFracGroupBF);
+nDataPoints = cellfun(@length, volFracGroup);
+volFracMean_SEM = volFracMean_Std ./ sqrt(nDataPoints);
+partionFacInt_Std = cellfun(@nanstd, partionFacIntMean);
+partionFacIntMean = cellfun(@nanmean, partionFacIntMean);
+
+% combine all the data to make it more easy to work on it
+volumeFraction = cell2mat(cellfun(@(x)  x(:,1), combinedData, 'UniformOutput', false)');
+
+% routine to check for concentrations that are below c_sat and pass it to fit
+dataExclude = (nDataPoints < 3)';
+dataExcludeNew = (nDataPoints < 3)';
+moreThanOne = true;
+goFit = true;
+while goFit
+    if  sum(dataExclude) ~= length(unqProteinC) - 1
+        currWeight = nDataPoints(~dataExclude);
+        [FitVinVtot, ~] = fit(unqProteinC(~dataExclude), volFracMean(~dataExclude)',  linFit_Eq ,...
+                              'Lower', [0, -1], 'Upper', [1, 0],...
+                              'StartPoint', [0.01, -0.01], 'Weight', currWeight);
+        currCout = abs(FitVinVtot.b / FitVinVtot.a);
+        dataExcludeNew = unqProteinC - currCout;
+        dataExcludeNew = dataExcludeNew < 0;
+    else
+        moreThanOne = false;
+    end
+    if sum(dataExcludeNew) - sum(dataExclude) > 0
+        goFit = true;
+        dataExclude = dataExcludeNew;
+    else
+        goFit = false;
+    end
+end
+
+indxSort = proteinConcentrations >= min(unqProteinC(~dataExclude));
+
+% determine cout and cin via linear fit to Vfrac over ctot
+try
+    if moreThanOne
+        % fit to Vfrac data from all individual emulsion droplets
+        [FitVinVtot, ~] = fit(proteinConcentrations(indxSort), volumeFraction(indxSort),  linFit_Eq ,...
+                              'StartPoint', [0.01, -0.01], 'Robust', 'Bisquare');
+        % fit to mean Vfrac data weighted by number per cot
+        currWeight = nDataPoints(~dataExclude);
+        [~, GoF] = fit(unqProteinC(~dataExclude), volFracMean(~dataExclude)',  linFit_Eq ,...
+                       'Lower', [0, -1], 'Upper', [1, 0],...
+                       'StartPoint', [0.01, -0.01], 'Weight', currWeight);
+        % cout and cin via inPhase equation
+        currA = FitVinVtot.a;
+        currB = FitVinVtot.b;
+        currConf = confint(FitVinVtot, 0.95);
+        currCout = - currB / currA;
+        currCin = (1 - currB ) / currA;
+        % Standard deviation of a and b from 95% confidence interval
+        currA_err =  (currConf(2,1) - currConf(1,1)) / 3.92;
+        currB_err =  (currConf(2,2) - currConf(1,2)) / 3.92;
+        % Error propagation to cout and cin
+        currCout_err = sqrt( (-1/currA)^2 * currB_err^2 + ...
+            (currB/(currA^2))^2 * currA_err^2); 
+        currCin_err = sqrt( (-1/currA)^2 * currB_err^2 + ...
+            ((currB-1)/(currA^2))^2 * currA_err^2);
+        % this would mean a negative cout
+        if currB > 0
+            currA = NaN;
+            currA_err = NaN;
+            currB = NaN;
+            currB_err = NaN;
+            currCin = NaN;
+            currCin_err = NaN;
+            currCout = NaN;
+            currCout_err = NaN;
+        end
+    else
+        currA = NaN;
+        currA_err = NaN;
+        currB = NaN;
+        currB_err = NaN;
+        currCin = NaN;
+        currCin_err = NaN;
+        currCout = NaN;
+        currCout_err = NaN;
+        GoF = NaN;
+        FitVinVtot = NaN;
+    end
+catch
+    currA = NaN;
+    currA_err = NaN;
+    currB = NaN;
+    currB_err = NaN;
+    currCin = NaN;
+    currCin_err = NaN;
+    currCout = NaN;
+    currCout_err = NaN;
+    GoF = NaN;
+    FitVinVtot = NaN;
+    
+end
+
+% determine cout via linear fit to volume fraction via BF
+weightsBF = cell2mat(cellfun(@size, volFracGroupBF, 'UniformOutput', false)');
+weightsBF = weightsBF(:,1)./max(weightsBF(:,1));
+try
+    currX = unique(proteinConcentrations);
+    currY = volFracBFMean';
+    % fit to mean Vfrac data weighted by number per cot
+    [FitVinVtot_BF, GoF_BF] = fit(currX(~dataExclude), currY(~dataExclude),  linFit_Eq ,...
+                                  'Lower', [0, -1], 'Upper', [1, 0],...
+                                  'StartPoint', [0.01, -0.01],...
+                                  'Weights', weightsBF);
+    % cout and cin via inPhase equation
+    currABF = FitVinVtot_BF.a;
+    currBBF = FitVinVtot_BF.b;
+    currCoutBF = - currBBF / currABF;
+    currCinBF = (1 - currBBF ) / currABF;
+    % SD of a and b from 95% confidence interval of fit
+    currConf = confint(FitVinVtot_BF, 0.95);
+    currABF_err =  (currConf(2,1) - currConf(1,1)) / 3.92; 
+    currBBF_err =  (currConf(2,2) - currConf(1,2)) / 3.92;
+    % error propagation of SD for cout and cin
+    currCoutBF_err = sqrt( (-1/currABF)^2 * currBBF_err^2 + ...
+        (currBBF/(currABF^2))^2 * currABF_err^2);
+    currCinBF_err = sqrt( (-1/currABF)^2 * currBBF_err^2 + ...
+        ((currBBF-1)/(currABF^2))^2 * currABF_err^2);
+    % this would mean a negative cout
+    if currBBF > 0
+        currABF = NaN;
+        currABF_err = NaN;
+        currBBF = NaN;
+        currBBF_err = NaN;
+        currCinBF = NaN;
+        currCinBF_err = NaN;
+        currCoutBF = NaN;
+        currCoutBF_err = NaN;
+        
+    end
+    
+catch
+    currABF = NaN;
+    currABF_err = NaN;
+    currBBF = NaN;
+    currBBF_err = NaN;
+    currCinBF = NaN;
+    currCinBF_err = NaN;
+    currCoutBF = NaN;
+    currCoutBF_err = NaN;
+    GoF_BF = NaN;
+    FitVinVtot_BF = NaN;
+end
+
+% derived vadiables
+dataOut.a_slope = currA;
+dataOut.a_slope_err = currA_err;
+dataOut.b_offset = currB;
+dataOut.b_offset_err = currB_err;
+dataOut.Cin = currCin;
+dataOut.Cin_err = currCin_err;
+dataOut.Cout = currCout;
+dataOut.Cout_err = currCout_err;
+dataOut.FitVinVtot = FitVinVtot;
+dataOut.GoF = GoF;
+try
+    dataOut.rsquare = GoF.rsquare;
+    dataOut.adjrsquare = GoF.adjrsquare;
+catch
+    dataOut.rsquare = nan;
+    dataOut.adjrsquare = nan;
+end
+
+dataOut.a_slope_BF = currABF;
+dataOut.a_slope_BF_err = currABF_err;
+dataOut.b_offset_BF = currBBF;
+dataOut.b_offset_BF_err = currBBF_err;
+dataOut.Cin_BF = currCinBF;
+dataOut.Cin_BF_err = currCinBF_err;
+dataOut.Cout_BF = currCoutBF;
+dataOut.Cout_BF_err = currCoutBF_err;
+dataOut.FitVinVtot_BF = FitVinVtot_BF;
+dataOut.GoF_BF = GoF_BF;
+dataOut.PartFac = currCin / currCout;
+dataOut.PartFac_err = sqrt( (1/currCout)^2 * currCin_err^2 + ...
+                            (-currCin/(currCout^2))^2 * currCout_err^2); 
+dataOut.PartFacIntensity = mean(partionFacIntMean);
+dataOut.PartFacIntensity_err = std(partionFacIntMean);
+dataOut.usedConcentrationsFit = (~dataExclude);
+dataOut.usedConcentrationsFitInd = indxSort;
+
+% cleaned raw data
+dataOut.volFracMean = volFracMean;
+dataOut.volFracMean_Std = volFracMean_Std;
+dataOut.volFracMean_SEM = volFracMean_SEM;
+dataOut.unqProteinC = unqProteinC;
+dataOut.volumeFraction = volumeFraction;
+dataOut.proteinConcentrations = proteinConcentrations;
+dataOut.partionFacIntMean = partionFacIntMean;
+dataOut.partionFacInt_Std = partionFacInt_Std;
+dataOut.partionFacIntIndGroup = partionFacIntIndGroup;
+dataOut.dropletIntensityIndGroup = dropletIntensityIndGroup;
+dataOut.dropletVolumesIndGroup = dropletVolumesIndGroup;
+dataOut.volFracGroup = volFracGroup;
+dataOut.emulsionVolumeGroup = emulsionVolumeGroup;
+dataOut.dropletVolumeGroup = dropletVolumeGroup;
+dataOut.partionFacIntGroup = partionFacIntGroup;
+dataOut.proteinConcentrationsGroup = proteinConcentrationsGroup;
+
+% carry on experiment protocol data
+dataOut.saltConcentration = unique(cell2mat(structCrawler(dataIn, 'saltConcentration')'));
+dataOut.temperature = unique(cell2mat(structCrawler(dataIn, 'temperature')'));
+dataOut.labelingFraction = unique(cell2mat(structCrawler(dataIn, 'labelingFraction')'));
+dataOut.pH = unique(cell2mat(structCrawler(dataIn, 'pH')'));
+dataOut.proteinName = cell2mat(unique(structCrawler(dataIn, 'proteinName')'));
+try
+    dataOut.folderExperiment = cell2mat(unique(structCrawler(dataIn, 'folderExperiment')'));
+catch
+    dataOut.folderExperiment = nan;
+end
+
+end
\ No newline at end of file