struct genoLay *genoLayNew(struct genoLayChrom *chromList,
	MgFont *font, int picWidth, int betweenChromHeight,
	int minLeftLabelWidth, int minRightLabelWidth,
	char *how)
/* Figure out layout.  For human and most mammals this will be
 * two columns with sex chromosomes on bottom.  This is complicated
 * by the platypus having a bunch of sex chromosomes. */
{
int margin = 3;
struct slRef *refList = NULL, *ref, *left, *right;
struct genoLayChrom *chrom;
struct genoLay *gl;
int autoCount, halfCount, bases, chromInLine;
int leftLabelWidth=0, rightLabelWidth=0, labelWidth;
int spaceWidth = mgFontCharWidth(font, ' ');
int extraLabelPadding = 0;
int autosomeOtherPixels=0, sexOtherPixels=0;
int autosomeBasesInLine=0;	/* Maximum bases in a line for autosome. */
int sexBasesInLine=0;		/* Bases in line for sex chromsome. */
double sexBasesPerPixel, autosomeBasesPerPixel, basesPerPixel;
int pos = margin;
int y = 0;
int fontHeight = mgFontLineHeight(font);
int chromHeight = fontHeight;
int lineHeight = chromHeight + betweenChromHeight;
boolean allOneLine = FALSE;

refList = refListFromSlList(chromList);

/* Allocate genoLay object and fill in simple fields. */
AllocVar(gl);
gl->chromList = chromList;
gl->chromHash = hashNew(0);
gl->font = font;
gl->picWidth = picWidth;
gl->margin = margin;
gl->spaceWidth = spaceWidth;
gl->lineHeight = lineHeight;
gl->betweenChromHeight = betweenChromHeight;
gl->betweenChromOffsetY = 0;
gl->chromHeight = chromHeight;
gl->chromOffsetY = lineHeight - chromHeight;


/* Save chromosomes in hash too, for easy access */
for (chrom = chromList; chrom != NULL; chrom = chrom->next)
    hashAdd(gl->chromHash, chrom->fullName, chrom);

if (sameString(how, genoLayOnePerLine))
    {
    gl->leftList = refList;
    }
else if (sameString(how, genoLayAllOneLine))
    {
    gl->bottomList = refList;
    allOneLine = TRUE;
    }
else
    {
    /* Put sex chromosomes on bottom, and rest on left. */
    separateSexChroms(refList, &refList, &gl->bottomList);
    autoCount = slCount(refList);
    gl->leftList = refList;

    /* If there are a lot of chromosomes, then move later
     * (and smaller) chromosomes to a new right column */
    if (autoCount > 12)
	{
	halfCount = (autoCount+1)/2;
	ref = slElementFromIx(refList, halfCount-1);
	gl->rightList = ref->next;
	ref->next = NULL;
	slReverse(&gl->rightList);
	}
    }

if (allOneLine)
    {
    unsigned long totalBases = 0, bStart=0, bEnd;
    int chromCount = 0, chromIx=0;
    for (ref = gl->bottomList; ref != NULL; ref = ref->next)
        {
	chrom = ref->val;
	totalBases += chrom->size;
	chromCount += 1;
	}
    int availablePixels = picWidth - minLeftLabelWidth - minRightLabelWidth
       - 2*margin - (chromCount-1);
    double basesPerPixel = (double)totalBases/availablePixels;
    gl->picHeight = 2*margin + lineHeight + fontHeight;
    for (ref = gl->bottomList; ref != NULL; ref = ref->next)
        {
	chrom = ref->val;
	bEnd = bStart + chrom->size;
	int pixStart = round(bStart / basesPerPixel);
	int pixEnd = round(bEnd / basesPerPixel);
	chrom->width = pixEnd - pixStart;
	chrom->height = lineHeight;
	chrom->x = pixStart + margin + chromIx + minLeftLabelWidth;
	chrom->y = 0;
//.........这里部分代码省略.........

static bioSeq *nextSeqFromMem(char **pText, boolean isDna, boolean doFilter)
/* Convert fa in memory to bioSeq.  Update *pText to point to next
 * record.  Returns NULL when no more sequences left. */
{
char *name = "";
char *s, *d;
struct dnaSeq *seq;
int size = 0;
char c;
char *filter = (isDna ? ntChars : aaChars);
char *text = *pText;
char *p = skipLeadingSpaces(text);
if (p == NULL)
    return NULL;
dnaUtilOpen();
if (*p == '>')
    {
    char *end;
    s = strchr(p, '\n');
    if (s != NULL) ++s;
    name = skipLeadingSpaces(p+1);
    end = skipToSpaces(name);
    if (end >= s || name >= s)
        errAbort("No name in line starting with '>'");
    if (end != NULL)
        *end = 0;
    }
else
    {
    s = p; 
    if (s == NULL || s[0] == 0)
        return NULL;
    }
name = cloneString(name);
    
d = text;
if (s != NULL)
    {
    for (;;)
	{
	c = *s;
	if (c == 0 || c == '>')
	    break;
	++s;
	if (!isalpha(c))
	    continue;
	if (doFilter)
	    {
	    if ((c = filter[(int)c]) == 0) 
		{
		if (isDna)
		    c = 'n';
		else
		    c = 'X';
		}
	    }
	d[size++] = c;
	}
    }
d[size] = 0;

/* Put sequence into our little sequence structure. */
AllocVar(seq);
seq->name = name;
seq->dna = text;
seq->size = size;
*pText = s;
return seq;
}

void reportAlt3Prime(struct altGraphX *ag, bool **em, int vs, int ve1, int ve2, 
		    int altBpStart, int altBpEnd, int startV, int endV, FILE *out)
/* Write out an altGraphX record for an alt3Prime splicing
event. Variable names are consistent with the rest of the program, but
can be misleading. Specifically vs = start of alt splicing, ve1 =
first end of alt splicing, etc. even though "vs" is really the end of
an exon. For an alt5Prime splice the edges are:

 Name       Vertexes         Class
 ------     ----------       -----
exon1:      startV->vs       constituative (0)
junction1:  vs->ve1          alternative (1)
junction2:  vs->ve2          alternative (2)
exon2:      ve1->e2        alternative (1)
exon3:      ve2->endV        constituative (0)
*/
{
struct altGraphX *agLoc = NULL;  /* Local altGraphX. */
struct evidence *ev = NULL, *evLoc = NULL;
int *vPos = ag->vPositions;
unsigned char *vT = ag->vTypes;
int *vPosLoc = NULL;    /* Vertex Positions. */
int *eStartsLoc = NULL; /* Edge Starts. */
int *eEndsLoc = NULL;   /* Edge ends. */
unsigned char *vTLoc = NULL;      /* Vertex Types. */
int *eTLoc = NULL;      /* Edge Types. */
int vCLoc = 0;
int eCLoc = 0;
int edgeIx = 0, vertexIx = 0;
int i =0;
struct dyString *dy = NULL;

if(out == NULL)
    return;
AllocVar(agLoc);
agLoc->tName = cloneString(ag->tName);
agLoc->name = cloneString(ag->name);
agLoc->tStart = vPos[startV];
agLoc->tEnd = vPos[endV];
agLoc->strand[0] = ag->strand[0];
agLoc->vertexCount = vCLoc = 6;
agLoc->edgeCount = eCLoc = 5;
agLoc->id = alt3Prime;
/* Allocate some arrays. */
AllocArray(vPosLoc, vCLoc);
AllocArray(eStartsLoc, eCLoc);
AllocArray(eEndsLoc, eCLoc);
AllocArray(vTLoc, vCLoc);
AllocArray(eTLoc, eCLoc);

/* Fill in the vertex positions. */
vertexIx = 0;
vPosLoc[vertexIx++] = vPos[startV]; /* 0 */
vPosLoc[vertexIx++] = vPos[vs];     /* 1 */
vPosLoc[vertexIx++] = vPos[ve1];    /* 2 */
vPosLoc[vertexIx++] = vPos[ve2];    /* 3 */
vPosLoc[vertexIx++] = vPos[ve2];    /* 4 */
vPosLoc[vertexIx++] = vPos[endV];   /* 5 */

/* Fill in the vertex types. */
vertexIx = 0;
vTLoc[vertexIx++] = vT[startV];
vTLoc[vertexIx++] = vT[vs];
vTLoc[vertexIx++] = vT[ve1];
vTLoc[vertexIx++] = vT[vs]; /* Faking a separate exon for the alt spliced portion. */
vTLoc[vertexIx++] = vT[ve2];
vTLoc[vertexIx++] = vT[endV];

edgeIx = 0;

/* Constitutive first exon. */
eStartsLoc[edgeIx] = 0;
eEndsLoc[edgeIx] = 1;
eTLoc[edgeIx] = 0;
ev = evidenceForEdge(ag, startV, vs);
evLoc = CloneVar(ev);
evLoc->mrnaIds = CloneArray(ev->mrnaIds, ev->evCount);
slAddHead(&agLoc->evidence, evLoc);
edgeIx++;

/* Alternative1 junction (shorter). */
eStartsLoc[edgeIx] = 1;
eEndsLoc[edgeIx] = 2;
eTLoc[edgeIx] = 1;
ev = evidenceForEdge(ag, vs, ve1);
evLoc = CloneVar(ev);
evLoc->mrnaIds = CloneArray(ev->mrnaIds, ev->evCount);
slAddHead(&agLoc->evidence, evLoc);
edgeIx++;

/* Alt2 junction (longer). */
eStartsLoc[edgeIx] = 1;
eEndsLoc[edgeIx] = 4;
eTLoc[edgeIx] = 2;
ev = evidenceForEdge(ag, vs, ve2);
evLoc = CloneVar(ev);
evLoc->mrnaIds = CloneArray(ev->mrnaIds, ev->evCount);
slAddHead(&agLoc->evidence, evLoc);
edgeIx++;

//.........这里部分代码省略.........

struct dgNodeRef *dgFindPath(struct diGraph *dg, struct dgNode *a, struct dgNode *b)
/* Find shortest path from a to b.  Return NULL if can't be found. */
{
struct dgNodeRef *refList  = NULL, *ref;
struct dgConnection *con;
struct dgNode *node, *nNode;
struct dlList *fifo;
struct dlNode *ffNode;
struct dgNode endNode;
int fifoSize = 1;

/* Do some quick and easy tests first to return if have no way out
 * of node A, or if B directly follows A. */
if (a->nextList == NULL)
    return NULL;
if (a == b)
    {
    AllocVar(ref);
    ref->node = a;
    return ref;
    }
if ((con = dgFindNodeInConList(a->nextList, b)) != NULL)
    {
    AllocVar(refList);
    refList->node = a;
    node = con->node;
    AllocVar(ref);
    ref->node = node;
    slAddTail(&refList, ref);
    return refList;
    }

/* Set up for breadth first traversal.  Will use a doubly linked
 * list as a fifo. */
for (node = dg->nodeList; node != NULL; node = node->next)
    node->tempEntry = NULL;
fifo = newDlList();
dlAddValTail(fifo, a);
a->tempEntry = &endNode;

while ((ffNode = dlPopHead(fifo)) != NULL)
    {
    --fifoSize;
    node = ffNode->val;
    freeMem(ffNode);
    for (con = node->nextList; con != NULL; con = con->next)
	{
	nNode = con->node;
	if (nNode->tempEntry == NULL)
	    {
	    nNode->tempEntry = node;
	    if (nNode == b)
		{
		while (nNode != &endNode && nNode != NULL)
		    {
		    AllocVar(ref);
		    ref->node = nNode;
		    slAddHead(&refList, ref);
		    nNode = nNode->tempEntry;
		    }
		break;
		}
	    else
		{
		dlAddValTail(fifo, nNode);
		++fifoSize;
		if (fifoSize > 100000)
		    errAbort("Internal error in dgFindPath");
		}
	    }
	}
    }
freeDlList(&fifo);
return refList;
}

struct sufa *sufaRead(char *fileName, boolean memoryMap)
/* Read in a sufa from a file.  Does this via memory mapping if you like,
 * which will be faster typically for about 100 reads, and slower for more
 * than that (_much_ slower for thousands of reads and more). */
{
/* Open file (low level), read in header, and check it. */
int fd = open(fileName, O_RDONLY);
if (fd < 0)
    errnoAbort("Can't open %s", fileName);
struct sufaFileHeader h;
if (read(fd, &h, sizeof(h)) < sizeof(h))
    errnoAbort("Couldn't read header of file %s", fileName);
if (h.magic != SUFA_MAGIC)
    errAbort("%s does not seem to be a sufa file.", fileName);
if (h.majorVersion > SUFA_MAJOR_VERSION)
    errAbort("%s is a newer, incompatible version of sufa format. "
             "This program works on version %d and below. "
	     "%s is version %d.",  fileName, SUFA_MAJOR_VERSION, fileName, h.majorVersion);

struct sufa *sufa;
verbose(2, "sufa file %s size %lld\n", fileName, h.size);

/* Get a pointer to data in memory, via memory map, or allocation and read. */
struct sufaFileHeader *header ;
if (memoryMap)
    {
#ifdef MACHTYPE_sparc
    header = (struct sufaFileHeader *)mmap(NULL, h.size, PROT_READ, MAP_SHARED, fd, 0);
#else
    header = mmap(NULL, h.size, PROT_READ, MAP_FILE|MAP_SHARED, fd, 0);
#endif
    if (header == (void*)(-1))
	errnoAbort("Couldn't mmap %s, sorry", fileName);
    }
else
    {
    header = needHugeMem(h.size);
    if (lseek(fd, 0, SEEK_SET) < 0)
	errnoAbort("Couldn't seek back to start of sufa file %s.  "
		   "Splix files must be random access files, not pipes and the like"
		   , fileName);
    if (read(fd, header, h.size) < h.size)
        errnoAbort("Couldn't read all of sufa file %s.", fileName);
    }

/* Allocate wrapper structure and fill it in. */
AllocVar(sufa);
sufa->header = header;
sufa->isMapped = memoryMap;

/* Make an array for easy access to chromosome names. */
int chromCount = header->chromCount;
char **chromNames = AllocArray(sufa->chromNames, chromCount);
char *s = pointerOffset(header, sizeof(*header) );
int i;
for (i=0; i<chromCount; ++i)
    {
    chromNames[i] = s;
    s += strlen(s)+1;
    }

/* Keep track of where we are in memmap. */
bits64 mapOffset = sizeof(*header) + header->chromNamesSize;

/* Point into chromSizes array. */
bits32 *chromSizes = sufa->chromSizes 
	= pointerOffset(header, mapOffset);
mapOffset += sizeof(bits32) * chromCount;

verbose(2, "total dna size %lld in %d chromosomes\n", (long long)header->dnaDiskSize, header->chromCount);
sufa->allDna = pointerOffset(header, mapOffset);
mapOffset += header->dnaDiskSize;

/* Calculate chromOffset array. */
bits32 offset = 0;
bits32 *chromOffsets = AllocArray(sufa->chromOffsets, chromCount);
for (i=0; i<chromCount; ++i)
    {
    chromOffsets[i] = offset;
    offset += chromSizes[i] + 1;
    verbose(2, "sufa contains %s,  %d bases, %d offset\n", 
    	sufa->chromNames[i], (int)sufa->chromSizes[i], (int)chromOffsets[i]);
    }

/* Finally point to the suffix array!. */
sufa->array = pointerOffset(header, mapOffset);
mapOffset += header->arraySize * sizeof(bits32);


assert(mapOffset == header->size);	/* Sanity check */
return sufa;
}

void doChainScore(char *chainIn, char *tNibDir, char *qNibDir, char *chainOut)
{
char qStrand = 0, tStrand = 0;
struct dnaSeq *qSeq = NULL, *tSeq = NULL;
char *qName = "",  *tName = "";
FILE *f = mustOpen(chainOut, "w");
struct chain *chainList = NULL, *chain;
struct chain *inputChains, *next;
FILE *details = NULL;
struct lineFile *lf = NULL;
struct dnaSeq *seq, *seqList = NULL;
struct hash *faHash = newHash(0);
struct hash *chainHash = newHash(0);
char comment[1024];
FILE *faF;
struct seqPair *spList = NULL, *sp;
struct dyString *dy = newDyString(512);
struct lineFile *chainsLf = lineFileOpen(chainIn, TRUE);

while ((chain = chainRead(chainsLf)) != NULL)
    {
    dyStringClear(dy);
    dyStringPrintf(dy, "%s%c%s", chain->qName, chain->qStrand, chain->tName);
    sp = hashFindVal(chainHash, dy->string);
    if (sp == NULL)
        {
	AllocVar(sp);
	slAddHead(&spList, sp);
	hashAddSaveName(chainHash, dy->string, sp, &sp->name);
	sp->qName = cloneString(chain->qName);
	sp->tName = cloneString(chain->tName);
	sp->qStrand = chain->qStrand;
	}
    slAddHead(&sp->chain, chain);
    }
slSort(&spList, seqPairCmp);
lineFileClose(&chainsLf);

if (optionExists("faQ"))
    {
    faF = mustOpen(qNibDir, "r");
    while ( faReadMixedNext(faF, TRUE, NULL, TRUE, NULL, &seq))
        {
        hashAdd(faHash, seq->name, seq);
        slAddHead(&seqList, seq);
        }
    fclose(faF);
    }
for (sp = spList; sp != NULL; sp = sp->next)
    {
    if (optionExists("faQ"))
        {
        assert (faHash != NULL);
        loadFaSeq(faHash, sp->qName, sp->qStrand, &qName, &qSeq, &qStrand);
        }
    else
        loadIfNewSeq(qNibDir, sp->qName, sp->qStrand, &qName, &qSeq, &qStrand);
    loadIfNewSeq(tNibDir, sp->tName, '+', &tName, &tSeq, &tStrand);
    scorePair(sp, qSeq, tSeq, &chainList, sp->chain);
    }


slSort(&chainList, chainCmpScore);
for (chain = chainList; chain != NULL; chain = chain->next)
    {
    assert(chain->qStart == chain->blockList->qStart 
	&& chain->tStart == chain->blockList->tStart);
    chainWrite(chain, f);
    }

carefulClose(&f);
}

static struct dnaSeq *dnaLoadNextFromStack(struct dnaLoad *dl)
/* Load next piece of DNA from stack of files.  Return NULL
 * when stack is empty. */
{
struct dnaLoadStack *dls;
struct dnaSeq *seq = NULL;
while ((dls = dl->stack) != NULL)
    {
    if (dls->twoBit)
        {
	if (dls->tbi != NULL)
	    {
	    seq = twoBitReadSeqFrag(dls->twoBit, dls->tbi->name, 0, 0);
	    dls->tbi = dls->tbi->next;
	    return seq;
	    }
	else
	    {
	    dl->stack = dls->next;
	    dnaLoadStackFree(&dls);
	    }
	}
    else if (dls->textIsFa)
        {
	DNA *dna;
	char *name;
	int size;
	if (faMixedSpeedReadNext(dls->textFile, &dna, &size, &name))
	    {
	    AllocVar(seq);
	    seq->dna = needLargeMem(size+1);
	    memcpy((void *)seq->dna, (void *)dna, size);
	    seq->dna[size] = 0;
	    seq->size = size;
	    seq->name = cloneString(name);
	    dl->curStart = 0;
	    dl->curEnd = size;
	    dl->curSize = size;
	    return seq;
	    }
	else
	    {
	    dl->stack = dls->next;
	    dnaLoadStackFree(&dls);
	    }
	}
    else	/* It's a file full of file names. */
        {
	char *line;
	if (lineFileNextReal(dls->textFile, &line))
	    {
	    line  = trimSpaces(line);
	    if ((seq = dnaLoadSingle(line, &dl->curStart, &dl->curEnd, &dl->curSize)) != NULL)
	         return seq;
	    else
	         {
		 struct dnaLoadStack *newDls;
		 newDls = dnaLoadStackNew(line);
		 slAddHead(&dl->stack, newDls);
		 }
	    }
	else
	    {
	    dl->stack = dls->next;
	    dnaLoadStackFree(&dls);
	    }
	}
    }
dl->finished = TRUE;
return NULL;
}

int main (int argc, char **argv)
{
  LineStream ls;
  Texta tokens = NULL;
  char *line;

  int hasQual = 0;
  int hasSeqs = 0;
  int start=1;
 
  ls = ls_createFromFile ("-");
  while (line = ls_nextLine (ls)) {
    // Put all the lines of the SAM header in comments
    if (line[0] == '@') {
      printf ("# %s\n", line);
      continue;
    }
    // Parse each SAM entry and store into array   
    tokens = textFieldtokP (line, "\t");
    if (arrayMax (tokens) < 11) {
      textDestroy( tokens );
      ls_destroy (ls);
      die ("Invalid SAM entry: %s", line);
    }
    SamEntry *currSamE = NULL;
    SamEntry *mateSamE = NULL;
    AllocVar(currSamE ); 

    int ret = generateSamEntry( tokens, currSamE, &hasSeqs, &hasQual );
    textDestroy( tokens );
    if ( ret==0 ) {
      if ( isPaired ( currSamE ) )
	ls_nextLine( ls ); // discarding next entry too (the mate)
      destroySamEntry( currSamE );
      freeMem( currSamE );
      continue;
    }   
    if ( isPaired( currSamE ) )   {
      int hasQual2, hasSeq2;
      AllocVar( mateSamE );
      Texta secondEnd = NULL;
      secondEnd = textFieldtok (ls_nextLine( ls ) , "\t");
      ret = generateSamEntry( secondEnd, mateSamE, &hasSeq2, &hasQual2 );
      textDestroy( secondEnd );
      if( ret == 0 ) {
	destroySamEntry( currSamE );
	destroySamEntry( mateSamE );
	freeMem( currSamE );
	freeMem( mateSamE );
	continue;
      }
      if (strcmp (currSamE->qname, mateSamE->qname) != 0) {
        die ("Please note that for paired-end data, sam2mrf requires the mate pairs to be on subsequent lines. You may want to sort the SAM file first.\nEx: sort -r file.sam | sam2mrf > file.mrf\n");
      }
    } 

    // Print MRF headers
    if( start ) {
      printf ("%s", MRF_COLUMN_NAME_BLOCKS);
      if (hasSeqs) printf("\t%s", MRF_COLUMN_NAME_SEQUENCE);
      if (hasQual) printf("\t%s", MRF_COLUMN_NAME_QUALITY_SCORES);
      printf ("\t%s\n", MRF_COLUMN_NAME_QUERY_ID);
      start=0;
    }
    
    // Print AlignmentBlocks   
    printMrfAlignBlocks (currSamE, R_FIRST);
    if( isPaired ( currSamE ) ) {  
      printf ("|");
      printMrfAlignBlocks (mateSamE, R_SECOND);
    }

    seq_init();
    // Print Sequence
    if (hasSeqs) {
      if (!currSamE->seq)
        die ("Entry missing sequence column\n");
      if( currSamE->flags & S_QUERY_STRAND )
	seq_reverseComplement( currSamE->seq, strlen(currSamE->seq));
      printf ("\t%s", currSamE->seq);
      if (mateSamE) {
        if (!mateSamE->seq)
          die ("Entry missing sequence column\n");
        if( mateSamE->flags & S_MATE_STRAND )
	  seq_reverseComplement( mateSamE->seq, strlen(mateSamE->seq));
	printf ("|%s", mateSamE->seq);
      }
    }
    // Print quality scores
    if (hasQual) {
      if (!currSamE->qual)
        die ("Entry missing quality scores column\n");
      printf ("\t%s", currSamE->qual);
      if (mateSamE) {
        if (!mateSamE->qual)
          die ("Entry missing quality scores column\n");
        printf ("|%s", mateSamE->qual);
      }
    }

//.........这里部分代码省略.........

void *esmStartHandler(struct xap *xp, char *name, char **atts)
/* Called by expat with start tag.  Does most of the parsing work. */
{
struct xapStack *st = xp->stack+1;
int depth = xp->stackDepth;
int i;

if (sameString(name, "Motifs"))
    {
    struct esmMotifs *obj;
    AllocVar(obj);
    for (i=0; atts[i] != NULL; i += 2)
        {
        char *name = atts[i], *val = atts[i+1];
        if  (sameString(name, "SeqFile"))
            obj->SeqFile = cloneString(val);
        }
    if (obj->SeqFile == NULL)
        xapError(xp, "missing SeqFile");
    return obj;
    }
else if (sameString(name, "Motif"))
    {
    struct esmMotif *obj;
    AllocVar(obj);
    for (i=0; atts[i] != NULL; i += 2)
        {
        char *name = atts[i], *val = atts[i+1];
        if  (sameString(name, "Consensus"))
            obj->Consensus = cloneString(val);
        else if (sameString(name, "Source"))
            obj->Source = cloneString(val);
        else if (sameString(name, "Name"))
            obj->Name = cloneString(val);
        else if (sameString(name, "Description"))
            obj->Description = cloneString(val);
        }
    if (obj->Consensus == NULL)
        xapError(xp, "missing Consensus");
    if (obj->Source == NULL)
        xapError(xp, "missing Source");
    if (obj->Name == NULL)
        xapError(xp, "missing Name");
    if (depth > 1)
        {
        if  (sameString(st->elName, "Motifs"))
            {
            struct esmMotifs *parent = st->object;
            slAddHead(&parent->esmMotif, obj);
            }
        }
    return obj;
    }
else if (sameString(name, "Weights"))
    {
    struct esmWeights *obj;
    AllocVar(obj);
    for (i=0; atts[i] != NULL; i += 2)
        {
        char *name = atts[i], *val = atts[i+1];
        if  (sameString(name, "ZeroWeight"))
            obj->ZeroWeight = atof(val);
        }
    if (depth > 1)
        {
        if  (sameString(st->elName, "Motif"))
            {
            struct esmMotif *parent = st->object;
            slAddHead(&parent->esmWeights, obj);
            }
        }
    return obj;
    }
else if (sameString(name, "Position"))
    {
    struct esmPosition *obj;
    AllocVar(obj);
    for (i=0; atts[i] != NULL; i += 2)
        {
        char *name = atts[i], *val = atts[i+1];
        if  (sameString(name, "Num"))
            obj->Num = atoi(val);
        else if (sameString(name, "Weights"))
            obj->Weights = cloneString(val);
        }
    if (obj->Weights == NULL)
        xapError(xp, "missing Weights");
    if (depth > 1)
        {
        if  (sameString(st->elName, "Weights"))
            {
            struct esmWeights *parent = st->object;
            slAddHead(&parent->esmPosition, obj);
            }
        }
    return obj;
    }
else
    {
    xapSkip(xp);
//.........这里部分代码省略.........

void agpVsMap(char *agpName, char *infoName, char *gifName)
/* agpVsMap - Plot clones in agp vs. map coordinates. */
{
struct mapPos *mapList, *mp;
struct agpFrag *agpList, *bp;
struct hash *cloneHash = newHash(14);
struct hashEl *hel;
struct cloneInfo *cloneList = NULL, *clone;
struct memGfx *mg = NULL;
int pixWidth = 600;
int pixHeight = 600;
int rulerHeight = 20;
int maxMapPos = 0, maxAgpPos = 0;
double scaleMap, scaleAgp;
Color orange, green;

mapList = readInfoFile(infoName);
agpList = readAgpFile(agpName);

for (mp = mapList; mp != NULL; mp = mp->next)
    {
    if (mp->phase > 0)
        {
	AllocVar(clone);
	hel = hashAddUnique(cloneHash, mp->cloneName, clone);
	clone->name = hel->name;
	clone->mp = mp;
	slAddHead(&cloneList, clone);
	if (mp->pos > maxMapPos) maxMapPos = mp->pos;
	}
    }
slReverse(&cloneList);

for (bp = agpList; bp != NULL; bp = bp->next)
    {
    if (bp->chromStart > maxAgpPos) maxAgpPos = bp->chromStart;
    }

/* Draw scatterplot on bitmap. */
mg = mgNew(pixWidth, pixHeight);
mgClearPixels(mg);
orange = mgFindColor(mg, 210, 150, 0);
green = mgFindColor(mg, 0, 200, 0);
mgDrawRuler(mg, 0, pixHeight-rulerHeight, rulerHeight, pixWidth, MG_BLACK,
       mgSmallFont(), 0, maxMapPos+1);
scaleMap = (double)pixWidth/(double)(maxMapPos+1.0);
scaleAgp = (double)(pixHeight)/(double)(maxAgpPos+1.0);
for (bp = agpList; bp != NULL; bp = bp->next)
    {
    char cloneName[128];
    fragToCloneName(bp->frag, cloneName);
    clone = hashFindVal(cloneHash, cloneName);
    if (clone == NULL)
        warn("%s is in %s but not %s", cloneName, 
	    agpName, infoName);
    else
	{
	int x = round(scaleMap*clone->mp->pos);
	int y = pixHeight - round(scaleAgp*bp->chromStart);
	int phase = clone->mp->phase;
	int back;
	if (phase <= 1) back = green;
	else if (phase == 2) back = orange;
	else back = MG_RED;
	drawPlus(mg, x, y, back);
	}
    }

mgSaveGif(mg, gifName);
}

struct fullExperiment *getFullExperimentList(struct sqlConnection *conn,
    struct edwExperiment *eeList, char *assembly, struct hash **retHash)
/* Given list of edwExperiments, return list of ones replicated with full file sets on
 * both replicates.  If optional retHash is non-NULL then return a hash full of same 
 * experiments keyed by experiment accession */
{
/* Build up a list of fullExperiments and a hash keyed by name. */
struct hash *hash = hashNew(14);
struct fullExperiment *fullList = NULL;
struct edwExperiment *ee;
for (ee = eeList; ee != NULL; ee = ee->next)
    {
    struct fullExperiment *full = hashFindVal(hash, ee->accession);
    if (full == NULL)
        {
	AllocVar(full);
	full->name = cloneString(ee->accession);
	full->exp = ee;
	slAddHead(&fullList, full);
	hashAdd(hash, full->name, full);
	}
    }
uglyf("Got %d in eeList, %d in fullList, %d in hash\n", slCount(eeList), slCount(fullList), hash->elCount);

/* Build up SQL query to efficiently fetch all good files and valid files from our experiment */
struct dyString *q = dyStringNew(16*1024);
sqlDyStringPrintf(q, "select edwValidFile.*,edwFile.*,eapOutput.* "
		  " from edwValidFile,edwFile,eapOutput "
                  " where edwValidFile.fileId = edwFile.id and edwFile.id = eapOutput.fileId "
		  " and edwFile.deprecated='' and edwFile.errorMessage='' "
		  " and edwValidFile.ucscDb != 'centro.hg19' "
		  " and edwValidFile.ucscDb like '%%%s' and edwValidFile.experiment in ("
		  , assembly);
for (ee = eeList; ee != NULL; ee = ee->next)
    {
    dyStringPrintf(q, "'%s'", ee->accession);
    if (ee->next != NULL)
        dyStringAppendC(q, ',');
    }
dyStringAppendC(q, ')');

/* Loop through this making up vFiles that ultimately are attached to replicates. */
int vCount = 0;
struct sqlResult *sr = sqlGetResult(conn, q->string);
char **row;
while ((row = sqlNextRow(sr)) != NULL)
    {
    ++vCount;
    struct edwValidFile *valid = edwValidFileLoad(row);
    fixOutputType(valid);
    struct edwFile *file = edwFileLoad(row + EDWVALIDFILE_NUM_COLS);
    struct eapOutput *eapOutput = eapOutputLoad(row + EDWVALIDFILE_NUM_COLS + EDWFILE_NUM_COLS);
    struct vFile *vf = vFileNew(file, valid, eapOutput);
    struct fullExperiment *full = hashMustFindVal(hash, valid->experiment);
    struct replicate *rep = findOrMakeReplicate(valid->replicate, &full->repList);
    char *format = valid->format;
    if (sameString(format, "bam"))
        slAddHead(&rep->bamList, vf);
    else if (sameString(format, "bigWig"))
        slAddHead(&rep->bigWigList, vf);
    else if (sameString(format, "narrowPeak") && !sameString(valid->outputType, "replicated_narrowPeak"))
        slAddHead(&rep->narrowList, vf);
    else if (sameString(format, "broadPeak") && !sameString(valid->outputType, "replicated_broadPeak"))
        slAddHead(&rep->broadList, vf);
    }
sqlFreeResult(&sr);
uglyf("Got %d vFiles\n", vCount);

dyStringFree(&q);

/* Free hash or return it, and return list. */
if (retHash == NULL)
    hashFree(&hash);
else
    *retHash = hash;
return fullList;
}

//.........这里部分代码省略.........
	hTvDropDownClass("ruler", rulerMode, FALSE, rulerMode ? "normalText" : "hiddenText");
	hPrintf("</TD>");
	hPrintf("<TD>");
	hPrintf("Chromosome position in bases.  (Clicks here zoom in 3x)");
	hPrintf("</TD>");
#ifdef PRIORITY_CHANGES_IN_CONFIG_UI
        if (withPriorityOverride)
            {
            hPrintf("<TD>");
            hPrintf("</TD>");
            hPrintf("<TD>");
            hPrintf("</TD>");
            }
#endif///def PRIORITY_CHANGES_IN_CONFIG_UI
	hPrintf("</TR>\n");
	}
    if (differentString(group->name, "user"))
        isFirstNotCtGroup = FALSE;
    /* Scan track list to determine which supertracks have visible member
     * tracks, and to insert a track in the list for the supertrack.
     * Sort tracks and supertracks together by priority */
    groupTrackListAddSuper(cart, group);

    if (!withPriorityOverride)
        {
        /* sort hierarchically by priority, considering supertracks */
        struct trackRef *refList = NULL, *ref;
        for (tr = group->trackList; tr != NULL; tr = tr->next)
            {
            struct track *track = tr->track;
            if (tdbIsSuperTrackChild(track->tdb))
                /* ignore supertrack member tracks till supertrack is found */
                continue;
            AllocVar(ref);
            ref->track = track;
            slAddTail(&refList, ref);
            if (tdbIsSuper(track->tdb))
                {
                struct trackRef *tr2;
                for (tr2 = group->trackList; tr2 != NULL; tr2 = tr2->next)
                    {
                    char *parent = tr2->track->tdb->parentName;
                    if (parent && sameString(parent, track->track))
                        {
                        AllocVar(ref);
                        ref->track = tr2->track;
                        slAddTail(&refList, ref);
                        }
                    }
                }
            }
        group->trackList = refList;
        }

    /* Loop through this group and display */
    int rowCount=1;
    for (tr = group->trackList; tr != NULL; tr = tr->next)
	{
	struct track *track = tr->track;
        struct trackDb *tdb = track->tdb;

	hPrintf("<TR %sid='%s-%d'>",(isOpen ? "" : "style='display: none'"),group->name, rowCount++);
	hPrintf("<TD NOWRAP>");
        if (tdbIsSuperTrackChild(tdb))
            /* indent members of a supertrack */
            hPrintf("&nbsp;&nbsp;&nbsp;&nbsp;");

static struct grp *makeGroupList(char *db, struct trackDb *trackList, struct grp **pHubGrpList,
                                 boolean allTablesOk)
/* Get list of groups that actually have something in them. */
{
struct grp *groupsAll, *groupList = NULL, *group;
struct hash *groupsInTrackList = newHash(0);
struct hash *groupsInDatabase = newHash(0);
struct trackDb *track;

/* Stream through track list building up hash of active groups. */
for (track = trackList; track != NULL; track = track->next)
    {
    if (!hashLookup(groupsInTrackList,track->grp))
        hashAdd(groupsInTrackList, track->grp, NULL);
    }

/* Scan through group table, putting in ones where we have data. */
groupsAll = hLoadGrps(db);
for (group = slPopHead(&groupsAll); group != NULL; group = slPopHead(&groupsAll))
    {
    if (hashLookup(groupsInTrackList, group->name))
	{
	slAddTail(&groupList, group);
	hashAdd(groupsInDatabase, group->name, group);
	}
    else
        grpFree(&group);
    }

/* if we have custom tracks, we want to add the track hubs
 * after that group */
struct grp *addAfter = NULL;
if ((groupList != NULL) && sameString(groupList->name, "user"))
    addAfter = groupList;

/* Add in groups from hubs. */
for (group = slPopHead(pHubGrpList); group != NULL; group = slPopHead(pHubGrpList))
    {
    // if the group isn't represented in any track, don't add it to list
    if (!hashLookup(groupsInTrackList,group->name))
	continue;
    /* check to see if we're inserting hubs rather than
     * adding them to the front of the list */
    struct grp *newGrp = grpDup(group);
    if (addAfter != NULL)
	{
	newGrp->next = addAfter->next;
	addAfter->next = newGrp;
	}
    else
	slAddHead(&groupList, newGrp);
    hashAdd(groupsInDatabase, newGrp->name, newGrp);
    }

/* Do some error checking for tracks with group names that are
 * not in database.  Just warn about them. */
if (!trackHubDatabase(db))
    for (track = trackList; track != NULL; track = track->next)
    {
    if (!hashLookup(groupsInDatabase, track->grp))
         warn("Track %s has group %s, which isn't in grp table",
	 	track->table, track->grp);
    }

/* Create dummy group for all tracks. */
AllocVar(group);
group->name = cloneString("allTracks");
group->label = cloneString("All Tracks");
slAddTail(&groupList, group);

/* Create another dummy group for all tables. */
if (allTablesOk)
    {
    AllocVar(group);
    group->name = cloneString("allTables");
    group->label = cloneString("All Tables");
    slAddTail(&groupList, group);
    }

hashFree(&groupsInTrackList);
hashFree(&groupsInDatabase);
return groupList;
}

struct improbRunInfo * analyseOneMotifRun(char *runName, char *seqDir,
        char *motifDir, int controlCount, char *controls[])
/* Bundle up data on one improbizer run and associated control runs. */
{
    char fileName[512];
    char motifName[256];
    int seqCount, baseCount;
    struct improbRunInfo *iriList = NULL, *iri;
    struct lineFile *lf = NULL;
    struct motif motif;
    int motifIx = 0;
    int i;
    float acc, best, mean, x;

    printf("%s\n", runName);

    /* Count bases in sequences - this will be used in each iri. */
    sprintf(fileName, "%s/%s.fa", seqDir, runName);
    countSeq(fileName, &seqCount, &baseCount);

    /* Allocate iri and read the main run. */
    sprintf(fileName, "%s/%s", motifDir, runName);
    lf = lineFileOpen(fileName, TRUE);
    while (readMotif(lf, &motif))
    {
        AllocVar(iri);
        slAddTail(&iriList, iri);
        ++motifIx;
        snprintf(motifName, sizeof(motifName), "%s.%d", runName, motifIx);
        iri->name = cloneString(motifName);
        iri->seqCount = seqCount;
        iri->runScore = motif.score;
        iri->runPos = motif.pos;
        iri->runPosSd = motif.posSd;
        iri->columnCount = motif.size;
        iri->consensus = cloneString(motif.consensus);
        iri->aProb = CloneArray(motif.profile[0], motif.size);
        iri->cProb = CloneArray(motif.profile[1], motif.size);
        iri->gProb = CloneArray(motif.profile[2], motif.size);
        iri->tProb = CloneArray(motif.profile[3], motif.size);
        iri->controlCount = controlCount;
        AllocArray(iri->controlScores, controlCount);
    }
    lineFileClose(&lf);

    /* Read the control runs. */
    for (i=0; i<controlCount; ++i)
    {
        sprintf(fileName, "%s/%s", controls[i], runName);
        lf = lineFileOpen(fileName, TRUE);
        for (iri = iriList; iri != NULL; iri = iri->next)
        {
            if (!readMotif(lf, &motif))
                errAbort("%s doesn't contain the expected number of motifs", lf->fileName);
            iri->controlScores[i] = motif.score;
        }
        lineFileClose(&lf);
    }

    /* Calculate best and mean on control runs. */
    for (iri = iriList; iri != NULL; iri = iri->next)
    {
        acc = best = 0;
        for (i=0; i<controlCount; ++i)
        {
            x = iri->controlScores[i];
            acc += x;
            if (x > best)
                best = x;
        }
        iri->bestControlScore = best;
        iri->meanControlScore = acc/controlCount;
    }

    /* Calculate standard deviation of control runs. */
    for (iri = iriList; iri != NULL; iri = iri->next)
    {
        acc = 0;
        mean = iri->meanControlScore;
        for (i=0; i<controlCount; ++i)
        {
            x = iri->controlScores[i] - mean;
            acc += x*x;
        }
        if (controlCount > 1)
            acc /= controlCount;
        iri->sdControlScore = sqrt(acc);
    }

    return iriList;
}

struct bbiChromUsage *bbiChromUsageFromBedFile(struct lineFile *lf, 
	struct hash *chromSizesHash, int *retMinDiff, double *retAveSize, bits64 *retBedCount)
/* Go through bed file and collect chromosomes and statistics. */
{
char *row[3];
struct hash *uniqHash = hashNew(0);
struct bbiChromUsage *usage = NULL, *usageList = NULL;
int lastStart = -1;
bits32 id = 0;
bits64 totalBases = 0, bedCount = 0;
int minDiff = BIGNUM;

lineFileRemoveInitialCustomTrackLines(lf);

for (;;)
    {
    int rowSize = lineFileChopNext(lf, row, ArraySize(row));
    if (rowSize == 0)
        break;
    lineFileExpectWords(lf, 3, rowSize);
    char *chrom = row[0];
    int start = lineFileNeedNum(lf, row, 1);
    int end = lineFileNeedNum(lf, row, 2);
    if (start > end)
        {
	    errAbort("end (%d) before start (%d) line %d of %s",
	    	end, start, lf->lineIx, lf->fileName);
	}
    ++bedCount;
    totalBases += (end - start);
    if (usage == NULL || differentString(usage->name, chrom))
        {
	if (hashLookup(uniqHash, chrom))
	    {
	    errAbort("%s is not sorted at line %d.  Please use \"sort -k1,1 -k2,2n\" or bedSort and try again.",
	    	lf->fileName, lf->lineIx);
	    }
	hashAdd(uniqHash, chrom, NULL);
	struct hashEl *chromHashEl = hashLookup(chromSizesHash, chrom);
	if (chromHashEl == NULL)
	    errAbort("%s is not found in chromosome sizes file", chrom);
	int chromSize = ptToInt(chromHashEl->val);
	AllocVar(usage);
	usage->name = cloneString(chrom);
	usage->id = id++;
	usage->size = chromSize;
	slAddHead(&usageList, usage);
	lastStart = -1;
	}
    if (end > usage->size)
        errAbort("End coordinate %d bigger than %s size of %d line %d of %s", end, usage->name, usage->size, lf->lineIx, lf->fileName);
    usage->itemCount += 1;
    if (lastStart >= 0)
        {
	int diff = start - lastStart;
	if (diff < minDiff)
	    {
	    if (diff < 0)
		errAbort("%s is not sorted at line %d.  Please use \"sort -k1,1 -k2,2n\" or bedSort and try again.",
		    lf->fileName, lf->lineIx);
	    minDiff = diff;
	    }
	}
    lastStart = start;
    }
slReverse(&usageList);
*retMinDiff = minDiff;
*retAveSize = (double)totalBases/bedCount;
*retBedCount = bedCount;
freeHash(&uniqHash);
return usageList;
}

static struct blastBlock *nextBlock(struct blastFile *bf, struct blastQuery *bq,
                                    struct blastGappedAli *bga, boolean *skipRet)
/* Read in next blast block.  Return NULL at EOF or end of gapped
 * alignment. If an unparsable block is found, set skipRet to TRUE and return
 * NULL. */
{
struct blastBlock *bb;
char *line;
char *words[16];
int wordCount;
char *parts[3];
int partCount;
static struct dyString *qString = NULL, *tString = NULL;

verbose(TRACE_LEVEL,  "blastFileNextBlock\n");
*skipRet = FALSE;

/* Seek until get something like:
 *   Score = 8770 bits (4424), Expect = 0.0
 * or something that looks like we're done with this gapped
 * alignment. */
for (;;)
    {
    if (!nextBlockLine(bf, bq, &line))
	return NULL;
    if (startsWith(" Score", line))
	break;
    }
AllocVar(bb);
bb->gappedAli = bga;
wordCount = chopLine(line, words);
if (wordCount < 8 || !sameWord("Score", words[0]) 
    || !isdigit(words[2][0]) || !(isdigit(words[7][0]) || words[7][0] == 'e')
    || !startsWith("Expect", words[5]))
    {
    bfError(bf, "Expecting something like:\n"
             "Score = 8770 bits (4424), Expect = 0.0");
    }
bb->bitScore = atof(words[2]);
bb->eVal = evalToDouble(words[7]);

/* Process something like:
 *   Identities = 8320/9618 (86%), Gaps = 3/9618 (0%)
 *             or
 *   Identities = 8320/9618 (86%)
 *             or
 *   Identities = 10/19 (52%), Positives = 15/19 (78%), Frame = +2
 *     (wu-tblastn)
 *             or
 *   Identities = 256/400 (64%), Positives = 306/400 (76%)
 *   Frame = +1 / -2
 *     (tblastn)
 *
 *   Identities = 1317/10108 (13%), Positives = 2779/10108 (27%), Gaps = 1040/10108
 *   (10%)
 *      - wrap on long lines
 *
 * Handle weird cases where the is only a `Score' line, with no `Identities'
 * lines by skipping the alignment; they seem line small, junky alignments.
 */
line = bfNeedNextLine(bf);
wordCount = chopLine(line, words);
if (wordCount < 3 || !sameWord("Identities", words[0]))
    {
    if (wordCount > 1 || sameWord("Score", words[0]))
        {
        /* ugly hack to skip block with no identities */
        *skipRet = TRUE;