@@ -1,11 +1,13 @@

 /* bwgCreate - create big wig files.  Implements write side of bwgInternal.h module. 

  * See the comment in bwgInternal.h for a description of the file format. */

+/* Copyright (C) 2014 The Regents of the University of California 

+ * See README in this or parent directory for licensing information. */

+

 #include "common.h"

 #include "linefile.h"

 #include "hash.h"

-#include "localmem.h"

-#include "errabort.h"

+#include "errAbort.h"

 #include "sqlNum.h"

 #include "sig.h"

 #include "zlibFace.h"

@@ -15,7 +17,6 @@

 #include "bwgInternal.h"

 #include "bigWig.h"

-static char const rcsid[] = "$Id: bwgCreate.c,v 1.27 2010/06/10 20:13:29 braney Exp $";

 static int bwgBedGraphItemCmp(const void *va, const void *vb)

 /* Compare to sort based on query start. */

@@ -628,6 +629,78 @@ slFreeList(&uniqList);

 *retMaxChromNameSize = maxChromNameSize;

 }

+static int bwgStrcmp (const void * A, const void * B) {

+	char * stringA = *((char **) A);

+	char * stringB = *((char **) B);

+	return strcmp(stringA, stringB);

+}

+

+void bwgMakeAllChromInfo(struct bwgSection *sectionList, struct hash *chromSizeHash,

+	int *retChromCount, struct bbiChromInfo **retChromArray,

+	int *retMaxChromNameSize)

+/* Fill in chromId field in sectionList.  Return array of chromosome name/ids. 

+ * The chromSizeHash is keyed by name, and has int values. */

+{

+/* Build up list of unique chromosome names. */

+int maxChromNameSize = 0;

+

+/* Get list of values */

+int chromCount = chromSizeHash->elCount;

+char ** chromName, ** chromNames;

+AllocArray(chromNames, chromCount);

+chromName = chromNames;

+struct hashEl* el;

+struct hashCookie cookie = hashFirst(chromSizeHash);

+for (el = hashNext(&cookie); el; el = hashNext(&cookie)) {

+	*chromName = el->name;

+	if (strlen(el->name) > maxChromNameSize)

+		maxChromNameSize = strlen(el->name);

+	chromName++;

+}

+qsort(chromNames, chromCount, sizeof(char *), bwgStrcmp);

+

+/* Allocate and fill in results array. */

+struct bbiChromInfo *chromArray;

+AllocArray(chromArray, chromCount);

+int i;

+for (i = 0; i < chromCount; ++i)

+    {

+    chromArray[i].name = chromNames[i];

+    chromArray[i].id = i;

+    chromArray[i].size = hashIntVal(chromSizeHash, chromNames[i]);

+    }

+

+// Assign IDs to sections:

+struct bwgSection *section;

+char *name = "";

+bits32 chromId = 0;

+for (section = sectionList; section != NULL; section = section->next)

+    {

+    if (!sameString(section->chrom, name))

+        {

+        for (i = 0; i < chromCount; ++i)

+            {

+	    if (sameString(section->chrom, chromArray[i].name)) 

+	        {

+		    section->chromId = i;

+	    	    break;

+	        }

+	    }

+	if (i == chromCount)

+		errAbort("Could not find %s in list of chromosomes\n", section->chrom);

+	chromId = section->chromId;

+	name = section->chrom;

+	}

+    else 

+	section->chromId = chromId;

+    }

+

+/* Clean up, set return values and go home. */

+*retChromCount = chromCount;

+*retChromArray = chromArray;

+*retMaxChromNameSize = maxChromNameSize;

+}

+

 int bwgAverageResolution(struct bwgSection *sectionList)

 /* Return the average resolution seen in sectionList. */

 {

@@ -798,46 +871,18 @@ slReverse(&outList);

 return outList;

 }

-void bwgCreate(struct bwgSection *sectionList, struct hash *chromSizeHash, 

-	int blockSize, int itemsPerSlot, boolean doCompress, char *fileName)

-/* Create a bigWig file out of a sorted sectionList. */

