ceb51b269f971336c8e12de1ec653787e1242575 ceisenhart Mon Feb 1 10:16:04 2016 -0800 Adding in changes recommended by Angie, cleaning things up a bit refs #16341 diff --git src/hg/expMatrixToJson/expMatrixToJson.c src/hg/expMatrixToJson/expMatrixToJson.c index 742feff..bf6a925 100644 --- src/hg/expMatrixToJson/expMatrixToJson.c +++ src/hg/expMatrixToJson/expMatrixToJson.c @@ -42,42 +42,76 @@ " leaf node when hovered over.\n" ); } /* Command line validation table. */ static struct optionSpec options[] = { {"multiThreads", OPTION_BOOLEAN}, {"CSV", OPTION_BOOLEAN}, {"threads", OPTION_INT}, {"memLim", OPTION_INT}, {"descFile", OPTION_STRING}, {"attributeTable", OPTION_STRING}, {NULL, 0}, }; +struct slDoubleInt + { + struct slDoubleInt *next; + double val; + int index; + }; + struct bioExpVector /* Contains expression information for a biosample on many genes. */ { struct bioExpVector *next; char *name; // name of biosample. char *desc; // description of biosample. int count; // Number of genes we have data for. double *vector; // An array allocated dynamically. struct rgbColor color; // Color for this one int children; // Number of bioExpVectors used to build the current + struct slDoubleInt *topGeneIndeces; // The indeces for the top 10 genes that drove the clustering up to this point + int contGenes; //The number of contributing genes }; + +struct slDoubleInt *slDoubleIntNew(double x, int y) +/* Return a new double. */ +{ +struct slDoubleInt *a; +AllocVar(a); +a->val = x; +a->index = y; +return a; +} + +int slDoubleIntCmp(const void *va, const void *vb) +/* Compare two slDoubles. */ +{ +const struct slDoubleInt *a = *((struct slDoubleInt **)va); +const struct slDoubleInt *b = *((struct slDoubleInt **)vb); +double diff = a->val - b->val; +if (diff < 0) + return -1; +else if (diff > 0) + return 1; +else + return 0; +} + double stringToDouble(char *s) /* Convert string to a double. Assumes all of string is number * and aborts on an error. Errors on 'nan'*/ { char* end; double val = strtod(s, &end); if (val != val) errAbort("A value of %f was encountered. Please change this value then re run the program.", val); if ((end == s) || (*end != '\0')) errAbort("invalid double: %s", s); return val; } struct bioExpVector *bioExpVectorListFromFile(char *matrixFile) /* Read a tab-delimited file and return list of bioExpVectors */ @@ -92,30 +126,31 @@ if (vectorSize == 0) { // Detect first row. vectorSize = chopByWhite(line, NULL, 0); // counting up AllocArray(row, vectorSize); continue; } AllocVar(el); AllocArray(el->vector, vectorSize); el->count = chopByWhite(line, row, vectorSize); assert(el->count == vectorSize); int i; for (i = 0; i < el->count; ++i) el->vector[i] = stringToDouble(row[i]); el->children = 1; + el->contGenes = 0; slAddHead(&list, el); } lineFileClose(&lf); slReverse(&list); return list; } int fillInNames(struct bioExpVector *list, char *nameFile) /* Fill in name field from file. */ { struct lineFile *lf = lineFileOpen(nameFile, TRUE); char *line; struct bioExpVector *el = list; int maxSize = 0 ; @@ -155,111 +190,134 @@ rAddLeaf(tree->left, pList); rAddLeaf(tree->right, pList); } } struct slRef *getOrderedLeafList(struct hacTree *tree) /* Return list of references to bioExpVectors from leaf nodes * ordered by position in hacTree */ { struct slRef *leafList = NULL; rAddLeaf(tree, &leafList); slReverse(&leafList); return leafList; } -static void rPrintHierarchicalJson(FILE *f, struct hacTree *tree, int level, double distance) +static void rPrintHierarchicalJson(FILE *f, struct hacTree *tree, int level, double distance, struct slName *geneNames) /* Recursively prints out the elements of the hierarchical .json file. */ { struct bioExpVector *bio = (struct bioExpVector *)tree->itemOrCluster; char *tissue = bio->name; struct rgbColor colors = bio->color; if (tree->childDistance > clLongest) /* In practice the first distance will be the longest, and is used for normalization. */ clLongest = tree->childDistance; int i; for (i = 0; i < level; i++) fputc(' ', f); // correct spacing for .json format // *****LEAVES***** if (tree->left == NULL && tree->right == NULL) // Print the leaf nodes { fprintf(f, "{\"name\":\"%s\",\"kids\":\"0\",\"length\":\"%f\",\"colorGroup\":\"rgb(%i,%i,%i)\"}", tissue, tree->parent->childDistance, colors.r, colors.g, colors.b); return; } else if (tree->left == NULL || tree->right == NULL) errAbort("\nHow did we get a node with one NULL kid??"); + // NOTE: There are no leaves past this point for (i = 0; i < level + 1; i++) fputc(' ', f); distance = tree->childDistance/clLongest; double length = 0; if (tree->parent != NULL) { length = tree->parent->childDistance; } // *****NODES***** ++internalNodes; fprintf(f, "{\"name\":\" \", \"number\":\"%i\", \"kids\":\"%i\", \"tpmDistance\": \"%f\", \"length\": \"%f\", \"colorGroup\": \"rgb(%i," - "%i,%i)\",",internalNodes, bio->children , tree->childDistance, length, colors.r,colors.g,colors.b); + "%i,%i)\",\"contGenes\":\"%i\",\"geneList\": {",internalNodes, bio->children , tree->childDistance, length, colors.r,colors.g,colors.b, bio->contGenes); + + struct slDoubleInt *j; + for (j = bio->topGeneIndeces; j != NULL; j = j->next) + { + struct slName *geneName = slElementFromIx(geneNames, j->index); + fprintf(f,"\"%s\":\"%f\"", geneName->name, j->val); + if (j->next != NULL) fprintf(f, ", "); + + } + fprintf(f , "},"); if (distance != distance) distance = 0; struct rgbColor wTB; struct rgbColor wTBsqrt; struct rgbColor wTBquad; if (distance == 0) { wTB = whiteToBlackRainbowAtPos(0); wTBsqrt = whiteToBlackRainbowAtPos(0); wTBquad = whiteToBlackRainbowAtPos(0); } else { wTB = whiteToBlackRainbowAtPos(distance*.95); wTBsqrt = whiteToBlackRainbowAtPos(sqrt(distance*.95)); wTBquad = whiteToBlackRainbowAtPos(sqrt(sqrt(distance*.95))); } fprintf(f, "\"normalizedDistance\": \"%f\", \"whiteToBlack\":\"rgb(%i,%i,%i)\", \"whiteTo", distance, wTB.r, wTB.g, wTB.b); fprintf(f, "BlackSqrt\":\"rgb(%i,%i,%i)\", \"whiteToBlackQuad\":\"rgb(%i,%i,%i)\",\n", wTBsqrt.r, wTBsqrt.g, wTBsqrt.b, wTBquad.r, wTBquad.g, wTBquad.b); for (i = 0; i < level + 1; i++) fputc(' ', f); fprintf(f, "\"children\":[\n"); - rPrintHierarchicalJson(f, tree->left, level+1, distance); + rPrintHierarchicalJson(f, tree->left, level+1, distance, geneNames); fputs(",\n", f); - rPrintHierarchicalJson(f, tree->right, level+1, distance); + rPrintHierarchicalJson(f, tree->right, level+1, distance, geneNames); fputc('\n', f); // Closes the children block for node objects for (i=0; i < level + 1; i++) fputc(' ', f); fputs("]\n", f); for (i = 0; i < level; i++) fputc(' ', f); fputs("}", f); } -void printHierarchicalJson(FILE *f, struct hacTree *tree) +void printHierarchicalJson(FILE *f, struct hacTree *tree, char *geneNamesFile) /* Prints out the binary tree into .json format intended for d3 * hierarchical layouts */ { if (tree == NULL) { fputs("Empty tree.\n", f); return; } double distance = 0; - rPrintHierarchicalJson(f, tree, 0, distance); + + + struct lineFile *lf = lineFileOpen(geneNamesFile, TRUE); + char *line; + struct slName *geneNames; + AllocVar(geneNames); + while (lineFileNextReal(lf, &line)) + { + struct slName *geneName = newSlName(cloneString(line)); + slAddTail(&geneNames, geneName); + } + lineFileClose(&lf); + rPrintHierarchicalJson(f, tree, 0, distance, geneNames); fputc('\n', f); } double slBioExpVectorDistance(const struct slList *item1, const struct slList *item2, void *extraData) /* Return the absolute difference between the two kids' values. Weight based on how many nodes have been merged * to create the current node. Designed for HAC tree use*/ { verbose(3,"Calculating Distance...\n"); const struct bioExpVector *kid1 = (const struct bioExpVector *)item1; const struct bioExpVector *kid2 = (const struct bioExpVector *)item2; int j; double diff = 0, sum = 0; for (j = 0; j < kid1->count; ++j) { @@ -276,33 +334,54 @@ * value is the average of kids' values. * Designed for HAC tree use*/ { verbose(3,"Merging...\n"); const struct bioExpVector *kid1 = (const struct bioExpVector *)item1; const struct bioExpVector *kid2 = (const struct bioExpVector *)item2; float kid1Weight = kid1->children / (float)(kid1->children + kid2->children); float kid2Weight = kid2->children / (float)(kid1->children + kid2->children); struct bioExpVector *el; AllocVar(el); AllocArray(el->vector, kid1->count); assert(kid1->count == kid2->count); el->count = kid1->count; el->name = catTwoStrings(kid1->name, kid2->name); int i; + int gCount = 0; for (i = 0; i < el->count; ++i) { el->vector[i] = (kid1Weight*kid1->vector[i] + kid2Weight*kid2->vector[i]); + double diff = abs((kid1Weight*kid1->vector[i] - kid2Weight*kid2->vector[i])); + + if (diff > 0.0){ + ++el->contGenes; + ++gCount; + int index = i + 1; + if (gCount <= 10){ + struct slDoubleInt *newGene = slDoubleIntNew(diff, index); + slAddHead(&el->topGeneIndeces, newGene); + slSort(&el->topGeneIndeces, slDoubleIntCmp); + } + else{ + if (el->vector[i] > el->topGeneIndeces->val){ + slPopHead(el->topGeneIndeces); + struct slDoubleInt *newGene = slDoubleIntNew(diff, index); + slAddHead(&el->topGeneIndeces, newGene); + slSort(&el->topGeneIndeces, slDoubleIntCmp); + } + } + } } el->children = kid1->children + kid2->children; return (struct slList *)(el); } void colorLeaves(struct slRef *leafList) /* Assign colors of rainbow to leaves. */ { float total = 0.0; //double purplePos = 0.80; struct slRef *el, *nextEl; /* Loop through list once to figure out total, since we need to normalize */ for (el = leafList; el != NULL; el = nextEl) { @@ -342,40 +421,41 @@ else bio2->color = whiteToBlackRainbowAtPos(normalized*100); //bio2->color = saturatedRainbowAtPos(distance); soFar += normalized; } /* Set first color to correspond to 0, since not set in above loop */ struct bioExpVector *bio = leafList->val; //bio->color = saturatedRainbowAtPos(0); bio->color = whiteToBlackRainbowAtPos(.95); } void convertInput(char *expMatrix, char *descFile, bool csv) /* Takes in a expression matrix and makes the inputs that this program will use. * Namely a transposed table with the first column removed. Makes use of system calls * to use cut, sed, kent utility rowsToCols, and paste (for descFile option). */ { - char cmd1[1024], cmd2[1024]; + char cmd[1024],cmd1[1024], cmd2[1024]; if (csv) /* A sed one liner will convert comma separated values into a tab separated values*/ { char cmd3[1024]; safef(cmd3, 1024, "sed -i 's/,/\\t/g' %s ",expMatrix); verbose(2,"%s\n", cmd3); mustSystem(cmd3); } - + safef(cmd, 1024, "cut -f 1 %s | sed \'1d\' > %s.genes", expMatrix, expMatrix); + mustSystem(cmd); safef(cmd1, 1024, "cat %s | sed '1d' | rowsToCols stdin %s.transposedMatrix", expMatrix, expMatrix); /* Exp matrices are X axis of cell lines and Y axis of transcripts. This causes long Y axis and short * X axis, which are not handled well in C. The matrix is transposed to get around this issue. */ verbose(2,"%s\n", cmd1); mustSystem(cmd1); /* Pull out the cell names, and store them in a separate file. This allows the actual data matrix to * have the first row identify the transcript, then all following rows contain only expression values. * By removing the name before hand the computation was made faster and easier. */ if (descFile) { char cmd3[1024]; safef(cmd2, 1024, "rowsToCols %s stdout | cut -f1 | sed \'1d\' > %s.cellNamesTemp", expMatrix, expMatrix); safef(cmd3, 1024, "paste %s.cellNamesTemp %s > %s.cellNames", expMatrix, descFile, expMatrix); verbose(2,"%s\n", cmd2); mustSystem(cmd2); @@ -384,79 +464,78 @@ } else { safef(cmd2, 1024, "rowsToCols %s stdout | cut -f1 | sed \'1d\' > %s.cellNames", expMatrix, expMatrix); verbose(2,"%s\n", cmd2); mustSystem(cmd2); } } void generateHtml(FILE *outputFile, int nameSize, char* jsonFile) // Generates a new .html file for the dendrogram. Will do some size calculations as well. { char *pageName = cloneString(jsonFile); chopSuffix(pageName); int textSize = 12 - log(nodeCount); - int radius = 540 + 270*log10(nodeCount); + //int radius = 540 + 270*log10(nodeCount); int width = 10 * nodeCount; int height = 10 * nodeCount; int labelLength = 10+nameSize*(15-textSize); if (labelLength > 100) labelLength = 100; fprintf(outputFile,"<!DOCTYPE html>\n"); fprintf(outputFile,"<head>\n"); fprintf(outputFile,"<title>New dendrogram tests</title>\n"); fprintf(outputFile,"<link rel=\"stylesheet\" href=\"http://maxcdn.bootstrapcdn.com/bootstrap/3.3.5/css/bootstrap.min.css\">\n"); fprintf(outputFile,"<script src=\"https://ajax.googleapis.com/ajax/libs/jquery/1.11.3/jquery.min.js\"></script>\n"); fprintf(outputFile,"<script src=\"http://maxcdn.bootstrapcdn.com/bootstrap/3.3.5/js/bootstrap.min.js\"></script>\n"); fprintf(outputFile,"<script src=\"http://d3js.org/d3.v3.min.js\" type=\"text/javascript\"></script>\n"); - fprintf(outputFile,"<script src=\"d3.dendrograms1.js\" type=\"text/javascript\"></script>\n"); + fprintf(outputFile,"<script src=\"d3.dendrograms.js\" type=\"text/javascript\"></script>\n"); fprintf(outputFile,"<div class = \"dropdown\">\n"); fprintf(outputFile," <div id = dropdown>\n"); fprintf(outputFile,"</div>\n"); fprintf(outputFile,"<script>\n"); fprintf(outputFile," function load() {\n"); fprintf(outputFile," var data;\n\n"); - fprintf(outputFile," d3.json(\"testClustersWMeta.json\", function(error,json){\n"); + fprintf(outputFile," d3.json(\"%s\", function(error,json){\n", jsonFile); fprintf(outputFile," if (error) return console.warn(error);\n"); fprintf(outputFile," data = json;\n"); + fprintf(outputFile," d3.dendrogram.makeRadialDendrogram('#dendrogram', data,{\n"); + fprintf(outputFile," });\n"); fprintf(outputFile," d3.dendrogram.makeCartesianDendrogram('#phylogram', data, {\n"); fprintf(outputFile," width: %i,\n", width); fprintf(outputFile," height: %i,\n", height); fprintf(outputFile," });\n\n"); - fprintf(outputFile," d3.dendrogram.makeRadialDendrogram('#dendrogram', data,{\n"); - fprintf(outputFile," radius: %i,\n", radius); - fprintf(outputFile," });\n"); fprintf(outputFile," });\n"); fprintf(outputFile," }\n"); fprintf(outputFile,"</script>\n"); fprintf(outputFile,"<style type=\"text/css\" media=\"screen\">\n"); fprintf(outputFile," body { font-family: \"Helvetica Neue\", Helvetica, sans-serif; }\n"); fprintf(outputFile," td { vertical-align: top; }\n"); fprintf(outputFile,"</style>\n"); fprintf(outputFile,"</head>\n"); fprintf(outputFile,"<body onload=\"load()\">\n"); fprintf(outputFile,"<table>\n"); fprintf(outputFile," <tr>\n"); fprintf(outputFile," <td>\n"); - fprintf(outputFile," <h1>Phylogram</h2>\n"); - fprintf(outputFile," <div id='phylogram'></div>\n"); - fprintf(outputFile," </td>\n"); - fprintf(outputFile," <td>\n"); fprintf(outputFile," <h1>Dendrogram</h1>\n"); fprintf(outputFile," <div id='dendrogram'></div>\n"); fprintf(outputFile," </td>\n"); + fprintf(outputFile," <td>\n"); + fprintf(outputFile," <h1>Phylogram</h2>\n"); + fprintf(outputFile," <div id='phylogram'></div>\n"); + fprintf(outputFile," </td>\n"); fprintf(outputFile," </tr>\n"); fprintf(outputFile,"</table>\n"); fprintf(outputFile,"</body>\n"); fprintf(outputFile,"</html>\n"); carefulClose(&outputFile); } void expData(char *matrixFile, char *outDir, char *descFile) /* Read matrix and names into a list of bioExpVectors, run hacTree to * associate them, and write output. */ { verbose(2,"Start binary clustering of the expression matrix by euclidean distance (expMatrixToJson).\n"); @@ -475,40 +554,42 @@ if (clMultiThreads) { verbose(2,"Using %i threads. \n", clThreads); clusters = hacTreeMultiThread(clThreads, (struct slList *)list, localMem, slBioExpVectorDistance, slBioExpVectorMerge, NULL, NULL); } else { verbose(2,"Using 1 threads. \n"); clusters = hacTreeFromItems((struct slList *)list, localMem, slBioExpVectorDistance, slBioExpVectorMerge, NULL, NULL); } struct slRef *orderedList = getOrderedLeafList(clusters); colorLeaves(orderedList); - printHierarchicalJson(f, clusters); + printHierarchicalJson(f, clusters, catTwoStrings(matrixFile, ".genes")); FILE *htmlF = mustOpen(catTwoStrings(outDir,".html"),"w"); generateHtml(htmlF,size,catTwoStrings(outDir,".json")); // Remove temporary files - char cleanup[1024], cleanup2[1024]; + char cleanup[1024], cleanup2[1024], cleanup3[1024]; safef(cleanup, 1024, "rm %s.cellNames", matrixFile); safef(cleanup2, 1024, "rm %s.transposedMatrix", matrixFile); + safef(cleanup3, 1024, "rm %s.genes", matrixFile); mustSystem(cleanup); mustSystem(cleanup2); + mustSystem(cleanup3); end = clock(); verbose(2,"%lld allocated at end. The program took %f seconds to complete.\n", (long long)carefulTotalAllocated(), (double)(end-begin)/CLOCKS_PER_SEC); verbose(2,"Completed binary clustering of the expression matrix by euclidean distance (expMatrixToJson).\n"); } int main(int argc, char *argv[]) /* Process command line. */ { optionInit(&argc, argv, options); clCSV = optionExists("CSV"); clMultiThreads = optionExists("multiThreads"); clThreads = optionInt("threads", clThreads); clMemLim = optionInt("memLim", clMemLim); clDescFile = optionVal("descFile", clDescFile); if (argc != 3)