SCN Data Analysis Tools

by Duncan K. Foley and Nicholas C. Foley
Assisted by Alexis Hatzis

Read this first:

This notebook, December 22, 2010, documents the various Mathematica routines used in the cluster analysis of SCN images reported in Foley, et al, Bioluminescence.

This notebook presumes that the user has a basic knowledge of Mathematica, including particularly file handling and the Options mechanism. Mathematica provides several tutorials in the Help menu.

This notebook assumes that the user has already prepared a normalized .data file using the companion notebook SCNTIFFProcessingToolsPublic201012.

Start here with normalized data:

To process a file, choose Evaluate Initialization Cells from the Evaluation menu and execute the following ClusterPlot[] cell with the appropriate Options selected by the user. The ClusterPlot[] function prompts the user to choose a normalized data file for processing. The name of this file will be used to label the output.

The Options available are:
BinDim (which determines the size of the mega-pixels used in clustering);
DistanceFunction->CorrelationDistance (the default)
or
DistanceFunction->CosineDistance
(or any other built-in Mathematica DistanceFunction);
Method (as described in the Mathematica documentation for FindClusters[]);
ClusterNumber->"Free" (the default, which allows FindClusters[] to determine the number of clusters according to the Method chosen)
or
CusterNumber->number (which constrains the number of clusters);
ClusterLimit->number (default is 5, which determines how many of the clusters are plotted);
ClusterColorFunction->NicoSCNcf (the default)
or
ClusterColorFunction-> any other color function built-in to Mathematica or defined by the user;
WritePDF->True (the default, which writes the ClusterPlot[] to the hard drive in a file chosen by the user)
or
WritePDF->False (which simply displays the ClusterPlot[] in the notebook).

ClusterPlot[DistanceFunction->CosineDistance,ClusterNumber->5,Method->{Automatic, "SignificanceTest"->{"Gap",Tolerance->3}},BinDim->{10,10},ClusterLimit->5]

Cluster Analysis Functions and Documentation

Peaks and troughs

To find the peaks, select the frames that are brighter than their neighbors over a window, default length 15 frames:

SCNAvTS[] transposes the data and averages over space:

Clear[SCNAvTS];
SCNAvTS[data_]:=SCNAvTS[data]=Mean[Flatten[Transpose[data,{3,1,2}],1]]