-{

-bits64 sectionCount = slCount(sectionList);

-FILE *f = mustOpen(fileName, "wb");

-bits32 sig = bigWigSig;

-bits16 version = bbiCurrentVersion;

-bits16 summaryCount = 0;

-bits16 reserved16 = 0;

-bits32 reserved32 = 0;

-bits64 reserved64 = 0;

-bits64 dataOffset = 0, dataOffsetPos;

-bits64 indexOffset = 0, indexOffsetPos;

-bits64 chromTreeOffset = 0, chromTreeOffsetPos;

-bits64 totalSummaryOffset = 0, totalSummaryOffsetPos;

-bits32 uncompressBufSize = 0;

-bits64 uncompressBufSizePos;

-struct bbiSummary *reduceSummaries[10];

-bits32 reductionAmounts[10];

-bits64 reductionDataOffsetPos[10];

-bits64 reductionDataOffsets[10];

-bits64 reductionIndexOffsets[10];

-int i;

-

-/* Figure out chromosome ID's. */

-struct bbiChromInfo *chromInfoArray;

-int chromCount, maxChromNameSize;

-bwgMakeChromInfo(sectionList, chromSizeHash, &chromCount, &chromInfoArray, &maxChromNameSize);

-

+static void bwgComputeDynamicSummaries(struct bwgSection *sectionList, struct bbiSummary ** reduceSummaries, bits16 * summaryCount, struct bbiChromInfo *chromInfoArray, int chromCount, bits32 * reductionAmounts, boolean doCompress) {

 /* Figure out initial summary level - starting with a summary 10 times the amount

  * of the smallest item.  See if summarized data is smaller than half input data, if

  * not bump up reduction by a factor of 2 until it is, or until further summarying

  * yeilds no size reduction. */

+int i;

 int  minRes = bwgAverageResolution(sectionList);

 int initialReduction = minRes*10;

 bits64 fullSize = bwgTotalSectionSize(sectionList);

-bits64 maxReducedSize = fullSize/2;

-struct bbiSummary *firstSummaryList = NULL, *summaryList = NULL;

 bits64 lastSummarySize = 0, summarySize;

+bits64 maxReducedSize = fullSize/2;

+struct bbiSummary *summaryList = NULL;

 for (;;)

     {

     summaryList = bwgReduceSectionList(sectionList, chromInfoArray, initialReduction);

@@ -862,31 +907,92 @@ for (;;)

     else

         break;

     }

-summaryCount = 1;

-reduceSummaries[0] = firstSummaryList = summaryList;

+*summaryCount = 1;

+reduceSummaries[0] = summaryList;

 reductionAmounts[0] = initialReduction;

 /* Now calculate up to 10 levels of further summary. */

 bits64 reduction = initialReduction;

-for (i=0; i<ArraySize(reduceSummaries)-1; i++)

+for (i=0; i<9; i++)

     {

     reduction *= 4;

     if (reduction > 1000000000)

         break;

-    summaryList = bbiReduceSummaryList(reduceSummaries[summaryCount-1], chromInfoArray, 

+    summaryList = bbiReduceSummaryList(reduceSummaries[*summaryCount-1], chromInfoArray, 

     	reduction);

     summarySize = bbiTotalSummarySize(summaryList);

     if (summarySize != lastSummarySize)

         {

- 	reduceSummaries[summaryCount] = summaryList;

-	reductionAmounts[summaryCount] = reduction;

-	++summaryCount;

+ 	reduceSummaries[*summaryCount] = summaryList;

+	reductionAmounts[*summaryCount] = reduction;

+	++(*summaryCount);

 	}

     int summaryItemCount = slCount(summaryList);

     if (summaryItemCount <= chromCount)

         break;

     }

+}

+

+static void bwgComputeFixedSummaries(struct bwgSection * sectionList, struct bbiSummary ** reduceSummaries, bits16 * summaryCount, struct bbiChromInfo *chromInfoArray, bits32 * reductionAmounts) {

+// Hack: pre-defining summary levels, set off Ensembl default zoom levels

+// The last two values of this array were extrapolated following Jim's formula

+int i;

+#define REDUCTION_COUNT 10

+bits32 presetReductions[REDUCTION_COUNT] = {30, 65, 130, 260, 450, 648, 950, 1296, 4800, 19200}; 

+

+bits64 reduction = reductionAmounts[0] = presetReductions[0];

+reduceSummaries[0] = bwgReduceSectionList(sectionList, chromInfoArray, presetReductions[0]);

+

+for (i=1; i<REDUCTION_COUNT; i++)

+    {

+    reduction = reductionAmounts[i] = presetReductions[i];

+    reduceSummaries[i] = bbiReduceSummaryList(reduceSummaries[i-1], chromInfoArray, 

+    	reduction);

+    }

+

+*summaryCount = REDUCTION_COUNT;

+}

+

+void bwgCreate(struct bwgSection *sectionList, struct hash *chromSizeHash, 

+	int blockSize, int itemsPerSlot, boolean doCompress, boolean keepAllChromosomes,

+        boolean fixedSummaries, char *fileName)

+/* Create a bigWig file out of a sorted sectionList. */

+{

+bits64 sectionCount = slCount(sectionList);

+FILE *f = mustOpen(fileName, "wb");

+bits32 sig = bigWigSig;

+bits16 version = bbiCurrentVersion;

+bits16 summaryCount = 0;

+bits16 reserved16 = 0;

+bits32 reserved32 = 0;

+bits64 reserved64 = 0;

+bits64 dataOffset = 0, dataOffsetPos;

+bits64 indexOffset = 0, indexOffsetPos;

+bits64 chromTreeOffset = 0, chromTreeOffsetPos;

+bits64 totalSummaryOffset = 0, totalSummaryOffsetPos;

+bits32 uncompressBufSize = 0;

+bits64 uncompressBufSizePos;

+struct bbiSummary *reduceSummaries[10];

+bits32 reductionAmounts[10];

+bits64 reductionDataOffsetPos[10];

+bits64 reductionDataOffsets[10];

+bits64 reductionIndexOffsets[10];

+int i;

+

+/* Figure out chromosome ID's. */

+struct bbiChromInfo *chromInfoArray;

+int chromCount, maxChromNameSize;

+if (keepAllChromosomes)

+    bwgMakeAllChromInfo(sectionList, chromSizeHash, &chromCount, &chromInfoArray, &maxChromNameSize);

+else

+    bwgMakeChromInfo(sectionList, chromSizeHash, &chromCount, &chromInfoArray, &maxChromNameSize);

+

+if (fixedSummaries) 

+    bwgComputeFixedSummaries(sectionList, reduceSummaries, &summaryCount, chromInfoArray, reductionAmounts);

+else 

+    bwgComputeDynamicSummaries(sectionList, reduceSummaries, &summaryCount, chromInfoArray, chromCount, reductionAmounts, doCompress);

+

 /* Write fixed header. */

 writeOne(f, sig);

 writeOne(f, version);

@@ -904,8 +1010,8 @@ totalSummaryOffsetPos = ftell(f);

 writeOne(f, totalSummaryOffset);

 uncompressBufSizePos = ftell(f);

 writeOne(f, uncompressBufSize);

-for (i=0; i<2; ++i)

-    writeOne(f, reserved32);

+writeOne(f, reserved64);  /* nameIndexOffset */

+assert(ftell(f) == 64);

 /* Write summary headers */

 for (i=0; i<summaryCount; ++i)

@@ -964,7 +1070,7 @@ for (i=0; i<summaryCount; ++i)

     }

 /* Calculate summary */

-struct bbiSummary *sum = firstSummaryList;

+struct bbiSummary *sum = reduceSummaries[0];

 if (sum != NULL)

     {

     totalSum.validCount = sum->validCount;

@@ -1087,13 +1193,15 @@ for (section = sectionList; section != NULL; section = nextSection)

 return sectionList;

 }

-void bigWigFileCreate(

+void bigWigFileCreateEx(

 	char *inName, 		/* Input file in ascii wiggle format. */

 	char *chromSizes, 	/* Two column tab-separated file: <chromosome> <size>. */

 	int blockSize,		/* Number of items to bundle in r-tree.  1024 is good. */

 	int itemsPerSlot,	/* Number of items in lowest level of tree.  512 is good. */

 	boolean clipDontDie,	/* If TRUE then clip items off end of chrom rather than dying. */

 	boolean compress,	/* If TRUE then compress data. */

+	boolean keepAllChromosomes,	/* If TRUE then store all chromosomes in chromosomal b-tree. */

+	boolean fixedSummaries,	/* If TRUE then impose fixed summary levels. */

 	char *outName)

 /* Convert ascii format wig file (in fixedStep, variableStep or bedGraph format) 

  * to binary big wig format. */

@@ -1109,7 +1217,22 @@ struct lm *lm = lmInit(0);

 struct bwgSection *sectionList = bwgParseWig(inName, clipDontDie, chromSizeHash, itemsPerSlot, lm);

 if (sectionList == NULL)

     errAbort("%s is empty of data", inName);

-bwgCreate(sectionList, chromSizeHash, blockSize, itemsPerSlot, compress, outName);

+bwgCreate(sectionList, chromSizeHash, blockSize, itemsPerSlot, compress, keepAllChromosomes, fixedSummaries, outName);

 lmCleanup(&lm);

 }

+void bigWigFileCreate(

+	char *inName, 		/* Input file in ascii wiggle format. */

+	char *chromSizes, 	/* Two column tab-separated file: <chromosome> <size>. */

+	int blockSize,		/* Number of items to bundle in r-tree.  1024 is good. */

+	int itemsPerSlot,	/* Number of items in lowest level of tree.  512 is good. */

+	boolean clipDontDie,	/* If TRUE then clip items off end of chrom rather than dying. */

+	boolean compress,	/* If TRUE then compress data. */

+	char *outName)

+/* Convert ascii format wig file (in fixedStep, variableStep or bedGraph format) 

+ * to binary big wig format. */

+{

+bigWigFileCreateEx( inName, chromSizes, blockSize, itemsPerSlot, clipDontDie,

+	compress, FALSE, FALSE, outName);

+}

+

@@ -15,70 +15,7 @@

 #include "bwgInternal.h"

 #include "bigWig.h"

-static char const rcsid[] = "$Id: bwgCreate.c,v 1.22 2009/11/25 07:17:25 kent Exp $";

-

-struct bwgBedGraphItem

-/* An bedGraph-type item in a bwgSection. */

-    {

-    struct bwgBedGraphItem *next;	/* Next in list. */

-    bits32 start,end;		/* Range of chromosome covered. */

-    float val;			/* Value. */

-    };

-

-struct bwgVariableStepItem

-/* An variableStep type item in a bwgSection. */

-    {

-    struct bwgVariableStepItem *next;	/* Next in list. */

-    bits32 start;		/* Start position in chromosome. */

-    float val;			/* Value. */

-    };

-

-struct bwgVariableStepPacked

-/* An variableStep type item in a bwgSection. */

-    {

-    bits32 start;		/* Start position in chromosome. */

-    float val;			/* Value. */

-    };

-

-struct bwgFixedStepItem

-/* An fixedStep type item in a bwgSection. */

-    {

-    struct bwgFixedStepItem *next;	/* Next in list. */

-    float val;			/* Value. */

-    };

-

-struct bwgFixedStepPacked

-/* An fixedStep type item in a bwgSection. */

-    {

-    float val;			/* Value. */

-    };

-

-union bwgItem

-/* Union of item pointers for all possible section types. */

-    {

-    struct bwgBedGraphItem *bedGraphList;		/* A linked list */

-    struct bwgFixedStepPacked *fixedStepPacked;		/* An array */

-    struct bwgVariableStepPacked *variableStepPacked;	/* An array */

-    /* No packed format for bedGraph... */

-    };

-

-struct bwgSection

-/* A section of a bigWig file - all on same chrom.  This is a somewhat fat data

- * structure used by the bigWig creation code.  See also bwgSection for the

- * structure returned by the bigWig reading code. */

-    {

-    struct bwgSection *next;		/* Next in list. */

-    char *chrom;			/* Chromosome name. */

-    bits32 start,end;			/* Range of chromosome covered. */

-    enum bwgSectionType type;

-    union bwgItem items;		/* List/array of items in this section. */

-    bits32 itemStep;			/* Step within item if applicable. */

-    bits32 itemSpan;			/* Item span if applicable. */

-    bits16 itemCount;			/* Number of items in section. */

-    bits32 chromId;			/* Unique small integer value for chromosome. */

-    bits64 fileOffset;			/* Offset of section in file. */

-    };

-

+static char const rcsid[] = "$Id: bwgCreate.c,v 1.27 2010/06/10 20:13:29 braney Exp $";

 static int bwgBedGraphItemCmp(const void *va, const void *vb)

 /* Compare to sort based on query start. */

@@ -200,7 +137,7 @@ return bufSize;

 }

-static int bwgSectionCmp(const void *va, const void *vb)

+int bwgSectionCmp(const void *va, const void *vb)

 /* Compare to sort based on chrom,start,end.  */

 {

 const struct bwgSection *a = *((struct bwgSection **)va);

@@ -338,7 +275,7 @@ struct lm *lmLocal = lmInit(0);

  * adding values from single column to list. */

 char *words[2];

 char *line;

-struct bwgVariableStepItem *item, *itemList = NULL;

+struct bwgVariableStepItem *item, *nextItem, *itemList = NULL;

 int originalSectionSize = 0;

 while (lineFileNextReal(lf, &line))

     {

@@ -372,6 +309,20 @@ while (lineFileNextReal(lf, &line))

     }

 slSort(&itemList, bwgVariableStepItemCmp);

+/* Make sure no overlap between items. */

+if (itemList != NULL)

+    {

+    item = itemList;

+    for (nextItem = item->next; nextItem != NULL; nextItem = nextItem->next)

+        {

+	if (item->start + span > nextItem->start)

+	    errAbort("Overlap on %s between items starting at %d and %d.\n"

+	             "Please remove overlaps and try again",

+		    chrom, item->start, nextItem->start);

+	item = nextItem;

+	}

+    }

+

 /* Break up into sections of no more than items-per-slot size. */

 int sizeLeft = originalSectionSize;

 for (item = itemList; item != NULL; )

@@ -435,7 +386,7 @@ else

     errAbort("Unknown type %s\n", typeWord);

 /* Set up defaults for values we hope to parse out of rest of line. */

-int span = 1;

+int span = 0;

 bits32 step = 0;

 bits32 start = 0;

 char *chrom = NULL;

@@ -468,8 +419,8 @@ while ((varEqVal = nextWord(&initialLine)) != NULL)

  * rest of section. */

 if (chrom == NULL)

     errAbort("Missing chrom= setting line %d of %s\n", lf->lineIx, lf->fileName);

-bits32 chromSize = hashIntVal(chromSizeHash, chrom);

-if (start >= chromSize)

+bits32 chromSize = (chromSizeHash ? hashIntVal(chromSizeHash, chrom) : BIGNUM);

+if (start > chromSize)

     {

     warn("line %d of %s: chromosome %s has %u bases, but item starts at %u",

     	lf->lineIx, lf->fileName, chrom, chromSize, start);

@@ -482,6 +433,8 @@ if (type == bwgTypeFixedStep)

 	errAbort("Missing start= setting line %d of %s\n", lf->lineIx, lf->fileName);

     if (step == 0)

 	errAbort("Missing step= setting line %d of %s\n", lf->lineIx, lf->fileName);

+    if (span == 0)

+	span = step;

     parseFixedStepSection(lf, clipDontDie, lm, itemsPerSlot, 

     	chrom, chromSize, span, start-1, step, pSectionList);

     }

@@ -491,6 +444,8 @@ else

 	errAbort("Extra start= setting line %d of %s\n", lf->lineIx, lf->fileName);

     if (step != 0)

 	errAbort("Extra step= setting line %d of %s\n", lf->lineIx, lf->fileName);

+    if (span == 0)

+	span = 1;

     parseVariableStepSection(lf, clipDontDie, lm, itemsPerSlot, 

     	chrom, chromSize, span, pSectionList);

     }

@@ -546,7 +501,7 @@ while (lineFileNextReal(lf, &line))

         {

 	lmAllocVar(chromHash->lm, chrom);

 	hashAddSaveName(chromHash, chromName, chrom, &chrom->name);

-	chrom->size = hashIntVal(chromSizeHash, chromName);

+	chrom->size = (chromSizeHash ? hashIntVal(chromSizeHash, chromName) : BIGNUM);

 	slAddHead(&chromList, chrom);

 	}

@@ -574,10 +529,22 @@ while (lineFileNextReal(lf, &line))

     }

 slSort(&chromList, bedGraphChromCmpName);

