db8efa2db09244e76a2e9d41daf607ab310e69b9 hartera Thu Jan 22 13:25:58 2015 -0800 Remove extra comments. Changed printf statements to verbose. diff --git src/hg/pslCDnaGenomeMatch/pslCDnaGenomeMatch.c src/hg/pslCDnaGenomeMatch/pslCDnaGenomeMatch.c index a99877f..ac9ca71 100644 --- src/hg/pslCDnaGenomeMatch/pslCDnaGenomeMatch.c +++ src/hg/pslCDnaGenomeMatch/pslCDnaGenomeMatch.c @@ -34,31 +34,30 @@ static char na[3] = "NA"; struct axtScoreScheme *ss = NULL; /* blastz scoring matrix */ struct hash *snpHash = NULL, *mrnaHash = NULL, *faHash = NULL, *tHash = NULL, *species1Hash = NULL, *species2Hash = NULL; int maxGap = 100; int minDiff = MINDIFF; int notAlignPenalty = NOTALIGNEDFACTOR; int aliCount = 0; /* number of alignments read */ int mrnaCount = 0; /* number of mrna read */ int filterCount = 0; /* number of mrna where best hit can be determined */ int outputCount = 0; /* number of alignments written */ int verbosity = 1; int lociCounter = 0; bool passthru = FALSE; bool computeSS = FALSE; struct dlList *fileCache = NULL; -//struct hash *nibHash = NULL; FILE *outFile = NULL; /* output tab sep file */ FILE *bedFile = NULL; /* output bed file */ FILE *scoreFile = NULL; /* output score file */ struct twoBitFile *cdnaSeqFile = NULL; struct twoBitFile *genomeSeqFile = NULL; char *species1 = NULL; char *species2 = NULL; char *nibdir1 = NULL; char *nibdir2 = NULL; char *mrna1 = NULL; char *mrna2 = NULL; char *bedOut = NULL; char *scoreOut = NULL; char *snpFile = NULL; /* snp tab file (browser format)*/ int histogram[256][256]; /* histogram of counts for suff statistics, index is a,c,g,t,-,. */ @@ -181,31 +180,30 @@ tSeq = el->tSeq; qSeq = el->qSeq; freeDnaSeq(&tSeq); freeDnaSeq(&qSeq); freez(pEl); } void alignFreeList(struct alignment **pList) /* Free a list of dynamically allocated alignment's */ { struct alignment *el, *next; for (el = *pList; el != NULL; el = next) { next = el->next; - //pslFreeList(&el->psl); alignFree(&el); } *pList = NULL; } bool purine(char a) { if (toupper(a) == 'A' || toupper(a) == 'G') return(TRUE); else return(FALSE); } bool transversion(char a, char b) { @@ -230,41 +228,35 @@ { struct dnaSeq *seq; int size = end - start; int seqSize = twoBitSeqSize(tbf, seqName); assert(size >= 0); if (strand == '-') { verbose(6,"before start %d end %d size %d %s\n", start, end, size, seqName); reverseIntRange(&start, &end, seqSize); verbose(6,"after start %d end %d size %d %s\n", start, end, size, seqName); assert(start >= 0); seq = twoBitReadSeqFrag(tbf, seqName, start, end); - //seq = nibInfoLoadSeq(nib, start, size); -// seq = nibTwoCacheSeqPart(ntc, seqName, -// start, size, &retFullSeqSize); reverseComplement(seq->dna, seq->size); } else { seq = twoBitReadSeqFrag(tbf, seqName, start, end); - //seq = nibInfoLoadSeq(nib, start, size); -// seq = nibTwoCacheSeqPart(ntc, seqName, -// start, size, &retFullSeqSize); } return seq; } void addLoci(struct loci **lociList, struct psl *psl) /* add a new loci from psl, if not already in list */ { struct loci *loci = NULL; bool found = FALSE; if (!found) { AllocVar(loci); loci->chrom = cloneString(psl->tName); loci->chromStart = psl->tStart; @@ -544,31 +536,30 @@ char *chrom, int chromStart, int mrnaLoc, char strand, struct loci *lociList, struct snp125 *snpList) /* add new mismatch to list, figure out which loci we are mapped to */ { struct loci *l = NULL; struct misMatch *misMatch = NULL; for (l = lociList ; l != NULL ; l = l->next) { int index = getLoci(lociList, chrom, chromStart); if (index == l->index) misMatch = newMisMatch(name, offset, genomeBase, mrnaBase, chrom, chromStart, mrnaLoc, strand, l->index, lociList, snpList); else misMatch = newMisMatch(name, offset, '-', mrnaBase, na, -1, mrnaLoc, '.' , l->index, lociList, snpList); - //printMisMatch(&misMatch); slAddHead(misMatchList, misMatch); } } struct mrnaMisMatch *tabulateMisMatches(struct misMatch *mm, int seqCount, struct loci *lociList) /* get all mismatches for an mRNA and group all mismatches by position in mrna */ { struct misMatch *mme = NULL; struct mrnaMisMatch *mrnaMisMatch = NULL, *misMatchList = NULL; int prevLoc = -99; if (mm == NULL) return NULL; for (mme = mm ; mme != NULL ; mme = mme->next) { @@ -590,39 +581,36 @@ AllocArray(mrnaMisMatch->chroms, seqCount); AllocArray(mrnaMisMatch->loci, seqCount); for (j = 0 ; j<seqCount; j++) { /* default to empty alignment */ mrnaMisMatch->bases[j] = '.'; mrnaMisMatch->tStarts[j] = 0; mrnaMisMatch->chroms[j] = cloneString(na); mrnaMisMatch->loci[j] = j /* NOVALUE mme->loci*/ ; } mrnaMisMatch->bases[seqCount] = '\0'; mrnaMisMatch->snpCount = 0; AllocArray(mrnaMisMatch->snps, seqCount); if (mme->snpCount > 0) { -// for (j = 0 ; j<mrnaMisMatch->snpCount; j++) - { mrnaMisMatch->snps[i] = cloneString(mme->snps[0]); mrnaMisMatch->snpCount++; verbose(4,"mrnaloc %d mmsnp [%d] loci [%d] %s\n", mrnaMisMatch->mrnaLoc,i, mme->loci, mrnaMisMatch->snps[i]); } } - } mrnaMisMatch->chroms[i] = cloneString(mme->chrom); mrnaMisMatch->tStarts[i] = mme->chromStart; mrnaMisMatch->bases[i] = mme->genomeBase; mrnaMisMatch->strands[i] = mme->strand; mrnaMisMatch->loci[i] = mme->loci; assert(mrnaMisMatch->loci[i] != NOVALUE); prevLoc = mme->mrnaLoc; } slReverse(&misMatchList); return misMatchList; } int misMatchCmpMrnaLoc(const void *va, const void *vb) /* Compare to sort based on mrnaLoc. */ { @@ -839,38 +827,30 @@ if (scoreFile) fprintf(scoreFile, "##>> %s score %d %s:%d-%d [%d] mismatch %d good %d neither %d indel %d \ gaps %d snps %d total %d diff %d maxScore %d maxHITS %d 2nd best %d diff %d\n", name, score, psl->tName , chromStart, chromEnd, z, missCount[z], goodCount[z], neither[z], indel, missCount[z]+ goodCount[z]+ neither[z]+ indel, gapCount[z], snpCount[z], seqCount - slCount(lociList), maxScore, maxHits, nextBestScore, diff); if (psl->strand[0] == '+' && psl->strand[1] == '-') pslRc(psl); if (psl->strand[1] == '+' || psl->strand[1] == '-') psl->strand[1] = '\0'; pslTabOut(psl, outFile); outputCount++; filterCount++; qNameCount++; - /* - freez(&missCount); - freez(&goodCount); - freez(&gapCount); - freez(&neither); - freez(&snpCount); - return TRUE; - */ } else { verbose(2, "%s NO score %d diff %d spread %d max %d next %d maxHits %d index %d pos %s:%d diff %d < %d\n", name, score , diff, spread, maxScore, nextBestScore, maxHits, z, psl->tName , psl->tStart, diff, minDiff); } } freez(&missCount); freez(&goodCount); freez(&gapCount); freez(&neither); freez(&snpCount); if (qNameCount > 0) return TRUE; @@ -934,120 +914,109 @@ { unsigned int i, j; for (i = 0 ; i < 256 ; i++) for (j = 0 ; j < 256 ; j++) { histogram[i][j] = 0; histoNorm[i][j] = 0.0; } } void computeSuffStats(struct misMatch **misMatchList, struct alignment *align, struct loci *lociList) /* compute sufficient statistics from each alignment */ { int blockIx = 0; -//int i, j, sum = 0; struct psl *psl = align->psl; -//struct nibInfo *tNib = nibInfoFromCache(nibHash, align->nibDir, psl->tName); -//nibTwoCache static struct dnaSeq *tSeq = NULL; int tStart = psl->tStart; int tEnd = psl->tEnd; -//int misMatchCount = 0; char genomeStrand = psl->strand[1] == '-' ? '-' : '+'; if (genomeStrand == '-') reverseIntRange(&tStart, &tEnd, psl->tSize); for (blockIx=0; blockIx < psl->blockCount; ++blockIx) /* for each alignment block get sequence for both strands */ { struct snp125 *snp = NULL, *snpList = NULL; if (hashFindVal(cdnaSeqFile->hash, psl->qName) == NULL) { - printf("skipping %s not found \n",psl->qName); + verbose(5, "Skipping %s not found \n",psl->qName); return; } int size = twoBitSeqSize(cdnaSeqFile, psl->qName); int i = 0; int ts = psl->tStarts[blockIx]; int te = psl->tStarts[blockIx]+(psl->blockSizes[blockIx]); int qe = psl->qStarts[blockIx]+(psl->blockSizes[blockIx]); struct dnaSeq *qSeq = NULL; if (qe > size) { - printf("skipping %s %d > %d \n",psl->qName, qe, size ); + verbose(5, "Skipping %s %d > %d \n",psl->qName, qe, size ); return; } qSeq = twoBitReadSeqFrag(cdnaSeqFile, psl->qName, psl->qStarts[blockIx], psl->qStarts[blockIx]+(psl->blockSizes[blockIx])); if (genomeStrand == '-') reverseIntRange(&ts, &te, psl->tSize); verbose(5," xyz %s t %s:%d-%d q %d-%d %s strand %c\n", psl->qName, psl->tName, ts, te, psl->qStarts[blockIx], psl->qStarts[blockIx]+psl->blockSizes[blockIx], psl->strand, genomeStrand); tSeq = twoBitGetSeqForStrand(align->twoBitSeqFile, psl->tName, psl->tStarts[blockIx], psl->tStarts[blockIx]+(psl->blockSizes[blockIx]), genomeStrand); verbose(6,"tSeq %s len %d %d\n",tSeq->dna, tSeq->size, (int)strlen(tSeq->dna)); verbose(6,"qSeq %s\n",qSeq->dna); assert(psl->strand[0] == '+'); /* count match and mismatches in each alignment block */ for (i = 0 ; i < tSeq->size; i++) { char t = toupper(tSeq->dna[i]); char q = toupper(qSeq->dna[i]); int genomeStart = ts+i; - //int mrnaLoc = psl->qStarts[blockIx]+i; if (q == 'n' || q == 'N') - printf("N in sequence %s\n",psl->qName); + verbose(5, "N in sequence %s\n",psl->qName); updateCount(t,q); if (genomeStrand == '-') genomeStart = te-i-1; /* get overlapping snps */ snpList = getSnpList(psl->tName, genomeStart, genomeStart+1, genomeStrand) ; for (snp = snpList ; snp != NULL ; snp = snp->next) { int offset = (snp->chromStart)-ts; char snpBase = ' '; assert(sameString(psl->tName, snp->chrom)); if (genomeStrand == '-') offset = (tSeq->size - offset)-1; snpBase = toupper(tSeq->dna[offset]); if (snpBase != t) { verbose(4,"snp mismatch genome %c snpbase %c mrna %c ts %d snp %d valid %s\n", t, snpBase, q, ts, snp->chromStart, snp->valid); -// assert(snpBase == t); } else verbose(4," [%d] snp match %s %s %s:%d %s %s offset %d gs %c ts %d genomic%c %c valid %s\n", getLoci(lociList, psl->tName, psl->tStart), snp->name, psl->qName, snp->chrom, snp->chromStart, snp->observed, snp->strand, offset, genomeStrand, ts, genomeStrand, t, snp->valid); } } - /* add indel for portions of mrna that do not align at all */ -// for (i = 0 ; i < psl->qStarts[0] ; i++) -// { -// char q = toupper(qSeq->dna[i]); -// } freeDnaSeq(&tSeq); freeDnaSeq(&qSeq); } } void addOtherAlignments( struct alignment *alignList, struct hash *speciesHash, char *name, struct twoBitFile *twoBitFile, char *mrnaPath) { /* add alignemnts from other species to the list */ struct hashEl *speciesList = NULL , *el = NULL; struct alignment *align = NULL; if (speciesHash != NULL) speciesList = hashLookup(speciesHash,name); for (el = speciesList ; el != NULL ; el = el->next) { @@ -1056,31 +1025,30 @@ align->psl = psl; align->twoBitSeqFile = twoBitFile; align->mrnaPath = mrnaPath; slAddHead(&alignList, align); } slReverse(&alignList); } void fillinMatches(struct misMatch **misMatchList, struct alignment *align, struct loci *lociList) /* for each mismatch , lookup bases in other loci */ { int blockIx = 0; static struct dnaSeq *tSeq = NULL; struct psl *psl = align->psl; - //struct nibInfo *tNib = nibInfoFromCache(nibHash, align->nibDir, psl->tName); int tStart = psl->tStart; int tEnd = psl->tEnd; char genomeStrand = psl->strand[1] == '-' ? '-' : '+'; int index = getLoci(lociList, psl->tName, psl->tStart); if (genomeStrand == '-') reverseIntRange(&tStart, &tEnd, psl->tSize); for (blockIx=0; blockIx < psl->blockCount; ++blockIx) /* for each alignment block get sequence for both strands */ { struct dnaSeq *qSeq = twoBitReadSeqFrag(cdnaSeqFile, psl->qName, psl->qStarts[blockIx], psl->qStarts[blockIx]+(psl->blockSizes[blockIx])); int ts = psl->tStarts[blockIx]; int te = psl->tStarts[blockIx]+(psl->blockSizes[blockIx]); @@ -1089,94 +1057,81 @@ struct misMatch *mm = NULL; if (genomeStrand == '-') reverseIntRange(&ts, &te, psl->tSize); tSeq = twoBitGetSeqForStrand(align->twoBitSeqFile, psl->tName, psl->tStarts[blockIx], psl->tStarts[blockIx]+(psl->blockSizes[blockIx]), genomeStrand); for (mm = *misMatchList ; mm != NULL ; mm = mm->next) { /* i = offset within the block */ int i = mm->mrnaLoc - qs; int genomeStart = ts+i; if (genomeStrand == '-') genomeStart = te-i-1; - /* if genomic chr not filled in yet (only mismatches are filled in ), and we are within the block, lookup genomic loc and base*/ + /* if genomic chr not filled in yet (only mismatches are filled in ), and we are within the block, + * lookup genomic loc and base*/ if (index == mm->loci) { if (sameString(mm->chrom , na) && mm->mrnaLoc >= qs && mm->mrnaLoc < qe && i >= 0) { char t = toupper(tSeq->dna[i]); char q = toupper(qSeq->dna[i]); if (q=='\0') q='-'; if (t=='\0') t='-'; if (i > 0 || i <= qe-qs) verbose(6," i %d qs %d qe %d\n",i,qs,qe); assert (i > 0 || i <= qe-qs) ; -// verbose(6, "WAS %s not %s\n",mm->chrom, psl->tName); mm->chrom = cloneString(psl->tName); mm->chromStart = genomeStart; mm->genomeBase = t; mm->strand = genomeStrand; if (q!=(mm->mrnaBase) && q!='-') { verbose(2, "mismatch %s %s q %c != mmBase %c offset %d\n",psl->qName, psl->tName,q,(mm->mrnaBase), i ); assert(q==(mm->mrnaBase)); } verbose(5," fM() FILL %s %c/%c t %s:%d q %d i %d qs %d qe %d %s %c loci %d BLK %d ", psl->qName, t, q, psl->tName, genomeStart, mm->mrnaLoc, i, qs, qe ,psl->strand, genomeStrand, mm->loci, blockIx); if (t==q) verbose(5,"GOOD\n"); else verbose(5,"MISMATCH\n"); } /* if we are outside the block, assume an indel */ else if (mm->mrnaLoc >= qs && mm->mrnaLoc < qe) { mm->genomeBase = '-'; verbose(5," fM() INDEL %s t %s:%d mrnaLoc %d i %d not in %d-%d %s loci %d BLK %d\n", psl->qName, psl->tName, psl->tStart, mm->mrnaLoc, i, qs, qe , psl->strand, mm->loci, blockIx); } } } - /* - snpList = getSnpList(psl->tName, ts, te, genomeStrand) ; - for (snp = snpList ; snp != NULL ; snp = snp->next) - { - int offset = (snp->chromStart)-ts; - if (genomeStrand == '-') - offset = (tSeq->size - offset)-1; - verbose(4," snp %s %s %s:%d offset %d %s %s gs %c ts %d genomic%c valid %s\n", - snp->name, psl->qName, snp->chrom, snp->chromStart, offset, snp->observed, - snp->strand, genomeStrand, ts, genomeStrand, snp->valid); - } - */ freeDnaSeq(&tSeq); freeDnaSeq(&qSeq); } } void buildMisMatches(struct misMatch **misMatchList, struct alignment *align, struct loci *lociList) /* build list of mismatches for each alignment (loci) */ { int blockIx = 0; struct psl *psl = align->psl; -//struct nibInfo *tNib = nibInfoFromCache(nibHash, align->nibDir, psl->tName); static struct dnaSeq *tSeq = NULL; int tStart = psl->tStart; int tEnd = psl->tEnd; int misMatchCount = 0; int transitionCount = 0; int transversionCount = 0; char genomeStrand = psl->strand[1] == '-' ? '-' : '+'; int matchCount = 0; for (blockIx=0; blockIx < psl->blockCount; ++blockIx) /* for each alignment block get sequence for both strands */ { struct snp125 *snp = NULL, *snpList = NULL; struct dnaSeq *qSeq = twoBitReadSeqFrag(cdnaSeqFile, psl->qName, psl->qStarts[blockIx], psl->qStarts[blockIx]+(psl->blockSizes[blockIx])); int i = 0; @@ -1219,84 +1174,53 @@ /* get overlapping snps */ snpList = getSnpList(psl->tName, genomeStart, genomeStart+1, genomeStrand) ; for (snp = snpList ; snp != NULL ; snp = snp->next) { int offset = (snp->chromStart)-ts; char snpBase = ' '; assert(sameString(psl->tName, snp->chrom)); if (genomeStrand == '-') offset = (tSeq->size - offset)-1; snpBase = toupper(tSeq->dna[offset]); if (snpBase != t) { verbose(4,"snp mismatch genome %c snpbase %c mrna %c ts %d snp %d valid %s\n", t, snpBase, q, ts, snp->chromStart, snp->valid); -// assert(snpBase == t); } else verbose(4," [%d] snp match %s %s %s:%d %s %s offset %d gs %c ts %d genomic%c %c valid %s\n", getLoci(lociList, psl->tName, psl->tStart), snp->name, psl->qName, snp->chrom, snp->chromStart, snp->observed, snp->strand, offset, genomeStrand, ts, genomeStrand, t, snp->valid); } verbose(5," add mismatch genome %c mrna %c ts %d gn %d offset %d mrnaLoc %d\n", t, q, ts, genomeStart, i, mrnaLoc); newMisMatches(misMatchList, psl->qName, i ,t, q, psl->tName, genomeStart, mrnaLoc, genomeStrand, lociList, snpList); } else { matchCount++; } } -// /* add indel for portions of mrna that do not align at all */ -// for (i = 0 ; i < psl->qStarts[0] ; i++) -// { -// char q = toupper(qSeq->dna[i]); -// -// verbose(4, "mismatch gap %c offset %d \n", -// q, i); -// newMisMatches(misMatchList, psl->qName, i ,'-', q, psl->tName, -// -1, i, genomeStrand, lociList, snpList); -// } - - -// if (misMatchCount == 0) -// newMisMatches(misMatchList, psl->qName, -1 ,'.', '.', psl->tName, -// 0, -1, genomeStrand, lociList, snpList); printMisMatch(misMatchList); - /* - snpList = getSnpList(psl->tName, ts, te, genomeStrand) ; - for (snp = snpList ; snp != NULL ; snp = snp->next) - { - int offset = (snp->chromStart)-ts; - if (genomeStrand == '-') - offset = (tSeq->size - offset)-1; - verbose(4," snp %s %s %s:%d offset %d %s %s gs %c ts %d genomic%c valid %s\n", - snp->name, psl->qName, snp->chrom, snp->chromStart, offset, snp->observed, - snp->strand, genomeStrand, ts, genomeStrand, snp->valid); - } - */ freeDnaSeq(&tSeq); freeDnaSeq(&qSeq); } verbose(4,"buildMismatches final score %s:%d-%d %s %s #match %d misMatch %d \n", psl->tName, tStart, tEnd, psl->strand, psl->qName, matchCount, misMatchCount); -// fprintf(scoreFile,"%s\t%d\t%d\t%s\t%d\n", -// psl->tName, psl->tStart, psl->tEnd, psl->qName, -// misMatchCount); } void doOneMrna(char *name, struct alignment *alignList) /* build list of mismatches based on all alignments (defined by locilist) of mrna to genome */ { struct misMatch *misMatchList = NULL; struct alignment *align = NULL; int seqCount = slCount(alignList); struct loci *lociList = NULL; mrnaCount++; lociCounter = 0; if (seqCount == 1) { @@ -1311,32 +1235,30 @@ filterCount++; return; } /* one loci for each place the mrna alignments to the genome */ for (align = alignList ; align != NULL ; align= align->next) { struct psl *psl = align->psl; addLoci(&lociList, psl); verbose(5,"add loci for %s:%d-%d\n", psl->tName, psl->tStart, psl->tEnd); } verbose(5, "name %s alignList %d loci %d \n", name, slCount(alignList), slCount(lociList)); -//addOtherAlignments(&alignList, species1Hash, name, nibDir1, mrna1); -//addOtherAlignments(&alignList, species2Hash, name, nibDir2, mrna2); if (alignList != NULL) verbose(5,"buildMisMatch %s aList count %d \n",name, slCount(alignList)); for (align = alignList; align != NULL ; align = align->next) { /* for each alignment, build mismatch list */ verbose(5,"buildMisMatch %s tName %s:%d-%d\n", name, align->psl->tName, align->psl->tStart, align->psl->tEnd); if (computeSS) computeSuffStats(&misMatchList, align, lociList); else buildMisMatches(&misMatchList, align, lociList); } slReverse(&misMatchList); if (!computeSS ) { @@ -1417,31 +1339,30 @@ align->psl = newPsl; align->db = "db"; align->twoBitSeqFile = genomeSeqFile; align->mrnaPath = NULL; rm = slRemoveEl(&pslList, newPsl); if (rm) { slAddHead(&subList, align); verbose(5, "add %s %s\n",align->psl->qName, align->psl->tName); } } verbose(5,"last doOneMrna\n"); slSort(&subList, pslCmpMatch); doOneMrna(lastName, subList); alignFreeList(&subList); -//pslFreeList(&pslList); verbose(1,"Wrote %d alignments out of %d \n", outputCount, aliCount); verbose(1,"%d out of %d mRNAs passed strict criteria (pass thru not counted)\n", filterCount, mrnaCount); if (computeSS) { fprintf(outFile, "%8d %8d %8d %8d %8d %8d %8d %8d %8d %8d %8d %8d %8d %8d %8d %8d \n", getHistogram('a','a'), getHistogram('a','c'), getHistogram('a','g'), getHistogram('a','t'), getHistogram('c','a'), getHistogram('c','c'), getHistogram('c','g'), getHistogram('c','t'), getHistogram('g','a'), getHistogram('g','c'), getHistogram('g','g'), getHistogram('g','t'), getHistogram('t','a'), getHistogram('t','c'), getHistogram('t','g'), getHistogram('t','t') ); for (i = 0 ; i < 128 ; i++) for (j = 0 ; j < 128 ; j++) if (getHistogram(i,j) > 0) fprintf(outFile, "%d %d %c %c = %d \n",i, j, (char)i, (char)j, getHistogram(i,j)) ; normalizeSS(); @@ -1469,32 +1390,30 @@ getHistoNorm('n','a')+ getHistoNorm('n','c')+ getHistoNorm('n','g')+ getHistoNorm('n','t') ); for (i = 0 ; i < 128 ; i++) for (j = 0 ; j < 128 ; j++) if (getHistogram(i,j) > 0) fprintf(outFile, "%c %c = %f \n",(char)i, (char)j, getHistoNorm(i,j)) ; } } int main(int argc, char *argv[]) /* Process command line. */ { optionInit(&argc, argv, optionSpecs); if (argc != 6) usage(); -//if (nibHash == NULL) -// nibHash = hashNew(0); lociCounter = 0; fileCache = newDlList(); verbosity = optionInt("verbose", verbosity); verboseSetLogFile("stdout"); verboseSetLevel(verbosity); ss = axtScoreSchemeDefault(); cdnaSeqFile = twoBitOpen(argv[3]); genomeSeqFile = twoBitOpen(argv[4]); outFile = fopen(argv[5],"w"); initSS(); minDiff = optionInt("minDiff", MINDIFF); notAlignPenalty = optionInt("notAlignPenalty", NOTALIGNEDFACTOR); snpFile = optionVal("snp", NULL); if (snpFile != NULL)