Clear[SCNPeaksTroughs];
SCNPeaksTroughs[data_,exwin_:15]:=
Module[{l1,l2,maxl,minl},
maxl={};
minl={};
Do[
If[Max[data[[Range[Max[i-exwin,1],Min[i+exwin,Length[data]]]]]]==data[[i]],maxl=Append[maxl,i]];If[Min[data[[Range[Max[i-exwin,1],Min[i+exwin,Length[data]]]]]]==data[[i]],minl=Append[minl,i]],
{i,1,Length[data]}];
{Transpose[{maxl,data[[maxl]]//N}],Transpose[{minl,data[[minl]]//N}]}
]

To test, use the average data series, but the program works for any time series:

Dimensions[NormData]

SCNPeaksTroughs[SCNAvTS[NormData]]

Clear[SCNPlotPeaksTroughs];
SCNPlotPeaksTroughs[data_,exwin_:15,opts:OptionsPattern[{ListPlot}]]:=
Module[{maxl,minl},
{maxl,minl}=SCNPeaksTroughs[data,exwin];
ListPlot[{data,maxl,minl},FilterRules[{opts},Options[ListPlot]],Joined->{True,False,False},PlotStyle->{Black,Red,Blue}]
]

SCNPlotPeaksTroughs[SCNAvTS[NormData],15]

Cluster finding

Cluster analysis is an unsupervised classification algorithm that takes a group of objects and a "dissimilarity function" as its basic inputs. For example, in SCN research the objects may be time series (vectors) at individual pixels of the image, average time series over mega-pixels of the image, or average time series over "regions of interest" representing particular identified neurons. In its simplest form the number of clusters is specified by the user, and the cluster analysis algorithm classifies the objects into that number of clusters so as to maximize the similarity of objects in each cluster and the differences between clusters. In the SCN data this results in the identification of regions of the image which have similar time patterns over the circadian cycle, according to whatever "dissimilarity function" is specified.

It is also possible to allow the cluster analysis program to find the number of clusters, by specifying a "Tolerance" parameter. In this mode the program calculates the clusters for k=1,2,... clusters and for each of these calculates a statistic (the "gap" statistic described in Tibshirani et al) and increases the number of clusters until the gap statistic passes the Tolerance test statistically. (There are a variety of other algorithms for choosing the number of clusters in use.)

The default dissimilarity function is CorrelationDistance, which for two vectors u,v is 1-(u-Mean[u]).(v-Mean[v])/(Norm[u-Mean[u]]Norm[v-Mean[v]]). Because of the normalizations, CorrelationDistance is insensitive to the amplitude of the vectors; in the SCN context this means it measures only the phase relation between two vectors. We also use CosineDistance, 1-u.v/(Norm[u]Norm[v]), which is sensitive to amplitude. (Note that CorrelationDistiance is the CosineDistance between u-Mean[u] and v-Mean[v].) Another distance function we use is BrayCurtisDistance, Total[Abs[u-v]]/Total[Abs[u+v]], which emphasizes the amplitude of the vectors even more than CosineDistance.

This program finds clusters using CorrelationDistance as the default "dissimilarity function", which can be overridden by the option DistanceFunction->CosineDistance (or whatever distance function the user chooses). The ClusterNumber Option allows the programmer to specify the number of clusters, n. The default ClusterNumber->"Free" allows the program to determine the number of clusters given the Method option.

The program then sorts the clusters according to their maximum brightnesses over all the frames of the image.

Clear[SCNClusters];
Options[SCNClusters]={ClusterNumber->"Free"};
SCNClusters[clusterdata_,opts:OptionsPattern[{SCNClusters,FindClusters}]]:=SCNClusters[clusterdata,opts]=
{Sort[If[OptionValue[ClusterNumber]==="Free",FindClusters[Flatten[MapIndexed[(#1->#2)&,clusterdata,{2}],1],FilterRules[{opts},Options[FindClusters]],DistanceFunction->CorrelationDistance],FindClusters[Flatten[MapIndexed[(#1->#2)&,clusterdata,{2}],1],OptionValue[ClusterNumber],FilterRules[{opts},Options[FindClusters]],DistanceFunction->CorrelationDistance]],(Max[ClusterTS[clusterdata,#1]]>Max[ClusterTS[clusterdata,#2]])&],Flatten[{opts},{1}]}

Cluster mean and standard deviation time series, peaks and troughs

The SCNClusters[] and related display algorithms use routines that calculate the mean and standard deviation time series, peaks and troughs, and peak to trough differences for clusters. A cluster is defined as a list of locations in an image.

Clear[ClusterTS,ClusterPeaks,ClusterTroughs,ClusterAmplitude];
ClusterTS[clusterdata_,cluster_]:=Mean[Extract[clusterdata,cluster]];
ClusterTSSD[clusterdata_,cluster_]:=StandardDeviation[Extract[clusterdata,cluster]];
ClusterPeaks[clusterdata_,cluster_]:=SCNPeaksTroughs[ClusterTS[clusterdata,cluster]][[1]];
ClusterTroughs[clusterdata_,cluster_]:=SCNPeaksTroughs[ClusterTS[clusterdata,cluster]][[2]];
ClusterAmplitude[clusterdata_,cluster_]:=Abs[ClusterTS[clusterdata,cluster][[ClusterPeaks[clusterdata,cluster][[1]]]]-ClusterTS[clusterdata,cluster][[ClusterTroughs[clusterdata,cluster][[1]]]]];

Compressing the data

Clustering on a pixel-by-pixel basis is extremely computationally intensive. It is convenient to be able to compress the image data into "mega-pixels".

The first routine compresses the data in a file to the bindim = {binwidth, binheight}. It is much faster to do cluster analysis on binned ("coarse-grained") data. This version of SCNCompress[] assumes that there is no first frame with the average image.

Clear[SCNCompress];
SCNCompress[data_,bindim_]:=SCNCompress[data,bindim]=Map[N[Mean[Flatten[#,1]]]&,Partition[Transpose[data,{3,1,2}],bindim],{2}]

SCNCompress[NormData,{10,10}]//Dimensions

Example

As an example, this uses SCNClusters[] with CosineDistance, and the image compressed into 10 x 10 megapixels. Notice that the option "DistanceFunction->CosineDistance" overrides the default, which is CorrelationDistance.

SCNClusters[SCNCompress[NormData,{10,10}],DistanceFunction->CosineDistance,Method->{Automatic, "SignificanceTest"->{"Gap",Tolerance->3}}];//Timing

Each cluster is a list of mega-pixel locations:

SCNClusters[SCNCompress[NormData,{10,10}],DistanceFunction->CosineDistance,Method->{Automatic, "SignificanceTest"->{"Gap",Tolerance->3}}][[1,1]]

Displaying Clusters

Color Function

In order to display the clusters, we need to choose colors for them:

The following function chooses colors from the standard Mathematica ColorData[] collection in order to maximize contrast of the clusters. The convention is to use the color SCNcf[clusternumber] for all plots that refer to that clusternumber.

Clear[SCNcf];
SCNcf=ColorData[45]

We have made the following color spectrum for clusters the standard:

Clear[NicoSCNcf];
NicoSCNcf[n_]=If[n≤1,ColorData[3,1],If[1<n≤10,ColorData[3][n],If[10<n≤21,ColorData[15][n],ColorData[1][n]]]]

The ClusterDisplay[] function

ClusterDisplay[] creates a graphics figure showing the clusters in contrasting colors. The colors are controlled by the Option ClusterColorFunction. The Option ClusterLimit allows control of the number of clusters plotted.

Clear[ClusterDisplay];
Options[ClusterDisplay]={ClusterColorFunction->NicoSCNcf,ClusterLimit->6};
ClusterDisplay[clusters_,opts:OptionsPattern[{Graphics,ClusterDisplay}]]:=With[{clusterseries=First[clusters][[1;;Min[{Length[First[clusters]],OptionValue[ClusterLimit]}]]]},Graphics[Raster[Reverse[SparseArray[Flatten[MapIndexed[Thread[#1->#2[[1]]]&,clusterseries]],Automatic,White]],Automatic,ColorFunction->OptionValue[ClusterColorFunction]],FilterRules[{opts},Options[Graphics]],ImageSize->Medium]]

ClusterDisplay[SCNClusters[SCNCompress[NormData,{10,10}],DistanceFunction->CosineDistance,Method->{Automatic, "SignificanceTest"->{"Gap",Tolerance->3}}],ClusterLimit->4,ImageSize->Small]

Standard Deviation and Mean Time Series Plots

The SCNMeanTSPlot[] plots the time series for a group of clusters, with the color keyed to the number of the cluster.

Clear[SCNMeanTSPlot];
Options[SCNMeanTSPlot]={TSCurveThickness->.005};
SCNMeanTSPlot[data_,clusters_,opts:OptionsPattern[{ListPlot,FindClusters,SCNMeanTSPlot,ClusterDisplay,SCNClusters}]]:=With[{clusterseries=First[clusters][[1;;Min[{Length[First[clusters]],OptionValue[ClusterLimit]}]]]},ListPlot[Map[ClusterTS[data,#1]&,clusterseries],FilterRules[{opts},Options[ListPlot]],AxesLabel->{"frame","brightness"},PlotStyle->MapIndexed[{Thickness[OptionValue[TSCurveThickness]],OptionValue[ClusterColorFunction][#2[[1]]]}&,clusterseries],Joined->True,PlotRange->Full]]

SCNMeanTSPlot[SCNCompress[NormData,{10,10}],SCNClusters[SCNCompress[NormData,{10,10}],DistanceFunction->CosineDistance,Method->{Automatic, "SignificanceTest"->{"Gap",Tolerance->3}}],ClusterLimit->4]

SCNSDTSPlot[] plots the standard deviation of the cluster:

Clear[SCNSDTSPlot];
SCNSDTSPlot[data_,clusters_,opts:OptionsPattern[{ListPlot,FindClusters,SCNMeanTSPlot,ClusterDisplay,SCNClusters}]]:=With[{clusterseries=First[clusters][[1;;Min[{Length[First[clusters]],OptionValue[ClusterLimit]}]]]},ListPlot[Map[ClusterTSSD[data,#1]&,clusterseries],FilterRules[{opts},Options[ListPlot]],AxesLabel->{"frame","SDBrightness"},PlotStyle->MapIndexed[{Thickness[OptionValue[TSCurveThickness]],OptionValue[ClusterColorFunction][#2[[1]]]}&,clusterseries],Joined->True,PlotRange->Full]]

SCNSDTSPlot[SCNCompress[NormData,{10,10}],SCNClusters[SCNCompress[NormData,{10,10}],DistanceFunction->CosineDistance,Method->{Automatic, "SignificanceTest"->{"Gap",Tolerance->3}}],ClusterLimit->5]

The ClusterPlot[] function

ClusterPlot[] displays the map image of the clusters together with a time series plot for the brightness of the clusters and the standard deviation of brightness with colors keyed to the colors in the display of the image. The option ClusterLimit allows the programmer to specify the number of clusters plotted. The option BinDim controls the compression into megapixels. All options are passed automatically to the lower-level programs, so the user needs to specify options only at the time of executing the ClusterPlot[] function.

Clear[ClusterPlot];
Options[ClusterPlot]={BinDim->{10,10},WritePDF->True};
ClusterPlot[opts:OptionsPattern[{ListPlot,ClusterDisplay,SCNClusters,SCNMeanTSPlot,ClusterPlot,FindClusters}]]:=
With[{data =Flatten[ReadList[(FileName=SystemDialogInput["FileOpen",WindowTitle->"Please choose a .normdata file"])],1]},With[{clusterdata = SCNCompress[data,OptionValue[BinDim]],clusters=SCNClusters[SCNCompress[data,OptionValue[BinDim]],FilterRules[{opts},Options[SCNClusters]],FilterRules[{opts},Options[FindClusters]]],filename=FileBaseName[FileName]},
If[OptionValue[WritePDF],GraphicsColumn[{Export[SystemDialogInput["FileSave",filename<>".pdf"],#],#}],#]&[GraphicsGrid[{{GraphicsRow[{Style[filename<>"\n"<>"First "<>ToString[OptionValue[ClusterLimit]]<>" of "<>ToString[Dimensions[First[clusters]][[1]]]<>" Clusters\n"<>"BinDim->"<>ToString[OptionValue[BinDim]]<>"\n"<>ToString[Last[clusters]],FontFamily->"Times",Bold]}],SpanFromLeft,SpanFromLeft},{ClusterDisplay[clusters,FilterRules[{opts},Options[ClusterDisplay]],ImageSize->Small],
SCNMeanTSPlot[clusterdata,clusters,FilterRules[{opts},Options[ListPlot]],Joined->True,FilterRules[{opts},Options[SCNMeanTSPlot]],FilterRules[{opts},Options[ClusterDisplay]],ImageSize->Small],
SCNSDTSPlot[clusterdata,clusters,FilterRules[{opts},Options[ListPlot]],Joined->True,FilterRules[{opts},Options[SCNMeanTSPlot]],FilterRules[{opts},Options[ClusterDisplay]],ImageSize->Small]}}]]]]

Reading normalized data from the hard drive manually

The user can read in a .normdata file by executing the following cell, which will prompt for the file to read in, and allow the user to inspect the data.

NormData =Flatten[ReadList[FileName=SystemDialogInput["FileOpen"]],1];

Dimensions[NormData]

FileName

FileBaseName[FileName]

NormData[[15]]//MatrixPlot

Precision

We restrict computations to machine precision to save memory and time.

$MaxPrecision=$MachinePrecision