3fd32cf592eeadb726874bc41b754985eb099f3e
kent
  Wed Aug 10 21:28:07 2011 -0700
Adding bedOut option to make output bed that is echo of input bed but with mean column appended
diff --git src/utils/bigWigAverageOverBed/bigWigAverageOverBed.c src/utils/bigWigAverageOverBed/bigWigAverageOverBed.c
index 2097e98..0a1e820 100644
--- src/utils/bigWigAverageOverBed/bigWigAverageOverBed.c
+++ src/utils/bigWigAverageOverBed/bigWigAverageOverBed.c
@@ -1,47 +1,52 @@
 /* bigWigAverageOverBed - Compute average score of big wig over each bed, which may have introns. */
 
 #include "common.h"
 #include "linefile.h"
 #include "hash.h"
 #include "localmem.h"
 #include "options.h"
 #include "verbose.h"
 #include "basicBed.h"
 #include "bigWig.h"
 #include "bits.h"
 
 static char const rcsid[] = "$Id: newProg.c,v 1.30 2010/03/24 21:18:33 hiram Exp $";
 
+char *bedOut = NULL;
+
 void usage()
 /* Explain usage and exit. */
 {
 errAbort(
   "bigWigAverageOverBed - Compute average score of big wig over each bed, which may have introns.\n"
   "usage:\n"
   "   bigWigAverageOverBed in.bw in.bed out.tab\n"
   "The output columns are:\n"
   "   name - name field from bed, which should be unique\n"
   "   size - size of bed (sum of exon sizes\n"
   "   covered - # bases within exons covered by bigWig\n"
   "   sum - sum of values over all bases covered\n"
   "   mean0 - average over bases with non-covered bases counting as zeroes\n"
   "   mean - average over just covered bases\n"
+  "Options:\n"
+  "   -bedOut=out.bed - Make output bed that is echo of input bed but with mean column appended\n"
   );
 }
 
 static struct optionSpec options[] = {
+   {"bedOut", OPTION_STRING},
    {NULL, 0},
 };
 
 void checkUniqueNames(struct bed *bedList)
 /* Make sure all names in bedList are unique */
 {
 struct hash *hash = hashNew(16);
 struct bed *bed;
 for (bed = bedList; bed != NULL; bed = bed->next)
     {
     char *name = bed->name;
     if (hashLookup(hash, name) != NULL)
         errAbort("%s duplicated in input bed", name);
     else
         hashAdd(hash, name, NULL);
@@ -66,60 +71,72 @@
     }
 }
 	
 int countBlocks(struct bed *bedList, int fieldCount)
 /* Return the number of blocks in list, or if non-blocked beds, just number of beds. */
 {
 if (fieldCount < 12)
     return slCount(bedList);
 int blockCount = 0;
 struct bed *bed;
 for (bed = bedList; bed != NULL; bed = bed->next)
     blockCount += bed->blockCount;
 return blockCount;
 }
 
