0e9604766086b49eecef30226e15cea726cc14b6
kent
  Wed Nov 24 14:03:14 2021 -0800
Updated usage message, which was missing a new column.

diff --git src/hca/hcaUnifyMatrix/hcaUnifyMatrix.c src/hca/hcaUnifyMatrix/hcaUnifyMatrix.c
index 2346a8b..422917a 100644
--- src/hca/hcaUnifyMatrix/hcaUnifyMatrix.c
+++ src/hca/hcaUnifyMatrix/hcaUnifyMatrix.c
@@ -1,482 +1,483 @@
 /* hcaUnifyMatrix - Given a list of matrices on maybe different gene sets produce 
  * a unified matrix on one gene set. */
 #include "common.h"
 #include "linefile.h"
 #include "hash.h"
 #include "options.h"
 #include "sqlNum.h"
 #include "rangeTree.h"
 #include "vMatrix.h"
 #include "obscure.h"
 #include "../../hg/inc/txGraph.h"
 #include "../../hg/inc/bed12Source.h"
 #include "../../hg/inc/bed.h"
 
 #define BED12_NAME2_SIZE 13
 #define BED_NAME_FIELD 3
 
 double emptyVal = -1;
 
 void usage()
 /* Explain usage and exit. */
 {
 errAbort(
   "hcaUnifyMatrix - Given a list of matrices on maybe different gene sets produce a unified matrix on one gene set\n"
   "usage:\n"
   "   hcaUnifyMatrix graph.txg inList.tsv outMatrix.tsv\n"
   "where\n"
   "   graph.txg is made with txBedToGraph on the mapping.beds in the input\n"
   "   inList.tsv has one line per input we are unifying.  It is tab separated with columns:\n"
   "       1 - name of input to add to column labels\n"
   "       2 - relative priority, lower numbers get attended to first in gene mapping\n"
-  "       3 - input mapping bed file, often result of gencodeVersionForGenes -bed output\n"
-  "       4 - input expression matrix file in tsv format\n"
+  "       3 - scale - number used to scale this input relative to others\n"
+  "       4 - input mapping bed file, often result of gencodeVersionForGenes -bed output\n"
+  "       5 - input expression matrix file in tsv format\n"
   "options:\n"
   "   -bed=mapping.bed\n"
   "   -empty=val - store val (rather than default -1) to show missing data\n"
   "   -trackDb=stanza.ra - save trackDb stanza here\n"
   "   -stats=stats.tsv - save merged stats file here\n"
   );
 }
 
 /* Command line validation table. */
 static struct optionSpec options[] = {
    {"bed", OPTION_STRING},
    {"empty", OPTION_DOUBLE},
    {"trackDb", OPTION_STRING},
    {"stats", OPTION_STRING},
    {NULL, 0},
 };
 
 struct uniInput
 /* A single input to our unification process. */
     {
     struct uniInput *next;  /* Next in singly linked list. */
     char *name;	/* Short, will be reused as label prefix */
     int priority;	/* 1 is attended to first, then 2, etc.  */
     double scaleOut;	/* Everything gets multiplied by this in the output */
     char *mappingBed;	/* Mapping bed file name */
     char *matrixFile;	/* Matrix file name */
     char *colorFile;	/* File it first column sample, last color */
     char *statsFile;	/* File stats are in */
     struct hash *rowHash;	/* keyed by 4th field, values char**  */
     struct hash *bedHash;	/* Also keyed by 4th field */
     struct slRef *rowList;	/* List of all parsed out rows as char** */
     struct memMatrix *matrix;   /* Loaded matrixFile */
     struct hash *matrixHash;	/* Matrix Y value keyed by gene name */
     };
 
 
 struct uniInput *uniInputLoad(char **row)
 /* Load a uniInput from row fetched with select * from uniInput
  * from database.  Dispose of this with uniInputFree(). */
 {
 struct uniInput *ret;
 
 AllocVar(ret);
 ret->name = cloneString(row[0]);
 ret->priority = sqlSigned(row[1]);
 ret->scaleOut = sqlDouble(row[2]);
 ret->mappingBed = cloneString(row[3]);
 ret->matrixFile = cloneString(row[4]);
 ret->colorFile = cloneString(row[5]);
 ret->statsFile = cloneString(row[6]);
 return ret;
 }
 
 struct uniInput *uniInputLoadAll(char *fileName) 
 /* Load all uniInput from a whitespace-separated file.
  * Dispose of this with uniInputFreeList(). */
 {
 struct uniInput *list = NULL, *el;
 struct lineFile *lf = lineFileOpen(fileName, TRUE);
 char *row[7];
 
 while (lineFileRow(lf, row))
     {
     el = uniInputLoad(row);
     slAddHead(&list, el);
     }
 lineFileClose(&lf);
 slReverse(&list);
 return list;
 }
 
 void loadMappingBed(struct uniInput *uni)
 /* Load up a mapping bed file and make up row and bed hashes from it */
 {
 struct hash *rowHash = uni->rowHash = hashNew(0);
 struct hash *bedHash = uni->bedHash = hashNew(0);
 struct lm *lm = rowHash->lm;
 
 struct lineFile *lf = lineFileOpen(uni->mappingBed, TRUE);
 char *row[BED12_NAME2_SIZE];
 while (lineFileRow(lf, row))
     {
     char *name = row[BED_NAME_FIELD];
     struct bed12Source *bed = bed12SourceLoad(row);
     hashAdd(bedHash, name, bed);
     char **rowCopy = lmCloneRow(lm, row, ArraySize(row));
     hashAdd(rowHash, name, rowCopy);
     refAdd(&uni->rowList, rowCopy);
     }
 slReverse(&uni->rowList);
 lineFileClose(&lf);
 }
 
 void loadMatrix(struct uniInput *uni)
 /* Load up a mapping bed file and make up row and bed hashes from it */
 {
 struct memMatrix *matrix = uni->matrix = memMatrixFromTsv(uni->matrixFile);
 struct hash *hash = uni->matrixHash = hashNew(0);
 int y;
 for (y=0; y<matrix->ySize; ++y)
     hashAddInt(hash, matrix->yLabels[y], y);
 }
 
 int countDistinctTypes(struct txGraph *gList)
 /* Just count up how many distinct source types we got to track in all graphs*/
 {
 /* Figure up how many distinct sources we have */
 struct hash *uniq = hashNew(0);
 struct txGraph *g;
 for (g = gList; g != NULL; g = g->next)
     {
     int i;
     for (i=0; i<g->sourceCount; ++i)
         {
 	char *typeName = g->sources[i].type;
 	hashStore(uniq, typeName);
 	}
     }
 int typeCount = uniq->elCount;
 hashFree(&uniq);
 return typeCount;
 }
 
 int countDistinctSourceTypes(struct txGraph *g)
 /* Count up number of sources in one graph */
 {
 struct hash *typeHash = hashNew(0);
 struct txSource *sources = g->sources;
 int count = g->sourceCount;
 int i;
 for (i=0; i<count; ++i)
     {
     char *typeName = sources[i].type;
     hashStore(typeHash, typeName);
     }
 int retVal = typeHash->elCount;
 hashFree(&typeHash);
 return retVal;
 }
 
 
 struct rbTree *rangeTreeOnBed(struct bed *bed)
 /* Return an rbTree holding all the blocks of the bed */
 {
 struct rbTree *t = rangeTreeNew();
 int i;
 for (i=0; i<bed->blockCount; ++i)
     {
     int start = bed->chromStart + bed->chromStarts[i];
     int end = start + bed->blockSizes[i];
     rangeTreeAdd(t, start, end);
     }
 return t;
 }
 
 int bedBedOverlap(struct bed *a, struct bed *b)
 /* Return amount of overlap at exon level between a and b */
 {
 if (!sameString(a->chrom, b->chrom))
     return 0;
 if (a->strand[0] != b->strand[0])
     return 0;
 int overlap = rangeIntersection(a->chromStart, a->chromEnd, b->chromStart, b->chromEnd);
 if (overlap <= 0)
     return 0;
 
 /* Ok, go do the slow things of building up a range tree.  We could write some code
  * that is a little complex like a merge sort to avoid this, maybe some day if this 
  * is too slow. */
 struct rbTree *t = rangeTreeOnBed(a);
 int totalOverlap = 0;
 int i;
 for (i=0; i<b->blockCount; ++i)
     {
     int start = b->chromStart + b->chromStarts[i];
     int end = start + b->blockSizes[i];
     totalOverlap += rangeTreeOverlapSize(t, start, end);
     }
 rangeTreeFree(&t);
 return totalOverlap;
 }
 
 char *findClosestOfType(struct txGraph *g, struct uniInput *type, struct bed *bed)
 /* Find the accession of the source of given type that overlaps most with bed */
 {
 struct txSource *txSource = g->sources;
 int i;
 int bestOverlap = 0;
 char *bestAcc = NULL;
 for (i=0; i<g->sourceCount; ++i)
     {
     if (sameString(type->mappingBed, txSource->type))
         {
 	struct bed *other = hashMustFindVal(type->bedHash, txSource->accession);
 	int overlap = bedBedOverlap(bed, other);
 	if (overlap > bestOverlap)
 	    {
 	    bestOverlap = overlap;
 	    bestAcc = txSource->accession;
 	    }
 	}
     txSource += 1;
     }
 return bestAcc;
 }
 
 struct txSource *findBestSource(struct txGraph *g, struct bed *bestBed, struct uniInput *uni)
 /* Go through all transcripts belonging to txSource in g, and try to find for them
  * the best mate among the chosen type in the very same g.  This will output no
  * more than one mate for each element of bestType and of ourType */
 {
 char *ourType = uni->mappingBed;
 int i;
 struct txSource *bestSource = NULL;
 int bestOverlap = 0;
 struct txSource *txSource = g->sources;
 for (i=0; i<g->sourceCount; ++i)
     {
     if (sameString(txSource->type, ourType))
         {
 	char *ourAcc = txSource->accession;
 	struct bed *ourBed = hashMustFindVal(uni->bedHash, ourAcc);
 	int overlap = bedBedOverlap(bestBed, ourBed);
 	if (overlap > bestOverlap)
 	    {
 	    bestOverlap = overlap;
 	    bestSource = txSource;
 	    }
 	}
     txSource += 1;
     }
 return bestSource;
 }
 
 
 void outputTrackDb(char *fileName, struct uniInput *uniList)
 /* Write out a trackDb stanza for us */
 {
 FILE *f = mustOpen(fileName, "w");
 fprintf(f, "track uniBarChart\n");
 fprintf(f, "type bigBarChart\n");
 struct dyString *bars = dyStringNew(0), *colors = dyStringNew(0);
 dyStringAppend(bars, "barChartBars");
 dyStringAppend(colors, "barChartColors");
 struct uniInput *uni;
 for (uni = uniList; uni != NULL; uni = uni->next)
     {
     struct lineFile *lf = lineFileOpen(uni->colorFile, TRUE);
     char *line;
     while (lineFileNextReal(lf, &line))
         {
 	char *firstWord = nextTabWord(&line);
 	char *lastWord = NULL;
 	char *word;
 	while ((word = nextTabWord(&line)) != NULL)
 	     lastWord = word;
 	subChar(firstWord, ' ', '_');	// Oh no, underbars!
 	dyStringPrintf(bars, " %s_%s", uni->name, firstWord);
 	dyStringPrintf(colors, " %s", lastWord);
 	}
     lineFileClose(&lf);
     }
 fprintf(f, "%s\n", bars->string);
 fprintf(f, "%s\n", colors->string);
 carefulClose(&f);
 }
 
 void outputStats(char *fileName, struct uniInput *uniList)
 /* Write out a trackDb stanza for us */
 {
 FILE *f = mustOpen(fileName, "w");
 struct uniInput *uni;
 for (uni = uniList; uni != NULL; uni = uni->next)
     {
     struct lineFile *lf = lineFileOpen(uni->statsFile, TRUE);
     char *line;
     while (lineFileNext(lf, &line, NULL))
         {
 	char *start = skipLeadingSpaces(line);
 	if (start[0] == '#' || start[0] == 0)
 	    fprintf(f, "%s\n", line);
 	else
 	    {
 	    fprintf(f, "%s %s\n", uni->name, line);
 	    }
 	}
     lineFileClose(&lf);
     }
 carefulClose(&f);
 }
 
 
 
 void writeOneGraph(struct txGraph *g,  struct uniInput *uniList, struct hash *uniHash, 
     FILE *f, FILE *bedF)
 /* Here we write out one graph.  Normally just write in one piece, but
  * we will at least snoop for those that need breaking up */
 {
 verbose(3, "writeOneGraph %s with %d sources\n", g->name, g->sourceCount);
 /* Find best txSource and it's associated type. */
 int i;
 struct txSource *bestTxSource = NULL;
 int bestPriority = BIGNUM;  // priority is better at 1 than 2
 for (i=0; i<g->sourceCount; ++i)
     {
     struct txSource *txSource = &g->sources[i];
     struct uniInput *uni = hashMustFindVal(uniHash, txSource->type);
     if (bestPriority > uni->priority)
         {
 	bestPriority = uni->priority;
 	bestTxSource = txSource;
 	}
     txSource += 1;
     }
 char *bestType = bestTxSource->type;
 struct uniInput *bestUni = hashMustFindVal(uniHash, bestType);
 struct hash *rowHash = bestUni->rowHash;
 fprintf(f, "%s", bestTxSource->accession);
 
 /* Now go through and output all from that type */
 for (i=0; i<g->sourceCount; ++i)
     {
     struct txSource *txSource = &g->sources[i];
     if (sameString(txSource->type, bestType))
         {
 	if (bedF != NULL)
 	    {
 	    char **row = hashMustFindVal(rowHash, txSource->accession);
 	    writeTsvRow(bedF, BED12_NAME2_SIZE, row);
 	    }
 
 	/* Output matrix */
 	struct uniInput *uni;
 	struct bed *bestBed = hashMustFindVal(bestUni->bedHash, txSource->accession);
 	for (uni = uniList; uni != NULL; uni = uni->next)
 	    {
 	    struct txSource *ourSource = findBestSource(g, bestBed, uni);
 	    if (ourSource != NULL)
 	        {
 		char *geneName = ourSource->accession;
 		struct hashEl *hel = hashLookup(uni->matrixHash, geneName);
 		if (hel == NULL)
 		    errAbort("Can't find %s in %s but it is in %s\n", geneName, uni->matrixFile, ourSource->type);
 		int y = ptToInt(hel->val);
 		int x;
 		for (x=0; x<uni->matrix->xSize; ++x)
 		    {
 		    fprintf(f, "\t%g", uni->scaleOut*uni->matrix->rows[y][x]);
 		    }
 		}
 	    else
 	        {
 		int x;
 		for (x=0; x<uni->matrix->xSize; ++x)
 		    {
 		    fprintf(f, "\t%g", emptyVal);
 		    }
 		}
 	    }
 	fprintf(f, "\n");
 	}
     }
 }
 
 void hcaUnifyMatrix(char *graphInput, char *inList, char *outMatrix, char *outBed)
 /* hcaUnifyMatrix - Given a list of matrices on maybe different gene sets produce a 
  * unified matrix on one gene set. */
 {
 struct txGraph *g, *gList = txGraphLoadAll(graphInput);
 verbose(1, "Read %d graphs from %s\n", slCount(gList), graphInput);
 
 int typeCount = countDistinctTypes(gList);
 int minSources = typeCount/2;
 verbose(1, "%d distinct sources of evidence in %s\n", typeCount, graphInput);
 
 
 struct uniInput *uni, *uniList = uniInputLoadAll(inList);
 struct hash *uniHash = hashNew(0);
 verbose(1, "Read %d inputs from %s\n", slCount(uniList), inList);
 int totalColumns = 0;
 for (uni = uniList; uni != NULL; uni = uni->next)
     {
     hashAdd(uniHash, uni->mappingBed, uni);
     loadMappingBed(uni);
     verbose(2, "%d genes in %s\n", uni->bedHash->elCount, uni->mappingBed);
     loadMatrix(uni);
     verbose(2, "%d genes and %d samples in %s\n", 
 	uni->matrix->ySize, uni->matrix->xSize, uni->matrixFile);
     totalColumns += uni->matrix->xSize;
     }
 verbose(1, "Output will be %d columns\n", totalColumns);
 
 /* Open output and write out header row */
 FILE *f = mustOpen(outMatrix, "w");
 fprintf(f, "#hcaUnifyMatrix");
 for (uni = uniList; uni != NULL; uni = uni->next)
     {
     struct memMatrix *mm = uni->matrix;
     int x;
     for (x=0; x<mm->xSize; ++x)
         {
 	fprintf(f, "\t%s %s", uni->name, mm->xLabels[x]);
 	}
     }
 fprintf(f, "\n");
 
 FILE *bedF = NULL;
 if (outBed != NULL)
     bedF = mustOpen(outBed, "w");
 
 /* Go through outputting matrix and bed rows */
 for (g = gList; g != NULL; g = g->next)
     {
     /* We only want ones with at least a certain minimum of evidence */
     if (g->sourceCount >= minSources)
 	{
 	int typeCount = countDistinctSourceTypes(g);
 	if (typeCount >= minSources)
 	   {
 	   writeOneGraph(g, uniList, uniHash, f, bedF);
 	   }
 	}
     }
 carefulClose(&f);
 carefulClose(&bedF);
 
 char *trackDb = optionVal("trackDb", NULL);
 if (trackDb != NULL)
     {
     outputTrackDb(trackDb, uniList);
     }
 
 char *stats = optionVal("stats", NULL);
 if (stats != NULL)
     {
     outputStats(stats, uniList);
     }
 }
 
 int main(int argc, char *argv[])
 /* Process command line. */
 {
 optionInit(&argc, argv, options);
 if (argc != 4)
     usage();
 emptyVal = optionDouble("empty", -1);
 hcaUnifyMatrix(argv[1], argv[2], argv[3], optionVal("bed", NULL));
 return 0;
 }