5
$\begingroup$

I have some data plotted which is of bacterial growth curves, I have 96 sets of growth curves and I am wondering if there is a quick way to extract the straight line area of my growth curves in each case in order to determine their gradient (growth rate) Below I attach an example of one of these plots. The linear part I am concerned with is the growth rate, which is the initial sharp slope region.

Bacterial growth curve

So far my strategy for smaller sets of data is to simply remove the data before and after my linear region and then find the gradient. But seeing as I have quite a bit of data here, is there a way of essentially automating this?

I should say that different curves have the linear region at different points on the x-axis, so I cannot remove data before and after the same point for each one of my curves as different growth curves seem have different amounts of lag time.

I have only been using Mathematica for a short while now, so thank you for your help.

EDIT: here's a link to the data I'm working with: https://www.dropbox.com/s/l04bwn5hifw62ny/Growth180716.xls?dl=0 Each "sub-cell" within each table is a data point.

Here is how I am currently graphing the data, I import the relavent section from my .xls file defined by some Range, and then I generate the x-axis using a Range. Then I Thread to combine the two lists of numbers together and use ListPlot to plot my data:

 GrowthTrip1Range = Range[2, 1208, 9] ;(*Define the range in the .xls for the first line of triplicates*)

GrowthTime = Range[0, 4020, 30]; (*Generating the X-axis (30 minute per point for 72 hours)*)

C1T1 = Import["C:\\Users\\Georgeos\\Desktop\\Growth180716.xls", {"Data", 1,GrowthTrip1Range, 2}];

C1T2 = Import["C:\\Users\\Georgeos\\Desktop\\Growth180716.xls", {"Data", 1, GrowthTrip1Range, 3}];

C1T3 = Import["C:\\Users\\Georgeos\\Desktop\\Growth180716.xls", {"Data", 1, GrowthTrip1Range, 4}];

MeanC1 = Mean[{C1T1, C1T2, C1T3}];

C1Plot = Thread[{GrowthTime, MeanC1}];