-void averageFetchingEachBlock(struct bbiFile *bbi, struct bed *bedList, int fieldCount, FILE *f)
+void optionallyPrintBedPlus(FILE *f, struct bed *bed, int fieldCount, double extra)
+/* Print BED to tab separated file plus an extra double-format column. */
+{
+if (f != NULL)
+    {
+    bedOutputN(bed, fieldCount, f, '\t', '\t');
+    fprintf(f, "%g\n", extra);
+    }
+}
+
+void averageFetchingEachBlock(struct bbiFile *bbi, struct bed *bedList, int fieldCount, 
+	FILE *f, FILE *bedF)
 /* Do the averaging fetching each block from bedList from bigWig.  Fastest for short bedList. */
 {
 struct lm *lm = lmInit(0);
 struct bed *bed;
 for (bed = bedList; bed != NULL; bed = bed->next)
     {
     int coverage = 0;
     double sum = 0.0;
     int size = 0;
 
     if (fieldCount < 12)
 	addBigWigIntervalInfo(bbi, lm, bed->chrom, bed->chromStart, bed->chromEnd, 
 		&size, &coverage, &sum);
     else
 	{
 	int i;
 	for (i=0; i<bed->blockCount; ++i)
 	    {
 	    int start = bed->chromStart + bed->chromStarts[i];
 	    int end = start + bed->blockSizes[i];
 	    addBigWigIntervalInfo(bbi, lm, bed->chrom, start, end, &size, &coverage, &sum);
 	    }
 	}
 
     /* Print out result, fudging mean to 0 if no coverage at all. */
     double mean = 0;
     if (coverage > 0)
 	 mean = sum/coverage;
     fprintf(f, "%s\t%d\t%d\t%g\t%g\t%g\n", bed->name, size, coverage, sum, sum/size, mean);
+    optionallyPrintBedPlus(bedF, bed, fieldCount, mean);
     }
 }
 
 int bedCmpChrom(const void *va, const void *vb)
 /* Compare strings such as chromosome names that may have embedded numbers,
  * so that chr4 comes before chr14 */
 {
 const struct bed *a = *((struct bed **)va);
 const struct bed *b = *((struct bed **)vb);
 return cmpStringsWithEmbeddedNumbers(a->chrom, b->chrom);
 }
 
 struct bed *nextChromInList(struct bed *bedList)
 /* Return first bed in list that starts with another chromosome, or NULL if none. */
 {
@@ -134,31 +151,32 @@
 void addBufIntervalInfo(double *valBuf, Bits *covBuf, int start, int end,
     int *pSumSize, int *pSumCoverage, double *pSumVal)
 /* Look at interval in buffers and add result to sums. */
 {
 int size1 = end - start;
 *pSumSize += size1;
 int cov1 = bitCountRange(covBuf, start, size1);
 *pSumCoverage += cov1;
 int i;
 double sum1 = 0;
 for (i=start; i<end; ++i)
     sum1 += valBuf[i];
 *pSumVal += sum1;
 }
 
-void averageFetchingEachChrom(struct bbiFile *bbi, struct bed **pBedList, int fieldCount, FILE *f)
+void averageFetchingEachChrom(struct bbiFile *bbi, struct bed **pBedList, int fieldCount, 
+	FILE *f, FILE *bedF)
 /* Do the averaging by sorting bedList by chromosome, and then processing each chromosome
  * at once. Faster for long bedLists. */
 {
 /* Sort by chromosome. */
 slSort(pBedList, bedCmpChrom);
 
 struct bigWigValsOnChrom *chromVals = bigWigValsOnChromNew();
 
 struct bed *bed, *bedList, *nextChrom;
 verbose(1, "processing chromosomes");
 for (bedList = *pBedList; bedList != NULL; bedList = nextChrom)
     {
     /* Figure out which chromosome we're working on, and the last bed using it. */
     char *chrom = bedList->chrom;
     nextChrom = nextChromInList(bedList);
@@ -183,69 +201,78 @@
 		{
 		int i;
 		for (i=0; i<bed->blockCount; ++i)
 		    {
 		    int start = bed->chromStart + bed->chromStarts[i];
 		    int end = start + bed->blockSizes[i];
 		    addBufIntervalInfo(valBuf, covBuf, start, end, &size, &coverage, &sum);
 		    }
 		}
 
 	    /* Print out result, fudging mean to 0 if no coverage at all. */
 	    double mean = 0;
 	    if (coverage > 0)
 		 mean = sum/coverage;
 	    fprintf(f, "%s\t%d\t%d\t%g\t%g\t%g\n", bed->name, size, coverage, sum, sum/size, mean);
+	    optionallyPrintBedPlus(bedF, bed, fieldCount, mean);
 	    }
 	verboseDot();
 	}
     else
         {
 	/* If no bigWig data on this chromosome, just output as if coverage is 0 */
 	for (bed = bedList; bed != nextChrom; bed = bed->next)
+	    {
 	    fprintf(f, "%s\t%d\t0\t0\t0\t0\n", bed->name, bedTotalBlockSize(bed));
+	    optionallyPrintBedPlus(bedF, bed, fieldCount, 0);
+	    }
 	}
     }
 verbose(1, "\n");
 }
 
 void bigWigAverageOverBed(char *inBw, char *inBed, char *outTab)
 /* bigWigAverageOverBed - Compute average score of big wig over each bed, which may have introns. */
 {
 struct bed *bedList;
 int fieldCount;
 bedLoadAllReturnFieldCount(inBed, &bedList, &fieldCount);
 checkUniqueNames(bedList);
 
 struct bbiFile *bbi = bigWigFileOpen(inBw);
 FILE *f = mustOpen(outTab, "w");
+FILE *bedF = NULL;
+if (bedOut != NULL)
+    bedF = mustOpen(bedOut, "w");
 
 /* Count up number of blocks in file.  It takes about 1/100th of of second to
  * look up a single block in a bigWig.  On the other hand to stream through
  * the whole file setting a array of doubles takes about 30 seconds, so we change
  * strategy at 3,000 blocks. 
  *   I (Jim) usually avoid having two paths through the code like this, and am tempted
  * to always go the ~30 second chromosome-at-a-time  way.  On the other hand the block-way
  * was developed first, and it was useful to have both ways to test against each other.
  * (This found a bug where the chromosome way wasn't handling beds in chromosomes not
  * covered by the bigWig for instance).  Since this code is not likely to change too
  * much, keeping both implementations in seems reasonable. */
 int blockCount = countBlocks(bedList, fieldCount);
 verbose(2, "Got %d blocks, if >= 3000 will use chromosome-at-a-time method\n", blockCount);
 
 if (blockCount < 3000)
-    averageFetchingEachBlock(bbi, bedList, fieldCount, f);
+    averageFetchingEachBlock(bbi, bedList, fieldCount, f, bedF);
 else
-    averageFetchingEachChrom(bbi, &bedList, fieldCount, f);
+    averageFetchingEachChrom(bbi, &bedList, fieldCount, f, bedF);
 
+carefulClose(&bedF);
 carefulClose(&f);
 }
 
 int main(int argc, char *argv[])
 /* Process command line. */
 {
 optionInit(&argc, argv, options);
 if (argc != 4)
     usage();
+bedOut = optionVal("bedOut", bedOut);
 bigWigAverageOverBed(argv[1], argv[2], argv[3]);
 return 0;
 }