Repository 'run_flock'
hg clone https://toolshed.g2.bx.psu.edu/repos/immport-devteam/run_flock

Changeset 1:81f9b44f5242 (2017-02-27)
Previous changeset 0:c1e924e37080 (2017-02-27) Next changeset 2:b6b4d08b6858 (2020-07-17)
Commit message:
Uploaded
added:
run_flock/src/README
run_flock/src/cent_adjust.c
run_flock/src/find_connected.c
run_flock/src/flock1.c
run_flock/src/flock2.c
b
diff -r c1e924e37080 -r 81f9b44f5242 run_flock/src/README
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/run_flock/src/README Mon Feb 27 13:30:23 2017 -0500
[
@@ -0,0 +1,139 @@
+This package contains R code for converting and transforming a binary
+FCS file, using the FCSTrans software, and the C code for running the
+flock1 and flock2 population identification software.
+
+src - contains the C code for flock1, flock2 and cent_adjust
+bin - contains the compiled C code for flock1, flock2 and cent_adjust,
+      plus the R code for using the FCSTrans algorithm for file
+      conversion and data transformation.
+doc - documentation for FCSTrans and FLOCK algorithms
+example - sample data and output from FCSTrans and FLOCK
+
+To run this software sucessfully the code assumes the software was installed
+in the /usr/local/flock directory and that R, plus the Bioconductor flowCore
+module has been installed. If you use the ipconvert.sh shell script, you
+may need to adjust the location of the RScript executable and the location
+of the FCSTrans.R code by editing the shell script.
+
+#############################################################################
+# Overview
+#############################################################################
+ImmPort-FLOCK
+
+FLOCK (FLOw Clustering without K), an automated population discovery tool for
+multidimensional FCM data was designed to specifically take into account the
+unique feature of FCM data and produce objective segregation of cell
+populations.
+
+FLOCK parameter settings can be customized by defining Bins and Density
+Threshold.  The number of bins is an integer specifying the number of
+equal-sized regions the data will be partitioned into on each axis. Increasing
+the number of bins increases the sensitivity to detect rare populations but
+may also result in single populations being divided.  Density Threshold is the
+cut-off value to separate the dense regions from background. It is a floating
+point number that helps define population centers; increasing the threshold
+may help separate major populations but could cause the algorithm to overlook
+rare populations.
+
+#############################################################################
+# Compiling C code
+#############################################################################
+cd bin
+cc -o flock1 ../src/flock1.c ../src/find_connected.c -lm
+cc -o flock2 ../src/flock2.c -lm
+cc -o cent_adjust ../src/cent_adjust.c -lm
+
+#############################################################################
+# FCSTrans
+#############################################################################
+A shell script named ipconvert.sh is included that runs the FCSTrans R
+code for converting and transforming a binary FCS file. The output consists
+of one text file containing the transformed channel intensity values and
+another file containing a list of the FCS parameters.
+
+cd bin
+/usr/local/flock/bin/ipconvert.sh ../example/data/FCS2.fcs
+/usr/local/flock/bin/ipconvert.sh ../example/data/FCS3.fcs
+
+#############################################################################
+# Running flock1 or flock2
+#############################################################################
+Running the FLOCK1 and FLOCK2 algorithms generate 8 output files that have
+generic file names. For this reason, it is recommended that one output
+directory be created for one input file to the program.
+
+cd example/output/FCS2
+/usr/local/flock/bin/flock1 ../../data/FCS2.txt
+
+cd example/output/FCS3
+/usr/local/flock/bin/flock2 ../../data/FCS3.txt
+
+Files created: MFI.txt, percentage.txt, population_id.txt, profile.txt,
+               flock_results.txt, coordinates.txt, population_center.txt
+               and percentage.txt
+
+Usage Information for FLOCK1
+----------------------------------------------------------------------------
+basic mode: flock1 fcs.txt
+advanced_ mode: flock1 fcs.txt num_bin density_index max_num_pop
+
+Usage Information for FLOCK2
+----------------------------------------------------------------------------
+basic mode: flock data_file
+advanced mode 0 (specify maximum # of pops): flock data_file max_num_pop
+advanced mode 1 (without # of pops): flock data_file num_bin density_index
+advanced mode 2 (specify # of pops): flock data_file num_bin density_index
+                                           number_of_pop
+advanced mode 3 (specify both # of pops): flock data_file num_bin density_index
+                                           number_of_pop max_num_pop
+
+FLOCK Output Files
+----------------------------------------------------------------------------
+coordinates.txt:
+    Output is the intensity values for each marker and event
+
+flock_results.txt:
+    A combination of the input file, event identifiers and population
+    identifiers.
+
+MFI.txt:
+    Provides the mean fluorescence intensity for each population for each
+    marker/parameter
+
+population_id.txt:
+    Contains population identifiers (i.e, values from [1 to n] where n is
+    the population assigned to the  corresponding events in the input data
+    file, one identifier per row.)
+
+population_center.txt:
+    Contains the centroid coordinates for each identified population
+
+percentage.txt:
+    Includes the population identifiers and percentage of events within that
+    population (relative to the whole data file)
+
+profile.txt:
+    Displays an expression profile, where the approximate expression level
+    for each marker is assigned a numeric value from 1-4, for each identified
+    population
+
+fcs_properties.txt:
+    Contains the number of events, number of populations, and number of
+    markers, as well as the algorithm parameters used in the  analysis
+
+#############################################################################
+# Running cent_adjust 
+#############################################################################
+Running the cent_adjust algorithm generates 4 output files that have
+generic file names. For this reason, it is recommened that one output
+directory be created for one input file to the program. 
+
+mkdir example/output/FCS2/cent_adjust
+cd example/output/FCS2/cent_adjust
+/usr/local/flock/bin/cent_adjust ../population_center.txt ../coordinates.txt
+
+Files created: MFI.txt, percentage.txt, population_id.txt and profile.txt
+
+Usage Information for cent_adjust
+----------------------------------------------------------------------------
+basic mode: cent_adjust input_center input_data_file
b
diff -r c1e924e37080 -r 81f9b44f5242 run_flock/src/cent_adjust.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/run_flock/src/cent_adjust.c Mon Feb 27 13:30:23 2017 -0500
[
b'@@ -0,0 +1,1009 @@\n+/////////////////////////////////////////////////////////\n+//  Cent_adjust version number and modification history\n+//  ImmPort BISC project\n+//  Author: Yu "Max" Qian\n+//  v1.01: Oct 16, 2009\n+//         Line 899 of the main function:\n+//         Changed kmean_term=1 to kmean_term=2\n+//////////////////////////////////////////////////////////\n+\n+\n+#include <time.h>\n+#include <stdio.h>\n+#include <stdlib.h>\n+#include <math.h>\n+#include <string.h>\n+\n+#define DEBUG 0\n+#define LINE_LEN 1024\n+#define FILE_NAME_LEN 128\n+#define PARA_NAME_LEN 64\n+#define MAX_VALUE 1000000000\n+#define CUBE 0\n+\n+\n+void getctrfileinfo(FILE *f_src_ctr, long *num_clust)\n+{\n+\tint ch=\'\\n\';\n+\tint prev=\'\\n\';\n+\tlong num_rows=0;\n+\n+\twhile ((ch = fgetc(f_src_ctr))!= EOF )\n+    {\n+\t\tif (ch == \'\\n\')\n+        {\n+\t\t\t++num_rows;\n+        }\n+\t\tprev = ch;\n+    }\n+\tif (prev!=\'\\n\')\n+\t\t++num_rows;\n+\t\n+\t*num_clust=num_rows;\n+\t//printf("center file has %ld rows\\n", *num_clust);\n+}\n+\n+/************************************* Read basic info of the source file **************************************/\n+void getfileinfo(FILE *f_src, long *file_Len, long *num_dm, char *name_string, int *time_ID)\n+{\n+\tchar src[LINE_LEN];\n+\tchar current_name[64];\n+\tchar prv;\n+\n+\tlong num_rows=0;\n+\tlong num_columns=0;\n+\tint ch=\'\\n\';\n+\tint prev=\'\\n\';\n+\tlong time_pos=0;\n+\tlong i=0;\n+\tlong j=0;\n+\n+\t\n+\n+\tsrc[0]=\'\\0\';\n+\tfgets(src, LINE_LEN, f_src);\n+\n+\tname_string[0]=\'\\0\';\n+\tcurrent_name[0]=\'\\0\';\n+\tprv=\'\\n\';\n+\n+\twhile ((src[i]==\' \') || (src[i]==\'\\t\')) //skip space and tab characters\n+\t\ti++;\n+\n+\twhile ((src[i]!=\'\\r\') && (src[i]!=\'\\n\')) //repeat until the end of the line\n+\t{\n+\t\tcurrent_name[j]=src[i];\n+\t\t\n+\t\tif ((src[i]==\'\\t\') && (prv!=\'\\t\')) //a complete word\n+\t\t{\n+\t\t\tcurrent_name[j]=\'\\0\';\n+\t\t\t\n+          /* \n+           * Commented out John Campbell, June 10 2010\n+           * We no longer want to automatically remove Time column.\n+           * This column should have been removed by column selection\n+          if (0!=strcmp(current_name,"Time"))\n+            {\n+              num_columns++; //num_columns does not inlcude the column of Time\n+              time_pos++;\n+              strcat(name_string,current_name); \n+              strcat(name_string,"\\t");\n+            }\n+          else\n+            {\n+              *time_ID=time_pos;\n+            }\n+          */\n+\n+           num_columns++;\n+           strcat(name_string,current_name);\n+           strcat(name_string,"\\t");\n+\n+\t\t   current_name[0]=\'\\0\';\n+\t\t   j=0;\t\t\t\n+\t\t}\t\t\n+\t\t\n+\t\tif ((src[i]==\'\\t\') && (prv==\'\\t\')) //a duplicate tab or space\n+\t\t{\n+\t\t\tcurrent_name[0]=\'\\0\';\n+\t\t\tj=0;\n+\t\t}\n+\t\t\n+        if (src[i]!=\'\\t\')\n+\t        j++;\n+\n+\t\t\n+\t\tprv=src[i];\n+\t\ti++;\n+\t}\n+\t\n+\tif (prv!=\'\\t\') //the last one hasn\'t been retrieved\n+\t{\n+\t\tcurrent_name[j]=\'\\0\';\n+      /* \n+       * Commented out John Campbell, June 10 2010\n+       * We no longer want to automatically remove Time column.\n+       * This column should have been removed by column selection\n+      if (0!=strcmp(current_name,"Time"))\n+        {\n+          num_columns++;\n+          strcat(name_string,current_name);\n+          time_pos++;\n+        }\n+      else\n+        {\n+          *time_ID=time_pos;\n+        }\n+      */\n+\n+      num_columns++;\n+      strcat(name_string,current_name);\n+\t}\n+\n+\tif (DEBUG==1)\n+\t{\n+\t\tprintf("time_ID is %d\\n",*time_ID);\n+\t\tprintf("name_string is %s\\n",name_string);\n+\t}\n+\n+\t// # of rows\n+\n+\twhile ((ch = fgetc(f_src))!= EOF )\n+    {\n+\t\tif (ch == \'\\n\')\n+        {\n+\t\t\t++num_rows;\n+        }\n+\t\tprev = ch;\n+    }\n+\tif (prev!=\'\\n\')\n+\t\t++num_rows;\n+\t\n+\t*file_Len=num_rows;\n+\t*num_dm=num_columns; \n+\n+\t//printf("original file size is %ld; number of dimensions is %ld\\n", *file_Len, *num_dm);\n+}\n+\n+\n+\n+///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////\n+/************************************* Read the source file into uncomp_data **************************************/\n+void readsource(FILE *f_src, lo'..b'0;\n+\tlong j=0;\n+\n+\tlong *cluster_id;\n+\tlong *IDmapping; //this is to keep the original populationID of the center.txt\n+\n+\tdouble kmean_term=0;\n+\t\n+\tdouble **cluster_center;\n+\tdouble **orig_data;\n+\tdouble **normalized_data;\n+\t\t\n+/*\n+\t_strtime( tmpbuf );\n+    printf( "Starting time:\\t\\t\\t\\t%s\\n", tmpbuf );\n+\t_strdate( tmpbuf );\n+    printf( "Starting date:\\t\\t\\t\\t%s\\n", tmpbuf );\n+*/\n+\n+\tif (argc!=3)\n+\t{\n+\t\tprintf("usage: cent_adjust input_center input_data_file\\n");       \n+\t\texit(0);\n+\t}\t\n+\t\n+\n+\tf_src_ctr=fopen(argv[1],"r");\t\n+\t\n+\t//read source data\n+\tf_src=fopen(argv[2],"r");\n+\t\n+\tgetfileinfo(f_src, &file_Len, &num_dm, name_string, &time_id); //get the filelength, number of dimensions, and num/name of parameters\n+\n+\trewind(f_src); //reset data file pointer\t\n+\n+\torig_data = (double **)malloc(sizeof(double*)*file_Len);\n+\tmemset(orig_data,0,sizeof(double*)*file_Len);\n+\tfor (i=0;i<file_Len;i++)\n+\t{\n+\t\torig_data[i]=(double *)malloc(sizeof(double)*num_dm);\n+\t\tmemset(orig_data[i],0,sizeof(double)*num_dm);\n+\t}\n+\t\n+\treadsource(f_src, file_Len, num_dm, orig_data, time_id); //read the data;\n+\t\n+\tfclose(f_src);\n+\t/////////////////////////////////////////////////////////////////////////////\n+\tgetctrfileinfo(f_src_ctr, &num_clust); //get how many populations\n+\tnorm_used=0;\n+\tdist_used=0;\n+\tkmean_term=2;  //modified on Oct 16, 2009: changed kmean_term=1 to kmean_term=2\n+\n+\trewind(f_src_ctr); //reset center file pointer\n+\n+\t//read population center\n+\tcluster_center=(double **)malloc(sizeof(double*)*num_clust);\n+\tmemset(cluster_center,0,sizeof(double*)*num_clust);\n+\tfor (i=0;i<num_clust;i++)\n+\t{\n+\t\tcluster_center[i]=(double*)malloc(sizeof(double)*num_dm);\n+\t\tmemset(cluster_center[i],0,sizeof(double)*num_dm);\n+\t}\n+\tfor (i=0;i<num_clust;i++)\n+\t\tfor (j=0;j<num_dm;j++)\n+\t\t\tcluster_center[i][j]=0;\n+\n+\tIDmapping=(long *)malloc(sizeof(long)*num_clust);\n+\tmemset(IDmapping,0,sizeof(long)*num_clust);\n+\n+\treadcenter(f_src_ctr,num_clust,num_dm,cluster_center,IDmapping); //read population center\n+    fclose(f_src_ctr);\n+\n+\t/////////////////////////////////////////////////////////////////////////////\n+\tnormalized_data=(double **)malloc(sizeof(double*)*file_Len);\n+\tmemset(normalized_data,0,sizeof(double*)*file_Len);\n+\tfor (i=0;i<file_Len;i++)\n+\t{\n+\t\tnormalized_data[i]=(double *)malloc(sizeof(double)*num_dm);\n+\t\tmemset(normalized_data[i],0,sizeof(double)*num_dm);\n+\t}\n+\t\n+\ttran(orig_data, file_Len, num_dm, norm_used, normalized_data);\n+\t/************************************************* Compute number of clusters *************************************************/\n+\t\n+\tcluster_id=(long*)malloc(sizeof(long)*file_Len);\n+\tmemset(cluster_id,0,sizeof(long)*file_Len);\n+\n+\tassign_event(normalized_data,num_clust,dist_used,kmean_term,file_Len,num_dm,cluster_id,cluster_center,0);\n+\n+\t\n+\t//show(orig_data,cluster_id,file_Len,num_clust,num_dm,show_data,num_disp,name_string); \n+\tshow(orig_data, cluster_id, file_Len, num_clust, num_dm, name_string, IDmapping);\n+\n+\tf_cid=fopen("population_id.txt","w");\n+\n+\tfor (i=0;i<file_Len;i++)\n+\t\tfprintf(f_cid,"%ld\\n",IDmapping[cluster_id[i]]);\n+\t\t\n+\n+\tfclose(f_cid);\n+ \n+\t//added April 16, 2009\n+\tf_mfi=fopen("MFI.txt","w");\n+\n+\tfor (i=0;i<num_clust;i++)\n+\t{\n+\t\tfprintf(f_mfi,"%ld\\t",IDmapping[i]);\n+\n+\t\tfor (j=0;j<num_dm;j++)\n+\t\t{\n+\t\t\tif (j==num_dm-1)\n+\t\t\t\tfprintf(f_mfi,"%.0f\\n",cluster_center[i][j]);\n+\t\t\telse\n+\t\t\t\tfprintf(f_mfi,"%.0f\\t",cluster_center[i][j]);\n+\t\t}\n+\t}\n+\tfclose(f_mfi);\n+\n+\t//ended April 16, 2009\n+\n+\tfor (i=0;i<num_clust;i++)\n+\t\tfree(cluster_center[i]);\n+\tfree(cluster_center);\n+\t\t\n+\n+\t/********************************************** Release memory ******************************************/\n+  \n+\tfor (i=0;i<file_Len;i++)\n+\t{\n+\t\tfree(orig_data[i]);\t\t\n+\t\tfree(normalized_data[i]);\n+\t}\n+\t\n+\tfree(orig_data);\n+\tfree(normalized_data);\n+\tfree(cluster_id);\n+\tfree(IDmapping);\n+\n+/*\n+\t_strtime( tmpbuf );\n+    printf( "Ending time:\\t\\t\\t\\t%s\\n", tmpbuf );\n+\t_strdate( tmpbuf );\n+    printf( "Ending date:\\t\\t\\t\\t%s\\n", tmpbuf );\n+*/\n+\n+}\n'
b
diff -r c1e924e37080 -r 81f9b44f5242 run_flock/src/find_connected.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/run_flock/src/find_connected.c Mon Feb 27 13:30:23 2017 -0500
[
@@ -0,0 +1,176 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+
+//static const char *rcsid = "$Id: find_connected.c,v 1.1 2008/09/05 21:54:40 rpl Exp $";
+
+int find_connected(int **G, int num_dense_grids, int ndim, int *grid_clusterID);
+void depth_first(int startnode);
+
+void bail(const char *);        /* exits via abort */
+static void check_clusters(int *gcID, int ndense);
+static void merge_cluster(int from, int into);
+  
+
+/* Vars that will not change througout the depth-first recursion.  We
+   store them here to avoid endless replication on the stack. */
+static int **Gr=0;
+static int *gcID = 0;          /* grid cluster IDs */
+static int *cluster_count=0;   /* count of nodes per cluster */
+static int ndense=0;
+static int ndim=0;
+/* cid changes between depth-first searches, but is constant within a
+   single search, so it goes here. */
+static int cid=0;
+
+/* Find connected components in the graph of neighboring grids defined in G.
+ *
+ * Output:
+ *
+ * grid_clusterID[]  -- cluster to which each dense grid was assigned
+ * return value      -- number of clusters assigned.
+ */
+int find_connected(int **G, int n_dense_grids, int num_dm, int *grid_clusterID)
+{
+  int nclust=0;                  /* number of clusters found */
+  int i;
+  int *subfac;
+  int subval=0,nempty=0;
+    int clustid=0;
+  
+  size_t sz = n_dense_grids*sizeof(int); 
+  cluster_count = malloc(sz);
+  if(!cluster_count)
+    bail("find_connected: Unable to allocate %zd bytes.\n");
+  memset(cluster_count,0,sz);
+
+  /* set up the statics that will be used in the DFS */
+  Gr=G;
+  gcID = grid_clusterID;
+  ndense = n_dense_grids;
+  ndim = num_dm;
+
+  
+  for(i=0;i<ndense;++i)
+    grid_clusterID[i] = -1;
+
+  for(i=0;i<ndense;++i) {
+    if(grid_clusterID[i] < 0) {  /* grid hasn't been assigned yet */
+      cid = nclust++;
+      depth_first(i);
+    }
+  }
+
+#ifndef NDEBUG
+  check_clusters(gcID,ndense);
+#endif 
+  
+  /* At this point we probably have some clusters that are empty due to merging.
+     We want to compact the cluster numbering to eliminate the empty clusters. */
+
+  subfac = malloc(sz);
+  if(!subfac)
+    bail("find_connected: Unable to allocate %zd bytes.\n");
+ subval=0;
+ nempty=0;
+
+  /* cluster #i needs to have its ID decremented by 1 for each empty cluster with
+     ID < i.  Precaclulate the decrements in this loop: */
+  for(i=0;i<nclust;++i) {
+    //clustid = grid_clusterID[i];
+    if(cluster_count[i] == 0) {
+      subval++;
+      nempty++;
+    }
+    subfac[i] = subval;
+  }
+
+  /* Now apply the decrements to all of the dense grids */
+  for(i=0;i<ndense;++i) {
+   clustid = grid_clusterID[i];
+    grid_clusterID[i] -= subfac[clustid];
+  }
+
+#ifndef NDEBUG
+  //  check_clusters(gcID,ndense);
+#endif  
+  
+  /* correct the number of clusters found */
+  nclust -= nempty;
+
+  return nclust;
+}
+
+
+/* Do a depth-first search for a single connected component in graph
+ * G.  Start from node, tag the nodes found with cid, and record
+ * the tags in grid_clusterID.  Also, record the node count in
+ * cluster_count.  If we find a node that has already been assigned to
+ * a cluster, that means we're merging two clusters, so zero out the
+ * old cid's node count.
+ *
+ * Note that our graph is constructed as a DAG, so we can be
+ * guaranteed to terminate without checking for cycles.  
+ *
+ * Note2: this function can potentially recurse to depth = ndense.
+ * Plan stack size accordingly.  
+ *
+ * Output:
+ *
+ * grid_clusterID[] -- array where we tag the nodes.
+ * cluster_count[]  -- count of the number of nodes per cluster.
+ */
+   
+void depth_first(int node)
+{
+  int i;
+
+  if(gcID[node] == cid)         // we're guaranteed no cycles, but it is possible to reach a node 
+    return;                     // through two different paths in the same cluster.  This early
+                                // return saves us some unneeded work.
+
+  /* Check to see if we are merging a cluster */
+  if(gcID[node] >= 0) {
+    /* We are, so zero the count for the old cluster. */
+    cluster_count[ gcID[node] ] = 0;  
+    merge_cluster(gcID[node], cid);
+    return;
+  }
+
+  /* Update for this node */
+  gcID[node] = cid;
+  cluster_count[cid]++;
+
+  /* Recursively search the child nodes */
+  for(i=0; i<ndim; ++i)
+    if(Gr[node][i] >= 0)      /* This is a child node */
+      depth_first(Gr[node][i]);
+}
+
+void bail(const char *msg)
+{
+  fprintf(stderr,"%s",msg);
+  abort();
+}
+
+
+static void check_clusters(int *gcID, int ndense)
+{
+  int i;
+
+  for(i=0; i<ndense; ++i)
+    if(gcID[i] < 0) {
+      fprintf(stderr,"faulty cluster id at i= %d\n",i);
+      abort();
+    }
+}
+
+static void merge_cluster(int from, int into)
+{
+  int i;
+
+  for(i=0; i<ndense; ++i)
+    if(gcID[i] == from)
+      gcID[i] = into;
+}
b
diff -r c1e924e37080 -r 81f9b44f5242 run_flock/src/flock1.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/run_flock/src/flock1.c Mon Feb 27 13:30:23 2017 -0500
[
b'@@ -0,0 +1,2400 @@\n+/*****************************************************************************\n+\t\n+\tFLOCK: FLOw cytometry Clustering without K (Named by: Jamie A. Lee and Richard H. Scheuermann) \n+\t\n+\tAuthor: (Max) Yu Qian, Ph.D.\n+\t\n+\tCopyright: Scheuermann Lab, Dept. of Pathology, UTSW\n+\t\n+\tDevelopment: November 2005 ~ Forever\n+\n+\tAlgorithm Status: May 2007: Release 1.0\n+\n+\tUsage: flock data_file\n+\t\t    Note: the input file format must be channel values and the delimiter between two values must be a tab.\n+\n+    Changes made July 23, 2010: made errors to STDERR\n+\tChanges made Nov 4, 2010: added one more error (select_num_bin<min_grid) || (select_num_bin>max_grid) to STDERR;\n+\t                          MAX_GRID changed to 50 as larger than 50 seems not useful for any file we have got\n+\t\n+******************************************************************************/\n+#include <time.h>\n+#include <stdio.h>\n+#include <stdlib.h>\n+#include <math.h>\n+#include <string.h>\n+#include <sys/stat.h>\n+#include <unistd.h>\n+#include <assert.h>\n+\n+\n+#define DEBUG 0\n+#define LINE_LEN 1024\n+#define FILE_NAME_LEN 128\n+#define PARA_NAME_LEN 64\n+#define MAX_VALUE 1000000000\n+#define MIN_GRID 6\n+#define MAX_GRID 50\n+\n+#define NORM_METHOD 2 //2 if z-score; 0 if no normalization; 1 if min-max \n+#define KMEANS_TERM 100\n+//#define MAX_NUM_POP 30\n+\n+\n+int find_connected(int **G, int num_dense_grids, int ndim, int *grid_clusterID);\n+\n+/************* Read basic info of the source file ****************************/\n+void getfileinfo(FILE *f_src, int *file_Len, int *num_dm, char *name_string, int *time_ID)\n+{\n+  char src[LINE_LEN];\n+  char current_name[64];\n+  char prv;\n+\n+  int num_rows=0;\n+  int num_columns=0;\n+  int ch=\'\\n\';\n+  int prev=\'\\n\';\n+  int time_pos=0;\n+  int i=0;\n+  int j=0;\n+  int sw=0;\n+\n+  src[0]=\'\\0\';\n+  fgets(src, LINE_LEN, f_src);\n+\n+  if ((src[0]==\'F\') && (src[1]==\'C\') && (src[2]==\'S\'))\n+\t{\n+\t\tfprintf(stderr,"the correct input format is a tab-delimited txt file, instead of FCS file.\\n");\n+\t\tabort();\n+\t}\n+\n+  name_string[0]=\'\\0\';\n+  current_name[0]=\'\\0\';\n+  prv=\'\\n\';\n+\n+  // skip space and tab characters\n+  while ((src[i]==\' \') || (src[i]==\'\\t\'))\n+    i++;\n+\n+  // repeat until the end of line is reached\n+  while ((src[i]!=\'\\0\') && (src[i]!=\'\\n\') && (src[i]!=\'\\r\'))\n+    {\n+      current_name[j]=src[i];\n+\t\t\n+      if ((src[i]==\'\\t\') && (prv!=\'\\t\')) //a complete word\n+        {\n+          current_name[j]=\'\\0\';\n+\t\t\t\n+          if (0!=strcmp(current_name,"Time"))\n+            {\n+              num_columns++; //num_columns does not inlcude the column of Time\n+              time_pos++;\n+              if (sw) {\n+                  strcat(name_string,"\\t");\n+              }\n+              strcat(name_string,current_name); \n+              sw = 1;\n+            }\n+          else\n+            {\n+              *time_ID=time_pos;\n+            }\n+          \n+          \n+          current_name[0]=\'\\0\';\n+          j=0;\t\t\t\n+        }\t\t\n+\t\t\n+      if ((src[i]==\'\\t\') && (prv==\'\\t\')) //a duplicate tab or space\n+        {\n+          current_name[0]=\'\\0\';\n+          j=0;\n+        }\n+\t\t\n+      if (src[i]!=\'\\t\')\n+        j++;\n+\t\t\n+      prv=src[i];\n+      i++;\n+    }\n+\t\n+  if (prv!=\'\\t\') //the last one hasn\'t been retrieved\n+    {\n+      current_name[j]=\'\\0\';\n+      \n+      if (0!=strcmp(current_name,"Time"))\n+        {\n+          num_columns++;\n+          strcat(name_string,"\\t");\n+          strcat(name_string,current_name);\n+          time_pos++;\n+        }\n+      else\n+        {\n+          *time_ID=time_pos;\n+        }\n+      \n+      \n+    }\n+  if (DEBUG==1)\n+    {\n+      printf("time_ID is %d\\n",*time_ID);\n+      printf("name_string is %s\\n",name_string);\n+    }\n+\n+  //start computing # of rows\n+\n+  while ((ch = fgetc(f_src))!= EOF )\n+    {\n+      if (ch == \'\\n\')\n+        {\n+          ++num_rows;\n+        }\n+      prev = ch;\n+    }\n+  if (prev!=\'\\n\')\n+    ++num_rows;\n+\n+  //added on July 23, 2010\n+  if (num_rows<50)\n+  {\n+    fprintf(stderr,"Number of events '..b'lation_ID=(int*)malloc(sizeof(int)*file_Len);\n+  memset(all_population_ID,0,sizeof(int)*file_Len);\n+\n+  kmeans(normalized_data, num_population, KMEANS_TERM, file_Len, num_dm, all_population_ID, population_center);\n+  show(input_data, all_population_ID, file_Len, num_population, num_dm, para_name_string);\n+\n+  ID2Center_all(input_data,file_Len,num_dm,num_population,all_population_ID,population_center);\n+  \n+\n+  f_cid=fopen("population_id.txt","w");\n+  f_ctr=fopen("population_center.txt","w");\n+  f_out=fopen("coordinates.txt","w");\n+  f_results=fopen("flock_results.txt","w");\n+\n+/*\n+  f_parameters=fopen("parameters.txt","w");\n+  fprintf(f_parameters,"Number_of_Bins\\t%d\\n",num_bin);\n+  fprintf(f_parameters,"Density\\t%f\\n",aver_index);\n+  fclose(f_parameters);\n+*/\n+\n+  for (i=0;i<file_Len;i++)\n+\tfprintf(f_cid,"%d\\n",all_population_ID[i]+1); //all_population_ID[i] changed to all_population_ID[i]+1 to start from 1 instead of 0: April 16, 2009\n+\n+  /*\n+   * New to check for min/max to add to parameters.txt\n+   *\n+  */\n+  \n+  fprintf(f_out,"%s\\n",para_name_string);\n+  //fprintf(f_results,"%s\\tEvent\\tPopulation\\n",para_name_string);\n+  fprintf(f_results,"%s\\tPopulation\\n",para_name_string);\n+  for (i=0;i<file_Len;i++)\n+  {\n+\tfor (j=0;j<num_dm;j++)\n+\t{\n+\t\tif (input_data[i][j] < min) {\n+\t\t\tmin = (int)input_data[i][j];\n+\t\t}\n+\t\tif (input_data[i][j] > max) {\n+\t\t\tmax = (int)input_data[i][j];\n+\t\t}\n+\t\tif (j==num_dm-1)\n+\t\t{\n+\t\t\tfprintf(f_out,"%d\\n",(int)input_data[i][j]);\n+\t\t\tfprintf(f_results,"%d\\t",(int)input_data[i][j]);\n+\t\t}\n+\t\telse\n+\t\t{\n+\t\t\tfprintf(f_out,"%d\\t",(int)input_data[i][j]);\n+\t\t\tfprintf(f_results,"%d\\t",(int)input_data[i][j]);\n+\t\t}\n+\t}\n+\t//fprintf(f_results,"%d\\t",i + 1);\n+\tfprintf(f_results,"%d\\n",all_population_ID[i]+1); //all_population_ID[i] changed to all_population_ID[i]+1 to start from 1 instead of 0: April 16, 2009\n+  }\n+\n+/*\n+  f_parameters=fopen("parameters.txt","w");\n+  fprintf(f_parameters,"Number_of_Bins\\t%d\\n",num_bin);\n+  fprintf(f_parameters,"Density\\t%d\\n",den_t_event);\n+  fprintf(f_parameters,"Min\\t%d\\n",min);\n+  fprintf(f_parameters,"Max\\t%d\\n",max);\n+  fclose(f_parameters);\n+*/\n+\n+  f_properties=fopen("fcs.properties","w");\n+  fprintf(f_properties,"Bins=%d\\n",num_bin);\n+  fprintf(f_properties,"Density=%d\\n",den_t_event);\n+  fprintf(f_properties,"Min=%d\\n",min);\n+  fprintf(f_properties,"Max=%d\\n",max);\n+  fprintf(f_properties,"Populations=%d\\n",num_population);\n+  fprintf(f_properties,"Events=%d\\n",file_Len);\n+  fprintf(f_properties,"Markers=%d\\n",num_dm);\n+  fclose(f_properties);\n+\n+\n+  for (i=0;i<num_population;i++) {\n+\t/* Add if we want to include population id in the output\n+\t*/\n+\tfprintf(f_ctr,"%d\\t",i+1);  //i changed to i+1 to start from 1 instead of 0: April 16, 2009\n+\n+\tfor (j=0;j<num_dm;j++) {\n+\t\tif (j==num_dm-1)\n+\t\t\tfprintf(f_ctr,"%.0f\\n",population_center[i][j]);\n+\t\telse\n+\t\t\tfprintf(f_ctr,"%.0f\\t",population_center[i][j]);\n+\t}\n+  }\n+\n+  \t//added April 16, 2009\n+\tf_mfi=fopen("MFI.txt","w");\n+\n+\tfor (i=0;i<num_population;i++)\n+\t{\n+\t\tfprintf(f_mfi,"%d\\t",i+1);\n+\n+\t\tfor (j=0;j<num_dm;j++)\n+\t\t{\n+\t\t\tif (j==num_dm-1)\n+\t\t\t\tfprintf(f_mfi,"%.0f\\n",population_center[i][j]);\n+\t\t\telse\n+\t\t\t\tfprintf(f_mfi,"%.0f\\t",population_center[i][j]);\n+\t\t}\n+\t}\n+\tfclose(f_mfi);\n+\n+\t//ended April 16, 2009\n+\t\t\t\n+  fclose(f_cid);\n+  fclose(f_ctr);\n+  fclose(f_out);\n+  fclose(f_results);\n+\n+\n+  for (i=0;i<num_population;i++)\n+  {\n+\tfree(population_center[i]);\n+  }\n+  free(population_center);\n+ \n+\n+  for (i=0;i<file_Len;i++)\n+    free(normalized_data[i]);\n+  free(normalized_data);\t\n+\t\n+  free(grid_populationID);\n+\n+  free(cluster_populationID);\n+  free(grid_clusterID);\n+  free(cluster_ID);\n+\n+  for (i=0;i<file_Len;i++)\n+    free(input_data[i]);\n+  free(input_data);\n+\n+  free(grid_ID);\n+  free(population_ID);\n+  free(all_population_ID);\n+  free(eventID_To_denseventID);\n+\t\t\n+  ///////////////////////////////////////////////////////////\n+  printf("Ending time:\\t\\t\\t\\t");\n+  fflush(stdout);\n+  system("/bin/date");\n+\n+  return 0;\n+\n+}\n'
b
diff -r c1e924e37080 -r 81f9b44f5242 run_flock/src/flock2.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/run_flock/src/flock2.c Mon Feb 27 13:30:23 2017 -0500
[
b"@@ -0,0 +1,3405 @@\n+///////////\n+// Changes made:\n+// 1. Added another parameter: number of populations\n+// 2. Added a hierarchical merging step based on density change between centroids of two hyper-regions\n+// 3. Picked the longest dimensions between the population centroids to judge whether the two parts should be merged\n+// 4. Removed checking time parameter\n+// 5. Output error to stderr\n+// 6. Fixed the bug of density threshold always = 3\n+// 7. Added another error (select_num_bin<min_grid) || (select_num_bin>max_grid) to STDERR\n+// 8. Fixed a bug for 2D data by using K=e*K\n+// 9. Added some header files, may not be necessary\n+// 10. Added a lower bound (at least two) for number of populations\n+/***************************************************************************************************************************************\n+\t\n+\tFLOCK: FLOw cytometry Clustering without K (Named by: Jamie A. Lee and Richard H. Scheuermann) \n+\t\n+\tAuthor: (Max) Yu Qian, Ph.D.\n+\t\n+\tCopyright: Scheuermann Lab, Dept. of Pathology, UTSW\n+\t\n+\tDevelopment: November 2005 ~ Forever\n+\n+\tStatus: July 2010: Release 2.0\n+\n+\tUsage: flock data_file\n+\t\t    Note: the input file format must be channel values and the delimiter between two values must be a tab.\n+\n+\n+    \t\n+****************************************************************************************************************************************/\n+\n+#include <time.h>\n+#include <stdio.h>\n+#include <stdlib.h>\n+#include <math.h>\n+#include <string.h>\n+#include <sys/stat.h>\n+#include <unistd.h>\n+#include <assert.h>\n+\n+\n+\n+#define DEBUG 0\n+#define LINE_LEN 1024\n+#define FILE_NAME_LEN 128\n+#define PARA_NAME_LEN 64\n+#define MAX_VALUE 1000000000\n+#define MIN_GRID 6\n+#define MAX_GRID 50\n+#define E_T 1.0\n+\n+#define NORM_METHOD 2 //2 if z-score; 0 if no normalization; 1 if min-max \n+#define KMEANS_TERM 10 \n+#define MAX_POP_NUM 128\n+\n+#ifndef max\n+    #define max( a, b ) ( ((a) > (b)) ? (a) : (b) )\n+#endif\n+\n+#ifndef min\n+    #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )\n+#endif\n+\n+static long **Gr=0;\n+static long *gcID = 0;          /* grid cluster IDs */\n+static long *cluster_count=0;   /* count of nodes per cluster */\n+static long ndense=0;\n+static long ndim=0;\n+/* cid changes between depth-first searches, but is constant within a\n+   single search, so it goes here. */\n+static long cid=0;\n+/* Do a depth-first search for a single connected component in graph\n+ * G.  Start from node, tag the nodes found with cid, and record\n+ * the tags in grid_clusterID.  Also, record the node count in\n+ * cluster_count.  If we find a node that has already been assigned to\n+ * a cluster, that means we're merging two clusters, so zero out the\n+ * old cid's node count.\n+ *\n+ * Note that our graph is constructed as a DAG, so we can be\n+ * guaranteed to terminate without checking for cycles.  \n+ *\n+ * Note2: this function can potentially recurse to depth = ndense.\n+ * Plan stack size accordingly.  \n+ *\n+ * Output:\n+ *\n+ * grid_clusterID[] -- array where we tag the nodes.\n+ * cluster_count[]  -- count of the number of nodes per cluster.\n+ */\n+static void merge_cluster(long from, long into)\n+{\n+  int i;\n+\n+  for(i=0; i<ndense; ++i)\n+    if(gcID[i] == from)\n+      gcID[i] = into;\n+}\n+   \n+void depth_first(long node)\n+{\n+  long i;\n+\n+  if(gcID[node] == cid)         // we're guaranteed no cycles, but it is possible to reach a node \n+    return;                     // through two different paths in the same cluster.  This early\n+                                // return saves us some unneeded work.\n+\n+  /* Check to see if we are merging a cluster */\n+  if(gcID[node] >= 0) {\n+    /* We are, so zero the count for the old cluster. */\n+    cluster_count[ gcID[node] ] = 0;  \n+    merge_cluster(gcID[node], cid);\n+    return;\n+  }\n+\n+  /* Update for this node */\n+  gcID[node] = cid;\n+  cluster_count[cid]++;\n+\n+  /* Recursively search the child nodes */\n+  for(i=0; i<ndim; ++i)\n+    if(Gr[node][i] >= 0)      /* This is a child node */\n"..b'ring);\n+\n+  ID2Center_all(input_data,file_Len,num_dm,num_real_pop,all_population_ID,new_population_center);\n+  \n+\n+  f_cid=fopen("population_id.txt","w");\n+  f_ctr=fopen("population_center.txt","w");\n+  f_out=fopen("coordinates.txt","w");\n+  f_results=fopen("flock_results.txt","w");\n+\n+/*\n+  f_parameters=fopen("parameters.txt","w");\n+  fprintf(f_parameters,"Number_of_Bins\\t%d\\n",num_bin);\n+  fprintf(f_parameters,"Density\\t%f\\n",aver_index);\n+  fclose(f_parameters);\n+*/\n+\n+  for (i=0;i<file_Len;i++)\n+\tfprintf(f_cid,"%ld\\n",all_population_ID[i]+1); //all_population_ID[i] changed to all_population_ID[i]+1 to start from 1 instead of 0: April 16, 2009\n+\n+  /*\n+   * New to check for min/max to add to parameters.txt\n+   *\n+  */\n+  \n+  fprintf(f_out,"%s\\n",para_name_string);\n+  //fprintf(f_results,"%s\\tEvent\\tPopulation\\n",para_name_string);\n+  fprintf(f_results,"%s\\tPopulation\\n",para_name_string);\n+  for (i=0;i<file_Len;i++)\n+  {\n+\tfor (j=0;j<num_dm;j++)\n+\t{\n+\t\tif (input_data[i][j] < min) {\n+\t\t\tmin = (int)input_data[i][j];\n+\t\t}\n+\t\tif (input_data[i][j] > max) {\n+\t\t\tmax = (int)input_data[i][j];\n+\t\t}\n+\t\tif (j==num_dm-1)\n+\t\t{\n+\t\t\tfprintf(f_out,"%d\\n",(int)input_data[i][j]);\n+\t\t\tfprintf(f_results,"%d\\t",(int)input_data[i][j]);\n+\t\t}\n+\t\telse\n+\t\t{\n+\t\t\tfprintf(f_out,"%d\\t",(int)input_data[i][j]);\n+\t\t\tfprintf(f_results,"%d\\t",(int)input_data[i][j]);\n+\t\t}\n+\t}\n+\t//fprintf(f_results,"%ld\\t",i + 1);\n+\tfprintf(f_results,"%ld\\n",all_population_ID[i]+1); //all_population_ID[i] changed to all_population_ID[i]+1 to start from 1 instead of 0: April 16, 2009\n+  }\n+\n+/*\n+  f_parameters=fopen("parameters.txt","w");\n+  fprintf(f_parameters,"Number_of_Bins\\t%ld\\n",num_bin);\n+  fprintf(f_parameters,"Density\\t%d\\n",den_t_event);\n+  fprintf(f_parameters,"Min\\t%d\\n",min);\n+  fprintf(f_parameters,"Max\\t%d\\n",max);\n+  fclose(f_parameters);\n+*/\n+\n+  f_properties=fopen("fcs.properties","w");\n+  fprintf(f_properties,"Bins=%ld\\n",num_bin);\n+  fprintf(f_properties,"Density=%d\\n",den_t_event);\n+  fprintf(f_properties,"Min=%d\\n",min);\n+  fprintf(f_properties,"Max=%d\\n",max);\n+  fprintf(f_properties,"Populations=%ld\\n",num_real_pop);\n+  fprintf(f_properties,"Events=%ld\\n",file_Len);\n+  fprintf(f_properties,"Markers=%ld\\n",num_dm);\n+  fclose(f_properties);\n+\n+  for (i=0;i<num_real_pop;i++) {\n+\t/* Add if we want to include population id in the output\n+\t*/\n+\tfprintf(f_ctr,"%ld\\t",i+1);  //i changed to i+1 to start from 1 instead of 0: April 16, 2009\n+\n+\tfor (j=0;j<num_dm;j++) {\n+\t\tif (j==num_dm-1)\n+\t\t\tfprintf(f_ctr,"%.0f\\n",new_population_center[i][j]);\n+\t\telse\n+\t\t\tfprintf(f_ctr,"%.0f\\t",new_population_center[i][j]);\n+\t}\n+  }\n+\n+  \t//added April 16, 2009\n+\tf_mfi=fopen("MFI.txt","w");\n+\n+\tfor (i=0;i<num_real_pop;i++)\n+\t{\n+\t\tfprintf(f_mfi,"%ld\\t",i+1);\n+\n+\t\tfor (j=0;j<num_dm;j++)\n+\t\t{\n+\t\t\tif (j==num_dm-1)\n+\t\t\t\tfprintf(f_mfi,"%.0f\\n",new_population_center[i][j]);\n+\t\t\telse\n+\t\t\t\tfprintf(f_mfi,"%.0f\\t",new_population_center[i][j]);\n+\t\t}\n+\t}\n+\tfclose(f_mfi);\n+\n+\t//ended April 16, 2009\n+\t\t\t\n+  fclose(f_cid);\n+  fclose(f_ctr);\n+  fclose(f_out);\n+  fclose(f_results);\n+\n+\n+  for (i=0;i<num_population;i++)\n+  \tfree(population_center[i]);\n+  \n+  free(population_center);\n+ \n+  for (i=0;i<num_real_pop;i++)\n+\t  free(new_population_center[i]);\n+  \n+  free(new_population_center);\n+\n+  for (i=0;i<file_Len;i++)\n+    free(normalized_data[i]);\n+  free(normalized_data);\t\n+\t\n+  free(grid_populationID);\n+\n+  free(cluster_populationID);\n+  free(grid_clusterID);\n+  free(cluster_ID);\n+\n+  for (i=0;i<file_Len;i++)\n+    free(input_data[i]);\n+  free(input_data);\n+\n+  free(grid_ID);\n+  free(population_ID);\n+  free(all_population_ID);\n+  free(eventID_To_denseventID);\n+\t\t\n+  ///////////////////////////////////////////////////////////\n+  printf("Ending time:\\t\\t\\t\\t");\n+  fflush(stdout);\n+  system("/bin/date");\n+\n+  /*\n+   * Windows version\n+  _strtime( tmpbuf );\n+  printf( "Ending time:\\t\\t\\t\\t%s\\n", tmpbuf );\n+  _strdate( tmpbuf );\n+  printf( "Ending date:\\t\\t\\t\\t%s\\n", tmpbuf );\n+ */\n+  \n+  return 0;\n+}\n'