ListPlot[C1Plot]
$\endgroup$
  • $\begingroup$ Can you edit your post to include sample data, perhaps putting it on pastebin if it's too long? $\endgroup$ – J. M. will be back soon Jul 22 '16 at 15:28
  • $\begingroup$ @J.M. sure thing $\endgroup$ – Georgeos Hardo Jul 22 '16 at 15:29
  • 2
    $\begingroup$ Actually can you include at least two such datasets, to make sure that any answers are general enough to work on multiple ones? $\endgroup$ – Szabolcs Jul 22 '16 at 15:32
  • 1
    $\begingroup$ A very simple step to take next is to look at ListLinePlot[MeanC1 // Differences, PlotRange -> All]. The noise on the data is low and the slopes calculated point-by-point via Differences can be easily thresholded without fancy PieceWise fitting. $\endgroup$ – John Morganthau Jul 22 '16 at 15:53
  • 1
    $\begingroup$ @mikado RANSAC is in ImageLines btw $\endgroup$ – John Morganthau Jul 22 '16 at 21:13
5
$\begingroup$

Consider the data independent of the abscissa for convenience.

 ListLinePlot[MeanC1 // Differences]

Since the sloped portions are roughly linear the Differences appear simply as offsets when plotted as point-by-point differences. Obviously, we are looking for growth (pun intended) and the slope is positive and so is the offset in the difference plot. Some estimation may be added to algorithmically determine the threshold by I have used a fixed value in this example (0.005). Chosen by "eye."

Chose the points above the threshold.

 growing = Position[MeanC1 // Differences, x_ /; x > .005] // Flatten

Chose the points above the threshold that are "generally" contiguous. The pattern in Position could be changed to "x ==1" to force exact continuity.

 continuousgrowing = Position[Differences[growing], x_ /; x < 2] // Flatten

Plot the selected data with the selected "GrowthTime"

 ListPlot[{GrowthTime[[growing[[continuousgrowing]]]], MeanC1[[growing[[continuousgrowing]]]]}\[Transpose]]
$\endgroup$
  • $\begingroup$ Thank you, John. You choose 0.005 by eye, for more noisy data, do you think increasing the threshold value is an appropriate way to isolate linear regions? $\endgroup$ – Georgeos Hardo Jul 22 '16 at 22:39
  • $\begingroup$ @GeorgeosHardo, Yes. I selected "by-eye" but it may be possible to raise the threshold and to automatically select the level using sorting of the Differences and selecting the top 25% or so of the values since the growth appears to be the "fastest" portion of the trends and occurs for about 25% of the sampled data. $\endgroup$ – John Morganthau Jul 25 '16 at 14:21
6
$\begingroup$

Since measurements are taken at equal intervals, there's no need to incorporate time range into the calculation of slope, we'll just scale the result for the 30-minute duration. It seems that FindClusters responds well to the averaged data, that is if one applies Mean as you do.

(* different for you probably *)
file = FileNameJoin[{NotebookDirectory[], "Growth180716.xls"}];

(* all twelve columns, averaged *)
data = Mean@With[{range = Range[2, 1208, 9]},
    Table[Import[file, {"Data", 1, range, col + 1}], {col, 12}]];

slope[data_] :=
 Module[{clusters, slopes},
  clusters = FindClusters[
    Thread[Differences[data] -> Range[Length[data] - 1]], 3];
  slopes = Table[Module[{i, j},
     {i, j} = MinMax[cluster];
     (data[[j]] - data[[i]])/(j - i)], {cluster, clusters}];
  Max[slopes]]

slope[data]

0.0110031

(* in parts per minute *)
%/30

0.0003667

Added: ImageLines

Different approach utilizes ImageLines which can find the line segment and slope in a somewhat more direct way. For a quick walkthrough, plot is rendered first and then turned into an image.

plot = ListPlot[data,
  Joined -> True,
  Axes -> None,
  PlotRangePadding -> 0,
  ImageSize -> Small]

enter image description here

Plot range of this plot and image dimensions will be needed latter to rescale the slope value obtained from the image.

AbsoluteOptions[plot, PlotRange]

{PlotRange->{{0.,135.},{0.0892083,0.406208}}}

img = ColorNegate@Binarize@Image@plot

enter image description here

This seems to work well for you data. One should tweak the two parameters and MaxFeatures option generally though.

lines = ImageLines[img, .1, .5]

{{{67.3932,112.},{31.3302,0.}}}

Show[img,
 Graphics[{Thick, Orange, Line /@ lines}]]

enter image description here

Module[{x1, x2, y1, y2, w, h},
 {{x1, x2}, {y1, y2}} = PlotRange /. AbsoluteOptions[plot, PlotRange];
 {w, h} = ImageDimensions[img];
 lines[[1]] /. {{a_, b_}, {c_, d_}} :> Divide[
    (d - b)*(y2 - y1)/h,
    (c - a)*(x2 - x1)/w]]

0.0117202

(* in parts per minute *)
%/30

0.000390674

Note the slightly different result (5 % steeper). Different algorithm finds a different line. To a lesser degree, plot with higher image size should change the result minimally too.

Added: NMinimize

Acquiring the slope by optimizing the parameters of a piecewise linear function.

pw[{k1_, k2_, k3_}, {n1_, n2_, n3_}, {x1_, x2_}] :=
 Function[x, Piecewise[{
    {k1 x + n1, x <= x1},
    {k2 x + n2, x <= x2}}, k3 x + n3]]

error[data_, {
   kk_?(VectorQ[#, NumericQ] &),
   nn_?(VectorQ[#, NumericQ] &),
   xx_?(VectorQ[#, NumericQ] &)}] :=
 Module[{values},
  values = pw[kk, nn, xx] /@ Range[Length@data];
  Total[(data - values)^2]]

(* supply sensible starting bounds *)
(* method "RandomSearch" works here *)
(min = NMinimize[error[data, {{k1, k2, k3}, {n1, n2, n3}, {x1, x2}}],
    {{k1, -.001, .001}, {k2, .01, .1}, {k3, -.001, .001},
     {n1, 0, .2}, {n2, -.5, 0}, {n3, .4, .8}, {x1, 20, 40}, {x2, 50, 
      70}},
    MaxIterations -> 5,
    Method -> "RandomSearch"]) // Timing

... sensible output ... it took 6.3 seconds on my computer ...

Visual verification:

Show[
 ListPlot[data],
 Plot[Evaluate[pw[
      {k1, k2, k3}, {n1, n2, n3}, {x1, x2}][x] /. min[[2]]],
  {x, 0, Length@data},
  PlotStyle -> Directive[Thick, Black]],
 ImageSize -> Small]

enter image description here

(* slope value in parts per minute *)
k2/30 /. min[[2]]

0.000380854

$\endgroup$
1
$\begingroup$

Thanks to Mikado for the suggestion. You can get the lines and the slopes for "free" using the following. However work must be done to get relate the image coordinates to the plot axes.

 img = Image[ListLinePlot[C1Plot, Axes -> False]];
 edg = EdgeDetect[GaussianFilter[ColorNegate[Binarize[img]], 2], 
 Method -> {"Canny"}] // DeleteSmallComponents[#, 50] &;
 lines = ImageLines[edg, MaxFeatures -> 2, "Segmented" -> False];
 Show[img, Graphics[{Thickness[.005], Orange, Line /@ lines}]]

 ang[expr_] := ArcTan[(Abs[expr[[2, 2]] - expr[[1, 2]]])/(Abs[
 expr[[2, 1]] - expr[[1, 1]]])]

 angs =  ang /@ lines/Degree
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.