+/* Loop through each chromosome and output the item list, broken into sections

+ * for that chrom. */

 for (chrom = chromList; chrom != NULL; chrom = chrom->next)

     {

     slSort(&chrom->itemList, bwgBedGraphItemCmp);

+    /* Check to make sure no overlap between items. */

+    struct bwgBedGraphItem *item = chrom->itemList, *nextItem;

+    for (nextItem = item->next; nextItem != NULL; nextItem = nextItem->next)

+        {

+	if (item->end > nextItem->start)

+	    errAbort("Overlap between %s %d %d and %s %d %d.\nPlease remove overlaps and try again",

+	        chrom->name, item->start, item->end, chrom->name, nextItem->start, nextItem->end);

+	item = nextItem;

+	}

+

     /* Break up into sections of no more than items-per-slot size. */

     struct bwgBedGraphItem *startItem, *endItem, *nextStartItem = chrom->itemList;

     for (startItem = chrom->itemList; startItem != NULL; startItem = nextStartItem)

@@ -832,8 +799,7 @@ return outList;

 }

 void bwgCreate(struct bwgSection *sectionList, struct hash *chromSizeHash, 

-               int blockSize, int itemsPerSlot, boolean doCompress,

-               char *fileName)

+	int blockSize, int itemsPerSlot, boolean doCompress, char *fileName)

 /* Create a bigWig file out of a sorted sectionList. */

 {

 bits64 sectionCount = slCount(sectionList);

@@ -862,12 +828,12 @@ struct bbiChromInfo *chromInfoArray;

 int chromCount, maxChromNameSize;

 bwgMakeChromInfo(sectionList, chromSizeHash, &chromCount, &chromInfoArray, &maxChromNameSize);

-/* Figure out initial summary level - starting with a summary 20 times the amount

+/* Figure out initial summary level - starting with a summary 10 times the amount

  * of the smallest item.  See if summarized data is smaller than half input data, if

  * not bump up reduction by a factor of 2 until it is, or until further summarying

  * yeilds no size reduction. */

 int  minRes = bwgAverageResolution(sectionList);

-int initialReduction = minRes*20;

+int initialReduction = minRes*10;

 bits64 fullSize = bwgTotalSectionSize(sectionList);

 bits64 maxReducedSize = fullSize/2;

 struct bbiSummary *firstSummaryList = NULL, *summaryList = NULL;

@@ -878,8 +844,7 @@ for (;;)

     bits64 summarySize = bbiTotalSummarySize(summaryList);

     if (doCompress)

 	{

-        summarySize *= 4;	// Compensate for summary not compressing as well as primary data

-	initialReduction *= 4;

+        summarySize *= 2;	// Compensate for summary not compressing as well as primary data

 	}

     if (summarySize >= maxReducedSize && summarySize != lastSummarySize)

         {

@@ -1049,18 +1014,30 @@ for (i=0; i<summaryCount; ++i)

     writeOne(f, reductionIndexOffsets[i]);

     }

+/* Write end signature. */

+fseek(f, 0L, SEEK_END);

+writeOne(f, sig);

+

 /* Clean up */

 freez(&chromInfoArray);

 carefulClose(&f);

 }

-struct bwgSection *bwgParseWig(char *fileName, boolean clipDontDie, struct hash *chromSizeHash,

-	int maxSectionSize, struct lm *lm)

+struct bwgSection *bwgParseWig(

+	char *fileName,       /* Name of ascii wig file. */

+	boolean clipDontDie,  /* Skip items outside chromosome rather than aborting. */

+	struct hash *chromSizeHash,  /* If non-NULL items checked to be inside chromosome. */

+	int maxSectionSize,   /* Biggest size of a section.  100 - 100,000 is usual range. */

+	struct lm *lm)	      /* Memory pool to allocate from. */

 /* Parse out ascii wig file - allocating memory in lm. */

 {

 struct lineFile *lf = lineFileOpen(fileName, TRUE);

 char *line;

 struct bwgSection *sectionList = NULL;

+

+/* remove initial browser and track lines */

+lineFileRemoveInitialCustomTrackLines(lf);

+

 while (lineFileNextReal(lf, &line))

     {

     verbose(2, "processing %s\n", line);

@@ -1089,7 +1066,7 @@ while (lineFileNextReal(lf, &line))

     }

 slSort(&sectionList, bwgSectionCmp);

-/* Check for overlap. */

+/* Check for overlap at section level. */

 struct bwgSection *section, *nextSection;

 for (section = sectionList; section != NULL; section = nextSection)

     {

new file mode 100644

@@ -0,0 +1,1138 @@

+/* bwgCreate - create big wig files.  Implements write side of bwgInternal.h module. 

+ * See the comment in bwgInternal.h for a description of the file format. */

+

+#include "common.h"

+#include "linefile.h"

+#include "hash.h"

+#include "localmem.h"

+#include "errabort.h"

+#include "sqlNum.h"

+#include "sig.h"

+#include "zlibFace.h"

+#include "bPlusTree.h"

+#include "cirTree.h"

+#include "bbiFile.h"

+#include "bwgInternal.h"

+#include "bigWig.h"

+

+static char const rcsid[] = "$Id: bwgCreate.c,v 1.22 2009/11/25 07:17:25 kent Exp $";

+

+struct bwgBedGraphItem

+/* An bedGraph-type item in a bwgSection. */

+    {

+    struct bwgBedGraphItem *next;	/* Next in list. */

+    bits32 start,end;		/* Range of chromosome covered. */

+    float val;			/* Value. */

+    };

+

+struct bwgVariableStepItem

+/* An variableStep type item in a bwgSection. */

+    {

+    struct bwgVariableStepItem *next;	/* Next in list. */

+    bits32 start;		/* Start position in chromosome. */

+    float val;			/* Value. */

+    };

+

+struct bwgVariableStepPacked

+/* An variableStep type item in a bwgSection. */

+    {

+    bits32 start;		/* Start position in chromosome. */

+    float val;			/* Value. */

+    };

+

+struct bwgFixedStepItem

+/* An fixedStep type item in a bwgSection. */

+    {

+    struct bwgFixedStepItem *next;	/* Next in list. */

+    float val;			/* Value. */

+    };

+

+struct bwgFixedStepPacked

+/* An fixedStep type item in a bwgSection. */

+    {

+    float val;			/* Value. */

+    };

+

+union bwgItem

+/* Union of item pointers for all possible section types. */

+    {

+    struct bwgBedGraphItem *bedGraphList;		/* A linked list */

+    struct bwgFixedStepPacked *fixedStepPacked;		/* An array */

+    struct bwgVariableStepPacked *variableStepPacked;	/* An array */

+    /* No packed format for bedGraph... */

+    };

+

+struct bwgSection

+/* A section of a bigWig file - all on same chrom.  This is a somewhat fat data

+ * structure used by the bigWig creation code.  See also bwgSection for the

+ * structure returned by the bigWig reading code. */

+    {

+    struct bwgSection *next;		/* Next in list. */

+    char *chrom;			/* Chromosome name. */

+    bits32 start,end;			/* Range of chromosome covered. */

+    enum bwgSectionType type;

+    union bwgItem items;		/* List/array of items in this section. */

+    bits32 itemStep;			/* Step within item if applicable. */

+    bits32 itemSpan;			/* Item span if applicable. */

+    bits16 itemCount;			/* Number of items in section. */

+    bits32 chromId;			/* Unique small integer value for chromosome. */

+    bits64 fileOffset;			/* Offset of section in file. */

+    };

+

+

+static int bwgBedGraphItemCmp(const void *va, const void *vb)

+/* Compare to sort based on query start. */

+{

+const struct bwgBedGraphItem *a = *((struct bwgBedGraphItem **)va);

+const struct bwgBedGraphItem *b = *((struct bwgBedGraphItem **)vb);

+int dif = (int)a->start - (int)b->start;

+if (dif == 0)

+    dif = (int)a->end - (int)b->end;

+return dif;

+}

+

+static int bwgVariableStepItemCmp(const void *va, const void *vb)

+/* Compare to sort based on query start. */

+{

+const struct bwgVariableStepItem *a = *((struct bwgVariableStepItem **)va);

+const struct bwgVariableStepItem *b = *((struct bwgVariableStepItem **)vb);

+return (int)a->start - (int)b->start;

+}

+

+void bwgDumpSummary(struct bbiSummary *sum, FILE *f)

+/* Write out summary info to file. */

+{

+fprintf(f, "summary %d:%d-%d min=%f, max=%f, sum=%f, sumSquares=%f, validCount=%d, mean=%f\n",

+     sum->chromId, sum->start, sum->end, sum->minVal, sum->maxVal, sum->sumData,

+     sum->sumSquares, sum->validCount, sum->sumData/sum->validCount);

+}

+

+static int bwgSectionWrite(struct bwgSection *section, boolean doCompress, FILE *f)

+/* Write out section to file, filling in section->fileOffset. */

+{

+UBYTE type = section->type;

+UBYTE reserved8 = 0;

+int itemSize;

+switch (section->type)

+    {

+    case bwgTypeBedGraph:

+        itemSize = 12;

+	break;

+    case bwgTypeVariableStep:

+        itemSize = 8;

+	break;

+    case bwgTypeFixedStep:

+        itemSize = 4;

+	break;

+    default:

+        itemSize = 0;  // Suppress compiler warning

+	internalErr();

+	break;

+    }

+int fixedSize = sizeof(section->chromId) + sizeof(section->start) + sizeof(section->end) + 

+     sizeof(section->itemStep) + sizeof(section->itemSpan) + sizeof(type) + sizeof(reserved8) +

+     sizeof(section->itemCount);

+int bufSize = section->itemCount * itemSize + fixedSize;

+char buf[bufSize];

+char *bufPt = buf;

+

+section->fileOffset = ftell(f);

+memWriteOne(&bufPt, section->chromId);

+memWriteOne(&bufPt, section->start);

+memWriteOne(&bufPt, section->end);

+memWriteOne(&bufPt, section->itemStep);

+memWriteOne(&bufPt, section->itemSpan);

+memWriteOne(&bufPt, type);

+memWriteOne(&bufPt, reserved8);

+memWriteOne(&bufPt, section->itemCount);

+

+int i;

+switch (section->type)

+    {

+    case bwgTypeBedGraph:

+	{

+	struct bwgBedGraphItem *item = section->items.bedGraphList;

+	for (item = section->items.bedGraphList; item != NULL; item = item->next)

+	    {

+	    memWriteOne(&bufPt, item->start);

+	    memWriteOne(&bufPt, item->end);

+	    memWriteOne(&bufPt, item->val);

+	    }

+	break;

+	}

+    case bwgTypeVariableStep:

+	{

+	struct bwgVariableStepPacked *items = section->items.variableStepPacked;

+	for (i=0; i<section->itemCount; ++i)

+	    {

+	    memWriteOne(&bufPt, items->start);

+	    memWriteOne(&bufPt, items->val);

+	    items += 1;

+	    }

+	break;

+	}

+    case bwgTypeFixedStep:

+	{

+	struct bwgFixedStepPacked *items = section->items.fixedStepPacked;

+	for (i=0; i<section->itemCount; ++i)

+	    {

+	    memWriteOne(&bufPt, items->val);

+	    items += 1;

+	    }

+	break;

+	}

+    default:

+        internalErr();

+	break;

+    }

+assert(bufSize == (bufPt - buf) );

+

+if (doCompress)

+    {

+    size_t maxCompSize = zCompBufSize(bufSize);

+    char compBuf[maxCompSize];

+    int compSize = zCompress(buf, bufSize, compBuf, maxCompSize);

+    mustWrite(f, compBuf, compSize);

+    }

+else

+    mustWrite(f, buf, bufSize);

+return bufSize;

+}

+

+

+static int bwgSectionCmp(const void *va, const void *vb)

+/* Compare to sort based on chrom,start,end.  */

+{

+const struct bwgSection *a = *((struct bwgSection **)va);

+const struct bwgSection *b = *((struct bwgSection **)vb);

+int dif = strcmp(a->chrom, b->chrom);

+if (dif == 0)

+    {

+    dif = (int)a->start - (int)b->start;

+    if (dif == 0)

+	dif = (int)a->end - (int)b->end;

+    }

+return dif;

+}

+

+static struct cirTreeRange bwgSectionFetchKey(const void *va, void *context)

+/* Fetch bwgSection key for r-tree */

+{

+struct cirTreeRange res;

+const struct bwgSection *a = *((struct bwgSection **)va);

+res.chromIx = a->chromId;

+res.start = a->start;

+res.end = a->end;

+return res;

+}

+

+static bits64 bwgSectionFetchOffset(const void *va, void *context)

+/* Fetch bwgSection file offset for r-tree */

+{

+const struct bwgSection *a = *((struct bwgSection **)va);

+return a->fileOffset;

+}

+

+static boolean stepTypeLine(char *line)

+/* Return TRUE if it's a variableStep or fixedStep line. */

+{

+return (startsWithWord("variableStep", line) || startsWithWord("fixedStep", line));

+}

+

+static boolean steppedSectionEnd(char *line, int maxWords)

+/* Return TRUE if line indicates the start of another section. */

+{

+int wordCount = chopByWhite(line, NULL, 5);

+if (wordCount > maxWords)

+    return TRUE;

+return stepTypeLine(line);

+}

+

+static void parseFixedStepSection(struct lineFile *lf, boolean clipDontDie, struct lm *lm,

+	int itemsPerSlot, char *chrom, bits32 chromSize, bits32 span, bits32 sectionStart, 

+	bits32 step, struct bwgSection **pSectionList)

+/* Read the single column data in section until get to end. */

+{

+struct lm *lmLocal = lmInit(0);

+

+/* Stream through section until get to end of file or next section,

+ * adding values from single column to list. */

+char *words[1];

+char *line;

+struct bwgFixedStepItem *item, *itemList = NULL;

+int originalSectionSize = 0;

+bits32 sectionEnd = sectionStart;

+while (lineFileNextReal(lf, &line))

+    {

+    if (steppedSectionEnd(line, 1))

+	{

+        lineFileReuse(lf);

+	break;

+	}

+    chopLine(line, words);

+    lmAllocVar(lmLocal, item);

+    item->val = lineFileNeedDouble(lf, words, 0);

+    if (sectionEnd + span > chromSize)

+	{

+	warn("line %d of %s: chromosome %s has %u bases, but item ends at %u",

+	    lf->lineIx, lf->fileName, chrom, chromSize, sectionEnd + span);

+	if (!clipDontDie)

+	    noWarnAbort();

+	}

+    else

+	{

+	slAddHead(&itemList, item);

+	++originalSectionSize;

+	}

+    sectionEnd += step;

+    }

+slReverse(&itemList);

+

+/* Break up into sections of no more than items-per-slot size, and convert to packed format. */

+int sizeLeft = originalSectionSize;

+for (item = itemList; item != NULL; )

+    {

+    /* Figure out size of this section  */

+    int sectionSize = sizeLeft;

+    if (sectionSize > itemsPerSlot)

+        sectionSize = itemsPerSlot;

+    sizeLeft -= sectionSize;

+

+

+    /* Allocate and fill in section. */

+    struct bwgSection *section;

+    lmAllocVar(lm, section);

+    section->chrom = chrom;

+    section->start = sectionStart;

+    sectionStart += sectionSize * step;

+    section->end = sectionStart - step + span;

+    section->type = bwgTypeFixedStep;

+    section->itemStep = step;

+    section->itemSpan = span;

+    section->itemCount = sectionSize;

+